Application Center - Maplesoft
  Submit your work
Application Center Applications Application of the Modified Gram-Schmidt Algorithm Preview

App Preview:

Application of the Modified Gram-Schmidt Algorithm

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



Application of the Modified Gram-Schmidt Algorithm in Maple and how this can be applied to a least squares curve fitting problem.

by Douglas Edan Lewit, currently enrolled in Real Analysis and Numerical Linear Algebra at Illinois Institute of Technology in Chicago, IL, Fall Semester of 2013.

 

I'm currently enrolled in a graduate level numerical linear algebra class at IIT in Chicago, IL.  I often use Maple for various assignments and computing projects in this class.  Maple has two excellent packages for manipulating matrices--the older linalg package and the newer LinearAlgebra package, the latter being preferable when working on large matrices containing floating point terms rather than symbolic terms.  Maple's LinearAlgebra package contains a special command called QRDecomposition.  The 'QRDecomposition' command successfully factors any matrix into the product of two matrices, the first one, usually called the Q-matrix, having orthonormal column vectors, and the second matrix, usually called the R-matrix, being an upper triangular non-singular matrix.  In fitting curves to data one often encounters overdetermined systems of linear equations.  (An overdetermined system is where the number of linear equations is greater than the number of unknowns in the system.)  Solving such a system can be challenging.  The QR decomposition method offers one convenient approach because it breaks down the coefficient matrix into two matrices, the Q matrix, which has the convenient property that its inverse is the same as its transpose--assuming of course that the Q matrix is square--and the R matrix, which is easily inverted because of its upper triangular shape.  Using the Q and R matrices, we thus can avoid the problems of taking the inverse of A^(T)A,another approach to solving this problem.  A^T*A may not be an easily inverted matrix, and even if A^T*A is invertible, this matrix usually has a very large condition number, meaning that we may not necessarily be able to trust the accuracy of any numerical results found using this matrix.  (The condition number of a matrix is usually defined as the square root of the ratio of the largest singular value to the smallest singular value.  If a non-singular matrix has a very large condition number, then the matrix is said to be ill-conditioned or almost singular.)  The QR decomposition method provides a very powerful alternative to working directly with the A^T*A matrix.

 

Maple's QRDecomposition command basically utilizes one of two routines for generating the Q and R matrices.  If the matrix contains only integers and/or symbolic expressions, then Maple performs a QR decomposition using the Classical Gram-Schmidt algorithm.  If however, the matrix contains a mixture of integers and floating point decimals or only floating point decimals, then Maple carries out the QR decomposition of the matrix using Householder transformations.  My approach below uses a third alternative, the Modified Gram-Schmidt algorithm, which I read about in Chapter 8 of the textbook, NUMERICAL LINEAR ALGEBRA, by Lloyd N. Trefethen and David Bau III.

Computer problem 8.2 (using Maple).

> 

restart;

> 

mgs:=proc(A::Matrix)  # implement the Modified Gram-Schmidt algorithm in Maple where A is a matrix.
local m,n;  # these variables will store the dimensions of the matrix A.
local i,j,k,total;  # counter variables for the loops in the program.
global V,Q;  # these variables will store the column vectors of A and the orthonormalized vectors of A, respectively.
global R;  # an upper triangular matrix.
total:=time();  # initialize the variable that will calculate the program's execution time.
m,n:=LinearAlgebra:-Dimensions(A);
if m < n then error "The number of rows must be greater than or equal to the number of columns to use this algorithm."
else
for j to n do
V[j]:=convert(LinearAlgebra:-SubMatrix(A,1..m,j..j),Vector[column]);
end do;
unassign('Q','R');
R:=Matrix(n,n,shape=triangular[upper]);
R[1,1]:=LinearAlgebra:-VectorNorm(V[1],2);
Q[1]:=V[1]/R[1,1];
for j from 2 to n do
for i to j do
if i<>j then
R[i,j]:=LinearAlgebra:-DotProduct(Q[i],V[j]);
V[j]:=V[j]-R[i,j]*Q[i];
else
R[i,j]:=LinearAlgebra:-VectorNorm(V[j],2);
Q[j]:=V[j]/R[i,j];
end if;
end do;
end do;
Q:=[seq(entries(Q)[k][1],k=1..n)];
Q:=convert(Q,Matrix);
total:=time()-total;
printf("The QR factorization of the matrix based on the MGS algorithm required %0.8f CPU seconds.",total);
end if;
end proc;

