Lay Chapter 6,
Orthogonality and Least Squares

Orthogonal Matrices

In[1]:=

Clear[u] ; 

u = ({{-6, -3, 6, 1}, {-1, 2, 1, -6}, {3, 6, 3, -2}, {6, -3, 6, -1}, {2, -1, 2, 3}, {-3, 6, 3, 2}, {-2, -1, 2, -3}, {1, 2, 1, 6}}) ;

Note that the columns of u are orthogonal, and they are all of length 10.

In[3]:=

Transpose[u] . u ;

%//MatrixForm

Out[4]//MatrixForm=

( {{100, 0, 0, 0}, {0, 100, 0, 0}, {0, 0, 100, 0}, {0, 0, 0, 100}} )

Representation of Vectors in V == S S^Perp

Let's find the projection of a vector y onto S = span(u1,u2,u3).

In[5]:=

Clear[y, yHat, u1, u2, u3] ; 

y    = {4, 3, 3, -1} ; 

u1 = {1, 1, 0, 1} ;

u2 = {-1, 3, 1, -2} ;

u3 = {-1, 0, 1, 1} ;

Are the u's mutually orthogonal?

In[10]:=

u = Transpose[{u1, u2, u3}] ;

%//MatrixForm

Out[11]//MatrixForm=

( {{1, -1, -1}, {1, 3, 0}, {0, 1, 1}, {1, -2, 1}} )

In[12]:=

Transpose[u] . u ;

%//MatrixForm

Out[13]//MatrixForm=

( {{3, 0, 0}, {0, 15, 0}, {0, 0, 3}} )

Yes, the u's are mutually orthogonal

Find the projection of y onto span(u1,u2,u3).

In[14]:=

yHat = y . u1/u1 . u1u1 + y . u2/u2 . u2u2 + y . u3/u3 . u3u3 ;

%//MatrixForm

Out[15]//MatrixForm=

( {{2}, {4}, {0}, {0}} )

Then y - yHat is in U^perp

In[16]:=

y - yHat ;

%//MatrixForm

Out[17]//MatrixForm=

( {{2}, {-1}, {3}, {-1}} )

... and the required decomposition is

In[18]:=

Print["y = ", yHat, " + ", y - yHat]

y =  {2, 4, 0, 0}  +  {2, -1, 3, -1}

Let's check that.
Is yHat in span(u1,u2,u3)?

In[19]:=

RowReduce[{u1, u2, u3, yHat}] ;

%//MatrixForm

Out[20]//MatrixForm=

( {{1, 0, 0, 2}, {0, 1, 0, -1}, {0, 0, 1, 3}, {0, 0, 0, 0}} )

Is y-yHat in (span(u1, u2, u3))^perp?

In[21]:=

Transpose[u] . (y - yHat)

Out[21]=

{0, 0, 0}

Is y == yHat + (y - yHat)?
This one should be easy.

In[22]:=

y == yHat + (y - yHat)

Out[22]=

True

Projection of a Vector onto the Column Space of an Orthogonal Matrix

In[23]:=

Clear[u, y] ; 

u = ({{-6, -3, 6, 1}, {-1, 2, 1, -6}, {3, 6, 3, -2}, {6, -3, 6, -1}, {2, -1, 2, 3}, {-3, 6, 3, 2}, {-2, -1, 2, -3}, {1, 2, 1, 6}}) ;

y = {1, 1, 1, 1, 1, 1, 1, 1} ;

Note that the columns of u are orthogonal, and they are all of length 10.

In[26]:=

Transpose[u] . u ;

%//MatrixForm

Out[27]//MatrixForm=

( {{100, 0, 0, 0}, {0, 100, 0, 0}, {0, 0, 100, 0}, {0, 0, 0, 100}} )

Therefore, the closest point to y in the column space of u is yHat, where

In[28]:=

(* yHat = 1/100 ((y . u1) u1 + ... + (y . u4) u4) =   1/10 u . Transpose[u] . y *)

Let's formalize that result. Let W be the column space of U.
See Lay 6.3, Theorem 10, p399.

In[29]:=

projW[y_] := 1/100u . Transpose[u] . y

yHat = projW[y]

Out[30]=

{6/5, 2/5, 6/5, 6/5, 2/5, 6/5, 2/5, 2/5}

Check: Verify that yHat really is in the column space of u ...

In[31]:=

Solve[u . {a1, a2, a3, a4} yHat, {a1, a2, a3, a4}]

Out[31]=

{{a10, a22/25, a36/25, a40}}

...and that y-yHat is orthogonal to each of the columns of u.

In[32]:=

Transpose[u] . (y - yHat)

