LAPACK's equivalent to Matlab's lscov?

Open discussion regarding features, bugs, issues, vendors, etc.

LAPACK's equivalent to Matlab's lscov?

Postby dfabregat » Mon Apr 18, 2011 11:16 am

Hello everyone,

I would like to know if there is an equivalent call to Matlab's lscov in LAPACK, more concretely lscov in this form:

x = lscov(A,b,V), where V is an m-by-m real symmetric positive definite matrix, returns the generalized least squares solution to the linear system A*x = b with covariance matrix proportional to V, that is, x minimizes (b - A*x)'*inv(V)*(b - A*x).

I've been looking at LAPACK's documentation and this forum but I can't find an answer :S
Can anyone help me on this?

Thanks,
Diego
dfabregat
 
Posts: 2
Joined: Mon Apr 18, 2011 11:04 am

Re: LAPACK's equivalent to Matlab's lscov?

Postby Julien Langou » Wed Apr 27, 2011 11:23 pm

Hi,

( If V is diagonal, stop reading and say it !!! )

[
I initially wrote " is not diagonal " that was a typo. Corrected. I meant "is diagonal".
What I tried to say is that often V is diagonal SPD, but of course the problem is pretty much trivial.
If V is not diagonal, please read below.
]

What about
1) Cholesky factorization of V. (Compute L such that V=LL^T ). Use DPOTRF.
2) Then scale the system. So overwrite A with ( L \ A ) (use DTRSM), overwrite B with ( L \ B ) (use DTRSM or DTRSV).
3) call DGELS on ( A, B )

Another way to do this would be to
1) call LAPACK DPOTRF on the matrix V to obtain its Cholesky factor, (upper version,)
2) zero out the lower part of the Cholesky factor
2) call LAPACK DGGGLM with B, the Cholesky factor of V.
I think this works. I never tried. I think this is more stable.

Julien.
Last edited by Julien Langou on Sat Apr 30, 2011 5:36 pm, edited 1 time in total.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: LAPACK's equivalent to Matlab's lscov?

Postby dfabregat » Sat Apr 30, 2011 4:28 pm

Hi Julien,

V is not diagonal, it is a full symmetric positive definite matrix. Did you mean ( If V is not SPD, stop reading and say it !!! )?

Anyway, both approaches you mentioned work! Thanks a lot! :) Only one comment: for the second one (DGGGLM based) to work, I had to factor V as L * L' instead of U' * U.

Thanks again,
Diego
dfabregat
 
Posts: 2
Joined: Mon Apr 18, 2011 11:04 am

Re: LAPACK's equivalent to Matlab's lscov?

Postby Julien Langou » Sat Apr 30, 2011 5:36 pm

Hi Diego,
Great that all this works. Have fun with it.
Julien.
PS: I corrected my previous message, there was indeed a typo. J.
Julien Langou
 
Posts: 727
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 1 guest