Boost logo

Ublas :

From: Karl Meerbergen (karl.meerbergen_at_[hidden])
Date: 2007-12-23 03:47:07

Alternatively, you can use the SVD (gesvd binding). This decomposes
the m x n matrix A matrix into
A = U S V^* where S is a diagonal m x n matrix, U is m x m unitary and
V is n x n unitary.
The diagonal elements of S are the singular values, in descending order.
If there are k nonzero singular values, you can write A as
A = U(:,1:k) * S(1:k,1:k) * V(:,1:k)^*

Solving A x = b in a least squares fashion is then done as
x = V(:,1:k) * inv(S(1:k,1:k)) * ( U(:1,k)^* * b )

More expensive than QR but more reliable.



Jens Seidel wrote:
> On Sat, Dec 22, 2007 at 09:32:07PM +0200, Sorkin Dima wrote:
>> I need to solve an undetermined system of equations
>> (minimum norm solution), several times, with the same
>> matrix. As far as I know Lapack does this via LQ
>> factorization.
>> In bindings, only QR factorization currently exists.
>> 1) Does someone have bindings for LQ ?
>> 2) Is there a simple way to solve undetermined system
>> using QR factorization
> Yep. The least square solution can be obtained as:
> x = R^(-1)*Q^T*y
> Please note that this decomposition often requires maximal column rank
> of A but allows also generalisations. If you don't have a maximal
> column rank you can simply dropping all dependent columns of A and
> obtain the QR decomposition for the reduced system. In this case you
> have to extent x by zeros for dropped columns.
> That's an easy algorithm to obtain a least square solution (there may
> exist multiple if the rank is not maximal). See
> and also
>> + some tricks with transposes,
> I think QR and LQ are indeed very similar and should be transformed into
> each other by using a transposed matrix (but I'm not sure, I use always
> QR).
>> using the existing bindings' functions ?
>> The matrix A can be pretty big and I would prefer not to
>> create copy of its transpose explicitly -
>> is it unavoidable if I use QR ?
> Please note that some QR decompositions (especially if Gram Schmidt
> orthogonalisation is used) are very unstable and could cause trouble with
> small systems as 5x5 (at least if you expect a high accuracy)!!!!
> Jens
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]