QR-Decomposition and Least Squares Approximation

(Part of Section 6.2 pp 343-344 and Section 6.5)

Introduction

If A is an mxn matrix, as a result of constructing an orthonormal basis for the column

space of matrix A, one can get a factorization of matrix A as a product QR where the

columns of matrix Q consist of an orthonormal basis and matrix R an upper triangular

matrix with positive entries on the main diagonal. This factorization is useful, for example,

in finding inverses, in solving systems of linear equations and in finding eigenvalues.

> with(linalg):with(linpdt);

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the names GramSchmidt, QRdecomp and leastsqrs have been redefined

[GramSchmidt, QRdecomp, leastsqrs, lsqrdemo]

>

Description of the QR-Algorithm

Example 1: Consider the matrix

> A:=matrix([[1,1,1],[1,0,2],[1,1,0]]);

A := matrix([[1, 1, 1], [1, 0, 2], [1, 1, 0]])

>

whose columns are the independent vectors

> v1:=vector([1,1,1]);
v2:=vector([1,0,1]);
v3:=vector([1,2,0]);

v1 := vector([1, 1, 1])

v2 := vector([1, 0, 1])

v3 := vector([1, 2, 0])

1. Convert the set of vectors { v[1], v[2], v[3] } into an orthonormal set.

Normalize the first vector v[1] to obtain the first orthonormal vector u[1]

> r11:=sqrt(innerprod(v1,v1));
u1:= evalm((1/r11)*v1);

r11 := sqrt(3)

