Page 1 of 1

Question about parameters for least squares

PostPosted: Tue Jul 11, 2006 3:00 pm
by jwwalker
I need to solve some linear least squares problems. The matrix is not necessarily full rank, so I guess I need to use xGELSY or xGELSD.

For xGELSY, I don't know what to pass for the JPVT and RCOND input parameters.

For xGELSD, I don't know what to pass for RCOND or what dimension to use for IWORK. How come the function can compute the size of WORK for you, but not IWORK?

I see that
NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) +1)
but what is SMLSIZ?

PostPosted: Tue Jul 11, 2006 6:37 pm
by jwwalker
By looking at the DGELSY example at, I can see that I can initialize JPVT to all zeros. I'm still not sure about RCOND. The example uses RCOND = 0.01 and comes up with the solution:

0.6344 0.9699 -1.4402 3.3678 3.3992

If I use RCOND = 0.01, or 0.001, I get about that answer. But if I use RCOND = 0.0001, I get

-0.799781 -3.288069 -7.475135 4.939312 0.767767

Now, I wouldn't expect to always get exactly the same solution, but since this algorithm is supposed to be finding a minimum-norm solution, I'd expect the norms to be fairly close... and they're not.

PostPosted: Wed Jul 12, 2006 3:55 am
by sven
RCOND is used to determine the numerical rank of the matrix A, and should be chosen to reflect the accuracy of your data.

Consider the example where

A = ( 1 0 ), b = ( 1 ).
( 0 eps ) ( 1 )

If eps is small, but treated as nonzero, so the rank of A is taken to be 2, then the solution is

x = 1, y = 1/eps.

On the other hand, if we treat eps as zero, so that we regard the rank of A as 1, then the solution is

x = 1, y = 0.

So you see that we can get quite a difference in solution.

When I run your example, I found that the value returned in RANK changed from 4 to 5 when RCOND = 0.0001.

I hope that helps,

Sven Hammarling.

PostPosted: Wed Jul 12, 2006 3:50 pm
by jwwalker
Thanks, that helps.