QR function

[Ted Harding]

( Re Message From: John W. Eaton )
> 
> On 14-Jul-1997, Moritz Harteneck <address@hidden> wrote:
> 
> | I am trying to make a QR decomposition of a fairly large matrix 
> | (ca. 1584x252 real elements). However, I am not interested in the 
> | rotation matrix q and I want to avoid it due to memory reasons. Does 
> | anybody know how to do this?
> 
> With Octave 2.0.9, you can get R by running qr(x) and asking for only
> one output:
> 
>   octave:1> help qr
>   qr is a builtin function
> 
>     ...
> 
>   qr (X) alone returns the output of the LAPACK routine dgeqrf, such
>   that R = triu (qr (X))
> 
> jwe

That's what "help qr" says for 2.0.8 too; but you get Q instead.

However, all is not entirely lost. Since with [Q,R] = qr(A) we have
A = Q*R (Q orthogonal, R upper-triangular), we also have A'*A = R'*R.
Hence if you do

    [Qs,Rs] = qr(sqrtm(A'*A))

you should find that Rs is the upper part (the upper-triangular square
block) of R, except that some rows of Rs will have opposite signs to the
corresponding rows of R.

You should still be able to use Rs for the same purposes as R, however.
In fact, the general QR decomposition is not unique, but the one with the
diagonal elements of R positive and decreasing is unique; so if you go
down Rs changing the sign of each row so that Rs(i,i)>0 then you will get
the same as if you did this to R:

   function s=sgn(x) s=x./abs(x); endfunction
   [r,c] = size(A); [Qs,Rs,ps] = qr(sqrtm(A'*A));
   for i=1:c, Ss = Rs(i,:)*sgn(R(i,i)); endfor
   R = [Ss ; zeros(r-c, c)]

(if you really need a full-sized "R").

This way you waste space on Qs, but in your case it is 252x252 instead of
1584x1584 (which looks like 508,032 bytes versus 20,072,448 bytes, to me).

You can always do "clear Q" afterwards.

Maybe this is helpful?

Ted.                                    (address@hidden)