Re: Linear least squares using SVD

From: Jerry W. Lewis (post_a_reply_at_no_e-mail.com)
Date: 06/28/04 

Date: Mon, 28 Jun 2004 12:00:42 GMT

See
   http://lib.stat.cmu.edu/apstat/274
   http://lib.stat.cmu.edu/apstat/75
for a Fortran implementation

Jerry

Dave Dodson wrote:

> spasmous@yahoo.com (spasmous) wrote in message news:<2f066762.0406251603.643ca28a@posting.google.com>...
>
>>kmcgarvey@key.net (Ken McGarvey) wrote in message news:<d4e32a1.0406251001.4358e4c4@posting.google.com>...
>>
>>>I have a linear system of 256 unknowns and up to 300,000 data samples.
>>> I would like to do a least squares fit and compute the weights of 256
>>>basis functions.
>>>
>>2. Yes. NR is saying SVD is the most accurate method, albeit slow.
>>3. Yes. Transforming Ax=b to At.Ax=At.b causes a loss of precision
>>(you have basically squared A). Doing the SVD on A rather than At.A
>>avoids this.
>>
>
> If you don't have enough memory to store that 300,000 by 256 matrix,
> but you have access to the matrix one row at a time, you can do a
> least squares fit via a QR factorization with only a 256 by 256 array
> and a few vectors by using Givens rotations. This avoids the loss of
> precision that you would incur by forming the normal equations.
>
> The algorithm goes like this:
>
> 1. Put the first 256 rows of the 300,000 by 256 matrix in the array,
> say A.
> 2. For i = 1 to 255 do
> For j = i+1 to 255 do
> Use a Givens rotation with row i and row j to annihilate
> A(i,j).
> [as you do this, you convert the data in A into upper
> triangular form]
> 3. For i = 257 to 300000 do
> Put row i into a one-dimensional array, say r.
> For j = 1 to 256 do
> Use a Givens rotation with row j and r to annihilate r(j).
> 4. When you have processed all 300,000 rows, A contains the upper
> triangular matrix R of the QR factorization. If you have been
> processing the right hand side vector with the Givens rotations at the
> same time, say, in a 257th column of A, then you can back substitute
> to solve Rx=b, where b is that 257th column of A. The vector x will be
> the solution to your least squares problem.
>
> Dave
>