Out[32]=

{0, 0, 0, 0}

QR Decomposition

In[33]:=

Clear[a, q, r] ; 

a = ({{3, 8}, {0, 5}, {-1, -6}}) ;

Mathematica has built-in support for QR factorization, but its definition of the matrices involved differs somewhat from that of Lay. Consult Mathematica's Help Browser for the function QRDecomposition, and see that entry's "Further Examples."

In[35]:=

{q, r} = QRDecomposition[a] ; 

Map[MatrixForm, {q, r}]

Out[36]=

{( {{3/10^(1/2), 0, -1/10^(1/2)}, {-1/35^(1/2), 5/7^(1/2), -3/35^(1/2)}} ), ( {{10^(1/2), 3 10^(1/2)}, {0, 35^(1/2)}} )}

Note that the matrix Q of Lay is the matrix Conjugate[Transpose[q]] of Mathematica.

In[37]:=

Q = Conjugate[Transpose[q]] ;

%//MatrixForm

Out[38]//MatrixForm=

( {{3/10^(1/2), -1/35^(1/2)}, {0, 5/7^(1/2)}, {-1/10^(1/2), -3/35^(1/2)}} )

Check.

In[39]:=

a == Q . r

Out[39]=

True

QR Decomposition

In[40]:=

Clear[a, q, r] ; 

a = ({{-10, 13, 7, -11}, {2, 1, -5, 3}, {-6, 3, 13, -3}, {16, -16, -2, 5}, {2, 1, -5, -7}}) ;

Mathematica has built-in support for QR factorization, but its definition of the matrices involved differs somewhat from that of Lay. Consult Mathematica's Help Browser for the function QRDecomposition, and see that entry's "Further Examples."

In[42]:=

{q, r} = QRDecomposition[a] ; 

Map[MatrixForm, {q, r}]

Out[43]=