proc (A::Matrix) local m, n, i, j, k, total; global V, Q, R; total := time(); m, n := LinearAlgebra:-Dimensions(A); if m < n then error "The number of rows must be greater than or equal to the number of columns to use this algorithm." else for j to n do V[j] := convert(LinearAlgebra:-SubMatrix(A, 1 .. m, j .. j), Vector[column]) end do; unassign('Q', 'R'); R := Matrix(n, n, shape = triangular[upper]); R[1, 1] := LinearAlgebra:-VectorNorm(V[1], 2); Q[1] := V[1]/R[1, 1]; for j from 2 to n do for i to j do if i <> j then R[i, j] := LinearAlgebra:-DotProduct(Q[i], V[j]); V[j] := -Q[i]*R[i, j]+V[j] else R[i, j] := LinearAlgebra:-VectorNorm(V[j], 2); Q[j] := V[j]/R[i, j] end if end do end do; Q := [seq(entries(Q)[k][1], k = 1 .. n)]; Q := convert(Q, Matrix); total := time()-total; printf("The QR factorization of the matrix based on the MGS algorithm required %0.8f CPU seconds.", total) end if end proc

(1)
> 

showstat(mgs);  # This output for a program is more traditional because the program lines are numbered in order.


mgs := proc(A::Matrix)
local m, n, i, j, k, total;
global V, Q, R;
   1   total := time();
   2   m, n := LinearAlgebra:-Dimensions(A);
   3   if m < n then
   4     error "The number of rows must be greater than or equal to the number of columns to use this algorithm."
       else
   5     for j to n do
   6       V[j] := convert(LinearAlgebra:-SubMatrix(A,1 .. m,j .. j),Vector[column])
         end do;
   7     unassign('Q','R');
   8     R := Matrix(n,n,shape = triangular[upper]);
   9     R[1,1] := LinearAlgebra:-VectorNorm(V[1],2);
  10     Q[1] := V[1]/R[1,1];
  11     for j from 2 to n do
  12       for i to j do
  13         if i <> j then
  14           R[i,j] := LinearAlgebra:-DotProduct(Q[i],V[j]);
  15           V[j] := -Q[i]*R[i,j]+V[j]
             else
  16           R[i,j] := LinearAlgebra:-VectorNorm(V[j],2);
  17           Q[j] := V[j]/R[i,j]
             end if
           end do
         end do;
  18     Q := [seq(entries(Q)[k][1],k = 1 .. n)];
  19     Q := convert(Q,Matrix);
  20     total := time()-total;
  21     printf("The QR factorization of the matrix based on the MGS algorithm required %0.8f CPU seconds.",total)
       end if
end proc

Let's apply the program to a few examples to see how well it works!

> 

f:=rand(1..10);  # create a short program that randomly generates any integer from 1 to 10, inclusive.

proc () (proc () option builtin = RandNumberInterface; end proc)(6, 10, 4)+1 end proc

(2)
> 

M:=Matrix(f,8,4);

M := Matrix(8, 4, {(1, 1) = 7, (1, 2) = 10, (1, 3) = 6, (1, 4) = 2, (2, 1) = 4, (2, 2) = 6, (2, 3) = 5, (2, 4) = 1, (3, 1) = 8, (3, 2) = 5, (3, 3) = 10, (3, 4) = 2, (4, 1) = 2, (4, 2) = 4, (4, 3) = 8, (4, 4) = 3, (5, 1) = 9, (5, 2) = 10, (5, 3) = 2, (5, 4) = 8, (6, 1) = 10, (6, 2) = 9, (6, 3) = 1, (6, 4) = 6, (7, 1) = 7, (7, 2) = 7, (7, 3) = 3, (7, 4) = 3, (8, 1) = 5, (8, 2) = 3, (8, 3) = 8, (8, 4) = 10})

(3)
> 

LinearAlgebra:-Map(evalf,M);

