minimum norm LS - different solns from different routines

[tadmor]

Hi-

I am using LAPACK to obtain the minimum norm least squares solution to a rank-deficient double precision (real) system of equations, see:

http://www.netlib.org/lapack/lug/node27.html#llsq

I used two of the routines provided for this in LAPACK: DGELSS and DGELSD, and obtained different answers from each. I also tried using the SVD routine in Numerical Recipes and obtained a third answer (different from the above two). All three solutions satisfy the original system of equations.

There are of course and infinite number of solutions to a rank-deficient system of equations, but there is only one that has a minimum norm. Of the three solutions I obtained the Numerical Recipes code gave the result with the lowest norm. The norms of the two LAPACK solutions were larger and different from each other. I don't understand why the LAPACK routines do not give me the same result as the Numerical Recipes code.

What am I missing?

[Julien Langou]

James Demmel wrote

[-- Mon Mar 7 00:02:09 EST 2005 --]


When you solve a rank deficient least squares problem with the SVD,
it works as follows:
(1) compute the SVD, and set the singular values smaller than a user supplied
threshold to zero;
(2) compute the minimum norm solution of the resulting exactly rank deficient problem

So depending on the threshold, you might get very different answers.
Furthermore, if the threshold is chosen very close to a true singular value,
roundoff might cause the computed singular value to be slightly greater
than the threshold (and so kept unchanged) or slightly less than the threshold
(and so zeroed out).

So the threshold should ideally be chosen as follows:
(1) at least a small multiple of macheps*norm(A), where macheps = machine epsilon
is the roundoff error, because there is this much uncertainty in any
computed singular value
(2) at least equal to the inherent uncertainty in the input data, which
only the user can know. For example, if the data is only good to
6 decimal digits, the threshold should be at least 1e-6 * norm(A).
(3) in a "gap" in the spectrum, i.e. far from any true singular value,
so roundoff won't change the computed rank. This is hard to
tell until the routine has actually computed the singular values.

So differences in roundoff and thresholds could cause different
implementations of the routine to get very different answers.
Typically the routines report rank and the singular values themselves,
so you can check to see what happened.

Jim Demmel

[tadmor]

Hi Jim,

Thank you for your timely and detailed response.

For both routines, I set rcond=-1, which sets the threshold to machine precision, so I assumed both would behave the same. However following your message, I took a closer look at the results and saw that one of the singular values was slightly larger than machine precision in one routine and lower than machine precision in the other. This difference was enough to cause the discrepancy. I changed rtol to 1d-10 (which hits a gap in the spectrum) and both routines give the same result as the Numerical Recipes code now. Thanks!

regards,