{( {{-1/2, 1/10, -3/10, 4/5, 1/10}, {1/2, 1/2, -1/2, 0, 1/2}, {1/3^(1/2), 0, 1/3^(1/ ... , -20, -10, 10}, {0, 6, -8, -6}, {0, 0, 6 3^(1/2), -3 3^(1/2)}, {0, 0, 0, 5 2^(1/2)}} )}

Note that the matrix Q of Lay is the matrix Conjugate[Transpose[q]] of Mathematica.

In[44]:=

Q = Conjugate[Transpose[q]] ;

%//MatrixForm

Out[45]//MatrixForm=

( {{-1/2, 1/2, 1/3^(1/2), 0}, {1/10, 1/2, 0, 1/2^(1/2)}, {-3/10, -1/2, 1/3^(1/2), 0}, {4/5, 0, 1/3^(1/2), 0}, {1/10, 1/2, 0, -1/2^(1/2)}} )

Check.

In[46]:=

a == Q . r

Out[46]=

True

Least-Squares Solution of an Inconsistent System

Find a least-squares solution to the inconsistent system ax=b, where

In[47]:=

Clear[a, b, x, x1, x2] ; 

a = ({{4, 0}, {0, 2}, {1, 1}}) ; 

b = {2, 0, 11} ;

Soln: Use the normal equations.

In[50]:=

Transpose[a] . a ;

%//MatrixForm

Out[51]//MatrixForm=

( {{17, 1}, {1, 5}} )

In[52]:=

Transpose[a] . b ;

%//MatrixForm

Out[53]//MatrixForm=

( {{19}, {11}} )

Now solve the system a^Tax =  a^Tb.

In[54]:=

x = {x1, x2} ;

soln = Solve[Transpose[a] . a . x == Transpose[a] . b, {x1, x2}][[1]]

Out[55]=

{x11, x22}

In[56]:=

x/.soln

Out[56]=

{1, 2}

Least-Squares Solution of an Inconsistent System

Find a least-squares solution to the inconsistent system ax=b, where

In[57]:=

Clear[a, b, x, x1, x2, x3, x4] ; 

a = ({{1, 1, 0, 0}, {1, 1, 0, 0}, {1, 0, 1, 0}, {1, 0, 1, 0}, {1, 0, 0, 1}, {1, 0, 0, 1}}) ; 

b = {-3, -1, 0, 2, 5, 1} ;

Soln: Use the normal equations.

In[60]:=

Transpose[a] . a ;

%//MatrixForm

Out[61]//MatrixForm=

( {{6, 2, 2, 2}, {2, 2, 0, 0}, {2, 0, 2, 0}, {2, 0, 0, 2}} )

In[62]:=

Transpose[a] . b ;

%//MatrixForm

Out[63]//MatrixForm=

( {{4}, {-4}, {2}, {6}} )

Now solve the system a^Tax =  a^Tb.

In[64]:=

x = {x1, x2, x3, x4} ;

soln = Solve[Transpose[a] . a . x == Transpose[a] . b, {x1, x2, x3, x4}]

Solve :: svars : Equations may not give solutions for all \"solve\" variables.  More…

Out[65]=

{{x13 - x4, x2 -5 + x4, x3 -2 + x4}}

Here is a clearer view of the solution.

In[66]:=

{x1, x2, x3, x4}/.soln

Out[66]=

{{3 - x4, -5 + x4, -2 + x4, x4}}

and, better still, we have

In[67]:=

(* x = {3, -5, -2, 0} + x4 {-1, 1, 1, 1} *)

Least-Squares Error

(1) Find a least-squares solution to the inconsistent system ax=b, where

In[68]:=

Clear[a, b, x, x1, x2, xHat, soln, d] ; 

a = ({{4, 0}, {0, 2}, {1, 1}}) ; 

b = {2, 0, 11} ;

Soln: Use the normal equations.

In[71]:=

Transpose[a] . a ;

%//MatrixForm

Out[72]//MatrixForm=

( {{17, 1}, {1, 5}} )

In[73]:=

Transpose[a] . b ;

%//MatrixForm

Out[74]//MatrixForm=

( {{19}, {11}} )

Now solve the system a^Tax =  a^Tb.

In[75]:=

x = {x1, x2} ;

soln = Solve[Transpose[a] . a . x == Transpose[a] . b, {x1, x2}]

Out[76]=

{{x11, x22}}

xHat is the least-squares solution. In this case it is unique because a^Ta  is invertible: xHat = (a^T a)^(-1)a^Tb.

In[77]:=

xHat = x/.soln[[1]]

Out[77]=

{1, 2}

(2) Find a least-squares error of this approximation.

In[78]:=

d = ((b - a . xHat) . (b - a . xHat))^(1/2)

Out[78]=

2 21^(1/2)

Least-Squares Solution via QR Decomposition

Find a least-squares solution to the system ax=b, where

In[79]:=

Clear[a, b, q, r, Q, xHat, x1, x2, x3, soln] ; 

a = ({{1, 3, 5}, {1, 1, 0}, {1, 1, 2}, {1, 3, 3}}) ; 

b = {3, 5, 7, -3} ;

Let's use the QR Decomposition of a to find xHat, as in Lay Theorem 15, page 414.
(1) Calculate the QR Decomposition of a.

In[82]:=

{q, r} = QRDecomposition[a] ;

Map[MatrixForm, {q, r}]

Out[83]=

{( {{1/2, 1/2, 1/2, 1/2}, {1/2, -1/2, -1/2, 1/2}, {1/2, -1/2, 1/2, -1/2}} ), ( {{2, 4, 5}, {0, 2, 3}, {0, 0, 2}} )}

Note that the matrix Q of Lay is the matrix Conjugate[Transpose[q]] of Mathematica.

In[84]:=

Q = Conjugate[Transpose[q]] ;

%//MatrixForm

Out[85]//MatrixForm=

( {{1/2, 1/2, 1/2}, {1/2, -1/2, -1/2}, {1/2, -1/2, 1/2}, {1/2, 1/2, -1/2}} )

Check.

In[86]:=

a == Q . r

Out[86]=

True

(2) Use the QR Decomposition of a to calculate xHat.

In[87]:=

xHat = Inverse[r] . Transpose[Q] . b

Out[87]=

{10, -6, 2}

Note that Lay suggests that, for numerical reasons, it is preferable to solve for xHat directly from the equation
     r.xHat==Transpose[Q].b.

In[88]:=

Clear[x1, x2, x3, soln] 

                   T     T Print[r.{x1,x2,x3 }  == Q  .b] <br />

Print[r//MatrixForm, ".", {x1, x2, x3}//MatrixForm, "==", Transpose[Q]//MatrixForm, ".", b//MatrixForm]

             T     T r.{x1,x2,x3 }  == Q  .b

( {{2, 4, 5}, {0, 2, 3}, {0, 0, 2}} )  .  ( {{x1}, { ...  {1/2, -1/2, 1/2, -1/2}} )  .  ( {{3}, {5}, {7}, {-3}} )

Let's try that.

In[91]:=

soln = Solve[r . {x1, x2, x3} Transpose[Q] . b, {x1, x2, x3}][[1]] ; 

{x1, x2, x3} = {x1, x2, x3}/.soln

Out[92]=

{10, -6, 2}


Created by Mathematica  (February 25, 2005) Valid XHTML 1.1!