Matrix(8, 4, {(1, 1) = 7., (1, 2) = 10., (1, 3) = 6., (1, 4) = 2., (2, 1) = 4., (2, 2) = 6., (2, 3) = 5., (2, 4) = 1., (3, 1) = 8., (3, 2) = 5., (3, 3) = 10., (3, 4) = 2., (4, 1) = 2., (4, 2) = 4., (4, 3) = 8., (4, 4) = 3., (5, 1) = 9., (5, 2) = 10., (5, 3) = 2., (5, 4) = 8., (6, 1) = 10., (6, 2) = 9., (6, 3) = 1., (6, 4) = 6., (7, 1) = 7., (7, 2) = 7., (7, 3) = 3., (7, 4) = 3., (8, 1) = 5., (8, 2) = 3., (8, 3) = 8., (8, 4) = 10.})

(4)
> 

mgs(M);

The QR factorization of the matrix based on the MGS algorithm required 0.05300000 CPU seconds.

> 

Q;

Matrix(8, 4, {(1, 1) = .355371157890000, (1, 2) = .536795101025082, (1, 3) = .125725382685606, (1, 4) = -.206838826912019, (2, 1) = .203069233080000, (2, 2) = .357255822814332, (2, 3) = .201801251387011, (2, 4) = -.148218607853882, (3, 1) = .406138466160000, (3, 2) = -.523124598971335, (3, 3) = .408554107555746, (3, 4) = -.525502414461480, (4, 1) = .101534616540000, (4, 2) = .355433089207166, (4, 3) = .555754240176830, (4, 4) = .178485820208510, (5, 1) = .456905774430000, (5, 2) = .185007479032248, (5, 3) = -.309642132170864, (5, 4) = .369047498397399, (6, 1) = .507673082700000, (6, 2) = -.167691509764169, (6, 3) = -.443789247197424, (6, 4) = 0.139117060144198e-1, (7, 1) = .355371157890000, (7, 2) = 0.637956762508179e-2, (7, 3) = -.120640311849369, (7, 4) = -.143129017130711, (8, 1) = .253836541350000, (8, 2) = -.349053521582084, (8, 3) = .400385961743801, (8, 4) = .685830489856305})

(5)
> 

R;

Matrix(4, 4, {(1, 1) = 19.6977156035922079, (1, 2) = 19.5961809922199990, (1, 3) = 12.5395251426899996, (1, 4) = 12.3364559096099988, (2, 1) = 0, (2, 2) = 5.65594295601091446, (2, 3) = 0.483024454814409143e-1, (2, 4) = -1.54658962230206076, (3, 1) = 0, (3, 2) = 0, (3, 3) = 12.0730268015982993, (3, 4) = 1.43968910037254383, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 8.38732445450181530})

(6)
> 

with(LinearAlgebra):  # load Maple's library of Linear Algebra commands and routines.

> 

map(round,Transpose(Q).Q);

Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1})

(7)
> 

map(round,Q.R);

Matrix(8, 4, {(1, 1) = 7, (1, 2) = 10, (1, 3) = 6, (1, 4) = 2, (2, 1) = 4, (2, 2) = 6, (2, 3) = 5, (2, 4) = 1, (3, 1) = 8, (3, 2) = 5, (3, 3) = 10, (3, 4) = 2, (4, 1) = 2, (4, 2) = 4, (4, 3) = 8, (4, 4) = 3, (5, 1) = 9, (5, 2) = 10, (5, 3) = 2, (5, 4) = 8, (6, 1) = 10, (6, 2) = 9, (6, 3) = 1, (6, 4) = 6, (7, 1) = 7, (7, 2) = 7, (7, 3) = 3, (7, 4) = 3, (8, 1) = 5, (8, 2) = 3, (8, 3) = 8, (8, 4) = 10})

(8)

Let's apply my program to solving a practical least squares problem.  Consider the following problem taken from the public website, http://calculator.maconstate.edu/cubic_regression/.  x = { -3, -2, -1, 0, 1, 2, 3} and y = { 3, -8, -7, 0, 7, 8, -3}.  Use the MGS (modified Gram-Schmidt) algorithm to perform a QR factorization on the matrix for this problem, and then determine the unique least squares cubic polynomial that best fit these data.

> 

X:=[seq(n,n=-3..3)];

[-3, -2, -1, 0, 1, 2, 3]

(9)
> 

Y:=[3,-8,-7,0,7,8,-3];

[3, -8, -7, 0, 7, 8, -3]

(10)
> 

A:=VandermondeMatrix(X,7,4);  # 7 rows because there are 7 points, and 4 columns because a cubic polynomial normally has 4 terms.

