# FILE anton3QR.mpl

# A procedure to find the QR decomposition of a matrix using  GramSchmidt
with(linalg):
qrfac:= proc(A,Q)
    local n,Qcols,R,i;
    n:= coldim(A);
    Qcols:= GramSchmidt( [col(A,1..n)] );
    Q:= transpose( matrix(Qcols) );
#   Normalise Q by dividing each  col(Q,i)  by its norm.
    for i from 1 to n do
        Q:= mulcol( Q,i,1/norm(col(Q,i),2) )
    od;
#   Since  Q &* R = A  and  transpose(Q) &* Q = Id,
    R:= evalm(transpose(Q)&*A);
end:

# AR7 Section 6.3 Example 8 p323
# Find the QR decomposition of the matrix
A:= matrix([ [1,0,0],
             [1,1,0],
             [1,1,1] ]):
R:= qrfac( A,'Q'):                                 print ( Q,R);
equal( A, evalm(Q&*R) );

# The rest of the items in this file are from Rod James
# accompaniment to a different linear algebra text

# Hill Section 4.5 Example 5, p285
# Use QR decomposition to find the line  y = mx + b  of best fit for the points
#   (2,3), (3,4), (4,4), (7,7)
#
# We want the least-squares solution for   A c = y   where c = [b,m] and
A:= matrix([ [1,2], [1,3], [1,4], [1,7] ]):
y:= vector( [3, 4, 4, 7] ):
R:= qrfac(A,'Q'):            # See the function  qrfac  given above
c:= linsolve( R, evalm(transpose(Q)&*y) ):
#
# The line of best fit is
y = c[2]*x + c[1];
quit

