Question about parameters for least squares

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

Question about parameters for least squares

Postby jwwalker » Tue Jul 11, 2006 3:00 pm

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
LIWORK >= 3 * MINMN * NLVL + 11 * MINMN
where
MINMN = MIN( M,N)
and
NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) +1)
but what is SMLSIZ?
jwwalker
 
Posts: 6
Joined: Mon Jul 10, 2006 6:18 pm
Location: San Diego, CA

Postby jwwalker » Tue Jul 11, 2006 6:37 pm

By looking at the DGELSY example at http://www.nag.com/lapack-ex/node48.html, 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.
jwwalker
 
Posts: 6
Joined: Mon Jul 10, 2006 6:18 pm
Location: San Diego, CA

Postby sven » Wed Jul 12, 2006 3:55 am

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.
sven
 
Posts: 146
Joined: Wed Dec 22, 2004 4:28 am

Postby jwwalker » Wed Jul 12, 2006 3:50 pm

Thanks, that helps.
jwwalker
 
Posts: 6
Joined: Mon Jul 10, 2006 6:18 pm
Location: San Diego, CA


Return to User Discussion

Who is online

Users browsing this forum: Google [Bot] and 10 guests

cron