A := Matrix(7, 4, {(1, 1) = 1, (1, 2) = -3, (1, 3) = 9, (1, 4) = -27, (2, 1) = 1, (2, 2) = -2, (2, 3) = 4, (2, 4) = -8, (3, 1) = 1, (3, 2) = -1, (3, 3) = 1, (3, 4) = -1, (4, 1) = 1, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (5, 1) = 1, (5, 2) = 1, (5, 3) = 1, (5, 4) = 1, (6, 1) = 1, (6, 2) = 2, (6, 3) = 4, (6, 4) = 8, (7, 1) = 1, (7, 2) = 3, (7, 3) = 9, (7, 4) = 27})

(11)
> 

Map(evalf,A);  # convert A's elements to floating point numbers.

Matrix(7, 4, {(1, 1) = 1., (1, 2) = -3., (1, 3) = 9., (1, 4) = -27., (2, 1) = 1., (2, 2) = -2., (2, 3) = 4., (2, 4) = -8., (3, 1) = 1., (3, 2) = -1., (3, 3) = 1., (3, 4) = -1., (4, 1) = 1., (4, 2) = 0., (4, 3) = 0., (4, 4) = 0., (5, 1) = 1., (5, 2) = 1., (5, 3) = 1., (5, 4) = 1., (6, 1) = 1., (6, 2) = 2., (6, 3) = 4., (6, 4) = 8., (7, 1) = 1., (7, 2) = 3., (7, 3) = 9., (7, 4) = 27.})

(12)
> 

A;   # checking to make sure that my changes were saved.

Matrix(7, 4, {(1, 1) = 1., (1, 2) = -3., (1, 3) = 9., (1, 4) = -27., (2, 1) = 1., (2, 2) = -2., (2, 3) = 4., (2, 4) = -8., (3, 1) = 1., (3, 2) = -1., (3, 3) = 1., (3, 4) = -1., (4, 1) = 1., (4, 2) = 0., (4, 3) = 0., (4, 4) = 0., (5, 1) = 1., (5, 2) = 1., (5, 3) = 1., (5, 4) = 1., (6, 1) = 1., (6, 2) = 2., (6, 3) = 4., (6, 4) = 8., (7, 1) = 1., (7, 2) = 3., (7, 3) = 9., (7, 4) = 27.})

(13)
> 

mgs(A);  # mgs is my own program for doing a QR factorization that implements the modified Gram-Schmidt algorithm.

The QR factorization of the matrix based on the MGS algorithm required 0.00500000 CPU seconds.

> 

Q;

Matrix(7, 4, {(1, 1) = .377964473000000, (1, 2) = -.566946709500000, (1, 3) = .545544725521309, (1, 4) = -.408248290449766, (2, 1) = .377964473000000, (2, 2) = -.377964473000000, (2, 3) = 0.2130922499e-10, (2, 4) = .408248290333490, (3, 1) = .377964473000000, (3, 2) = -.188982236500000, (3, 3) = -.327326835278691, (3, 4) = .408248290356745, (4, 1) = .377964473000000, (4, 2) = 0., (4, 3) = -.436435780378691, (4, 4) = 0., (5, 1) = .377964473000000, (5, 2) = .188982236500000, (5, 3) = -.327326835278691, (5, 4) = -.408248290356745, (6, 1) = .377964473000000, (6, 2) = .377964473000000, (6, 3) = 0.2130926161e-10, (6, 4) = -.408248290333490, (7, 1) = .377964473000000, (7, 2) = .566946709500000, (7, 3) = .545544725521309, (7, 4) = .408248290449766})

(14)
> 

R;

