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" 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