u1 := vector([1/3*sqrt(3), 1/3*sqrt(3), 1/3*sqrt(3)...

>

Find a vector w[2] , that is orthogonal to u[1]

> r12:=innerprod(v2,u1);
p1:=evalm(r12*u1):
w2:=evalm(v2-p1);

r12 := 2/3*sqrt(3)

w2 := vector([1/3, -2/3, 1/3])

Normalize the vector w[2] to obtain the second orthonormal vector u[2]

> r22:=sqrt(innerprod(w2,w2)) ;
u2:=evalm((1/r22)*w2);

r22 := 1/3*sqrt(6)

u2 := vector([1/6*sqrt(6), -1/3*sqrt(6), 1/6*sqrt(6...

Compute the third orthonormal vector by computing the inner products:

> r13:=innerprod(v3,u1);
r23:=innerprod(v3,u2);
p2:=evalm(r23*u2 + r13*u1):
w3:=evalm(v3-p2);

r13 := sqrt(3)

r23 := -1/2*sqrt(6)

w3 := vector([1/2, 0, -1/2])

>

Normalize the vector w[3] to obtain the third orthonormal vector u[3]

> r33:=sqrt(innerprod(w3,w3));
u3:=evalm((1/r33)*w3);

r33 := 1/2*sqrt(2)

u3 := vector([1/2*sqrt(2), 0, -1/2*sqrt(2)])

2. Construct the matrix Q whose columns are the orthonormal vectors u[1], u[2] , and u[3]

> Q:=augment(u1,u2,u3);

Q := matrix([[1/3*sqrt(3), 1/6*sqrt(6), 1/2*sqrt(2)...

>

3.Construct an upper triangular matrix R whose entries are the r[ij] obtained in step 1

> R:=matrix([[r11,r12,r13],[0,r22,r23],[0,0,r33]]);

R := matrix([[sqrt(3), 2/3*sqrt(3), sqrt(3)], [0, 1...

>

Perform the multiplication of matrix Q with matrix R.

> multiply(Q,R);

matrix([[1, 1, 1], [1, 0, 2], [1, 1, 0]])

>

which is the matrix A.

Least Squares Approximation

Least squares problems are introduced in Section 6.5. System Ax = b does not have a

solution if b is not in the column space of A. In this case an approximate solution is sought.

This problem can be stated as:

(LSP) Find an x that minimizes the quantity || b - Ax ||

Algorithm:

If the vector b is in the column space of A, then x is the solution to Ax = b and we are done.

If b is not in the column space of A, then find the projection b[1] of b upon the column space

of A. Now b[1] is the closest vector in the column space of A to b, i.e. || b[1] - b|| < || b[2] -b||

for any other vector b[2] . Let Ax[1] = b[1] and Ax[2] = b[2] and then we have

|| Ax[1] - b|| < || Ax[2] - b||. Hence, x[1] is the solution to the least squares problem. Note that,

since b-b[1] is orthogonal to the column space of A; that is, <b-Ax[1],A> = 0

where Ax[1] = b[1] . From this we may deduce that the least square solution must satisfy:

A^T*A*x = A^T*b

This is called the normal equation . The solution x of the normal equation is called the

least squares solution to the original system Ax = b .

Example 2: Given the matrices

> A:=matrix([[4,0],[0,2],[2,2]]);
b:=matrix([[2],[0],[5]]);
x:=matrix([[x1],[x2]]);

A := matrix([[4, 0], [0, 2], [2, 2]])

b := matrix([[2], [0], [5]])

x := matrix([[x1], [x2]])

find the least squares solution to:

> evalm(A)*evalm(x)=evalm(b);

matrix([[4, 0], [0, 2], [2, 2]])*matrix([[x1], [x2]...

The system Ax = b is inconsistent. We can verify this by writing the augmented matrix

> AUG:=matrix([[4,0,2],[0,2,0],[2,2,5]]);

AUG := matrix([[4, 0, 2], [0, 2, 0], [2, 2, 5]])

and apply Gauss elimination:

> ReducedMatrix:=gausselim(AUG);

ReducedMatrix := matrix([[4, 0, 2], [0, 2, 0], [0, ...

>

The reduced augmented matrix implies that the system is inconsistent since b is not in

the column space of A. Therefore we seek an approximate solution.

Least Squares Procedure :

Premultiply the original system Ax = b by the transpose(A) to get the new system

A^T*A*x = A^T*b .

Construct the matrix A^T*A by premultiplying A by its transpose:

> trAA:=multiply (transpose(A),A);

trAA := matrix([[20, 4], [4, 8]])

>

Construct the matrix A^T*b by premultiplying b by the transpose of A:

> trAb:=multiply(transpose(A),b);

trAb := matrix([[18], [10]])

>

The new system trAA*x = trAb is consistent.

> evalm(trAA)*evalm(x) = evalm(trAb);

matrix([[20, 4], [4, 8]])*matrix([[x1], [x2]]) = ma...

>

Construct the augmented matrix of the new system

> AUG1:=matrix([[20,4,18],[4,8,10]]);

AUG1 := matrix([[20, 4, 18], [4, 8, 10]])

apply Gauss elimination

> AUG2:=gausselim(AUG1);

AUG2 := matrix([[4, 8, 10], [0, -36, -32]])

apply backsubstitution to the reduced matrix AUG[2] to get the least square solution x[1] is:

> x1 := backsub(AUG2);

x1 := vector([13/18, 8/9])

>

******************************************************************

Notice that the soltuion to || Ax[1]-b || can be written as x[1] = (A^T*A)^`-1`*A^T*b . Replace

A =QR from the QR-factorization then x[1] = R^`-1`*Q^T*b . Thus, the least squares solution

can be written as solving R*x = Q^T*b .

******************************************************************

Compute the QR-factorization of A.

> Q := matrix([[2/5*sqrt(5), -2/15*sqrt(5)], [0, 1/3*sqrt(5)], [1/5*sqrt(5), 4/15*sqrt(5)]]);
R := matrix([[2*sqrt(5), 2/5*sqrt(5)], [0, 6/5*sqrt(5)]]);

Q := matrix([[2/5*sqrt(5), -2/15*sqrt(5)], [0, 1/3*...

R := matrix([[2*sqrt(5), 2/5*sqrt(5)], [0, 6/5*sqrt...

>

>

Solve the linear system R*x = Q^T*b .

> Qtb := multiply(transpose(Q),b);

Qtb := matrix([[9/5*sqrt(5)], [16/15*sqrt(5)]])

>

> Aug := augment(R,Qtb);

Aug := matrix([[2*sqrt(5), 2/5*sqrt(5), 9/5*sqrt(5)...

>

> rref(Aug);

matrix([[1, 0, 13/18], [0, 1, 8/9]])

>

Exercise

(Section 6.5) 16 - 20