From: israel@math.ubc.ca (Robert Israel)
Newsgroups: sci.math.symbolic
Subject: Re: Help on a QR decomposition. Anyone?
Date: 19 Apr 1995 19:44:11 GMT

In article <3n2rq2$1rg@epsilon.qmw.ac.uk>, "E.O.D.C.D.COVAS" <AH4137@QMWCC2.qmw.ac.uk> writes:
|> My problem is: first, how do I do the QR decompostion, I mean I could
|> eventualy do it a a algebraic system (well, in fact no cause I use
|> Maple instead of Mathematica and it has no QR there), but I must
|> do it in a algebraic form. It should be simple, I think.
|> 
|> What I know is that Q should be orthogonal and R upper triangular, but
|> what algorithm should I use, to construct Q and R? Every book talks 
|> about it but does not show the algorith! Is it similar to a SVD?

The QR factorization is basically the same as the Gram-Schmidt process.  It's 
curious that Maple has a "GramSchmidt" function but not a QR.  Even more
curious that "GramSchmidt" doesn't always output the vectors in the right
order, even if you give it a list rather than a set (which makes it just
about useless).

 So here's a QR procedure.  Note that you must first load in the "linalg" 
package.  The result of QR is a structure of the form Q &* R.  If you save 
it as F, you can get Q as op(1,F) and R as op(2,F).  No type checking is
done here, and you will get a "division by 0" error if the input matrix
does not have full column rank.


QR := proc(A)
       local QC, R, n, a, qp;
       n:= coldim(A);
       R:= matrix(n,n,0);
       for j from 1 to n do
         a:= col(A,j); qp:= a;
         for k from 1 to j-1 do 
           R[k,j]:= dotprod(QC[k],a);
           qp:= evalm(qp-"*QC[k]);
           od;
         R[j,j]:= norm(qp,2);
         QC[j]:= evalm(qp/");
         od;
       augment(seq(QC[j],j=1..n))&* eval(R);
       end;

For example:

> with(linalg):
> A:= matrix([[1,0,1,1],[1,1,0,-1],[1,0,-1,1],[1,1,0,1]]):
> QR(A);

         [                  1/2              ]
         [ 1/2  -1/2   1/2 2           0     ]    [ 2  1    0     1  ]
         [                                   ]    [                  ]
         [                               1/2 ]    [ 0  1    0    -1  ]
         [ 1/2   1/2       0      - 1/2 2    ]    [                  ]
         [                                   ] &* [        1/2       ]
         [                   1/2             ]    [ 0  0  2       0  ]
         [ 1/2  -1/2  - 1/2 2          0     ]    [                  ]
         [                                   ]    [              1/2 ]
         [                              1/2  ]    [ 0  0    0   2    ]
         [ 1/2   1/2       0       1/2 2     ]





 

-- 
Robert Israel                            israel@math.ubc.ca
Department of Mathematics             
University of British Columbia
Vancouver, BC, Canada V6T 1Y4