Matrix(4, 4, {(1, 1) = 2.64575131106459072, (1, 2) = 0., (1, 3) = 10.5830052440000024, (1, 4) = 0., (2, 1) = 0, (2, 2) = 5.29150262212918143, (2, 3) = -0.4440892099e-15, (2, 4) = 37.0405183539999996, (3, 1) = 0, (3, 2) = 0, (3, 3) = 9.16515138991167966, (3, 4) = 0., (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 14.6969384566990673})

(15)
> 

whattype(Y);

list

(16)
> 

Y:=convert(Y,Vector[column]);

Y := Vector(7, {(1) = 3, (2) = -8, (3) = -7, (4) = 0, (5) = 7, (6) = 8, (7) = -3})

(17)
> 

Y:=Transpose(Q).Y;

Y := Vector(4, {(1) = 0., (2) = 5.291502622, (3) = 0.2220446049e-15, (4) = -14.696938453028853})

(18)
> 

M:=<R|Y>;  # augment matrix R with the product Q*Y from the line above.

M := Matrix(4, 5, {(1, 1) = 2.64575131106459072, (1, 2) = 0., (1, 3) = 10.5830052440000024, (1, 4) = 0., (1, 5) = 0., (2, 1) = 0, (2, 2) = 5.29150262212918143, (2, 3) = -0.4440892099e-15, (2, 4) = 37.0405183539999996, (2, 5) = 5.29150262200000032, (3, 1) = 0, (3, 2) = 0, (3, 3) = 9.16515138991167966, (3, 4) = 0., (3, 5) = 0.2220446049e-15, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 14.6969384566990673, (4, 5) = -14.6969384530288529})

(19)
> 

M:=ReducedRowEchelonForm(M);

M := Matrix(4, 5, {(1, 1) = 1.0, (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = -0.9690821045e-16, (2, 1) = 0., (2, 2) = 1.0, (2, 3) = 0., (2, 4) = 0., (2, 5) = 7.9999999980566106, (3, 1) = 0., (3, 2) = 0., (3, 3) = 1.0, (3, 4) = 0., (3, 5) = 0.2422705261e-16, (4, 1) = 0., (4, 2) = 0., (4, 3) = 0., (4, 4) = 1.0, (4, 5) = -.9999999997502734})

(20)
> 

p0:=M[1,5];

HFloat(-9.690821044768376e-17)

(21)
> 

p1:=M[2,5];

HFloat(7.9999999980566106)

(22)
> 

p2:=M[3,5];

HFloat(2.422705261251239e-17)

(23)
> 

p3:=M[4,5];

HFloat(-0.9999999997502734)

(24)
> 

X:=(9);  # redefine my X and Y data points.

[-3, -2, -1, 0, 1, 2, 3]

(25)
> 

Y:=(10);

[3, -8, -7, 0, 7, 8, -3]

(26)
> 

with(plots):

> 

P:=x-> p0 + p1*x + p2*x^2 + p3*x^3;   # define the least squares cubic polynomial for this problem.

proc (x) options operator, arrow; p0+p1*x+p2*x^2+p3*x^3 end proc

(27)
> 

Data:=[seq([X[n],Y[n]],n=1..7)];

[[-3, 3], [-2, -8], [-1, -7], [0, 0], [1, 7], [2, 8], [3, -3]]

(28)
> 

display(listplot(Data,style=point,symbol=solidcircle,symbolsize=15,color=black,axesfont=[Arial,bold,9]),plot(P(x),x=-3..3,color=red,thickness=2,title="\nLeast Squares Cubic Polynomial\nfound through QR Decomposition\n",titlefont=[Arial,bold,12],caption="\nGraph by Douglas Lewit of Math-577",captionfont=[Arial,bold,10],labels=["x","y"],labelfont=[Times,11]));

What happens if we feed a matrix to the 'mgs' program where the matrix has fewer rows than columns?

> 

A:=Matrix(f,4,5);

A := Matrix(4, 5, {(1, 1) = 3, (1, 2) = 1, (1, 3) = 1, (1, 4) = 3, (1, 5) = 5, (2, 1) = 5, (2, 2) = 10, (2, 3) = 10, (2, 4) = 10, (2, 5) = 2, (3, 1) = 8, (3, 2) = 4, (3, 3) = 2, (3, 4) = 9, (3, 5) = 5, (4, 1) = 9, (4, 2) = 7, (4, 3) = 4, (4, 4) = 3, (4, 5) = 10})

(29)
> 

mgs(A);

Error, (in mgs) The number of rows must be greater than or equal to the number of columns to use this algorithm.

The modified Gram-Schmidt algorithm contains the assumption that the matrix has at least as many rows as columns.  For example, in the matrix above we have a sample of five vectors from R^4, but that doesn't make any sense.  Any basis of R^4 must contain no more than four linearly independent vectors.  Any subspace of R^4 has a basis composed of four or fewer linearly independent vectors.  If we are looking at five vectors in R^4, then we are assured that at least one of those vectors is dependent or in other words, is a linear combination of the other vectors in the basis.  The Gram Schmidt orthogonalization algorithm, whether in its classical or modified form, can only be successfully applied to a linearly independent set of vectors.  

>