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