Re: [ublas] LQ factorization with lapack bindings

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.

Best,

Karl

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
> http://en.wikipedia.org/wiki/Qr_decomposition
> and also
> http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse
>
>   
>> + 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]
> http://lists.boost.org/mailman/listinfo.cgi/ublas
>