Unit 6: Matrix Algebra

Chapter 33: Matrix Computations

Section 33.2: projections

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> with(linalg):
with(plots):
with(plottools):
read(`C:/Program Files/Maple 6.01/pvac.txt`):

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

Warning, the name changecoords has been redefined

>

Resolution of a Vector into Orthogonal Components

If neigher of the vectors B and A are 0 , the zero vector, we can express B in terms of its vector components along A and perpendicular to A . Each of these vector components are called projections , either B , the projection of B on, or along, A , and B , the projection (or component) of B orthogonal to A .

>

Projection of B along A

The vector (" B on A "), the projection of B on, or along, A , is given by (33.1), namely,

= A

(In words, read (33.1) as: the vector (" B on A ") is given by " B -dot- A over A -dot- A , times A ". There are four A 's on the right, and just one B . If what is "more" is "heavier," then the A 's sink to the bottom, and the lone B floats to the top. Thus, B is on the A 's, and the formula for " B on A " has but one B and multiple A 's. That is how the author himself continues to remember this particular formula.)

As an example, consider the two vectors

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

>

graphed, along with B in Figure 33.1.

> p1 := arrow([0,0], A, .05,.2,.1, color=cyan):
p2 := arrow([0,0], B, .04,.2,.1,color=black):
p3 := arrow([0,0],vector([1,0]), .1,.3,.2, color=red):
p4 := arrow([0,0],vector([0,1]), .05,.3,.2, color=green):
p5 := textplot({[1,1.1,`B`],[2,.1,`A`], [.5,-.2,`B`], [.57,-.3,`A`], [-.7,.72,`B`], [-.5,.6,`_|_A`]}, font=[TIMES,BOLD,12]):
p6 := plot([[1,0],[1,1],[0,1]], color=black, linestyle=2):
display([p||(1..6)], scaling=constrained, xtickmarks=[0,1,2], ytickmarks=[0,1]);

>

From the figure, the angle between A and B is , so the length of B is found to be

|| B ||

A vector of this length in the direction of A is 1 ( A / || A ||) = i . Alternatively,

= A = (2 i ) = i

Working in Maple, we can get the angle between A and B via

> angleAB := simplify(angle(A,B));

>

With simple right-triangle trig, the length of the red arrow B is found to be

|| B ||

Computed in Maple, this length of this red vector is

> red_length := norm(B,2)*cos(angleAB);

>

Run this length out along a normalized version of the vector A . This vector is found in Maple by

> unitizedA := normalize(A);

>

The red vector is then the length 1 run out along the normalized version of A . This is computed as

> shadow_vector := evalm(red_length * unitizedA);

>

Now, we'll use the "formula"

> evalm(dotprod(B,A)/dotprod(A,A) * A);

>

The operation of projecting B onto A can also be thought of as a function. In Maple, this function is constructed via

> projBA := (B,A) -> evalm(dotprod(B,A)/dotprod(A,A)*A);

>

This function is useful if many such projections were being computed. Otherwise, construct the occasional projection as needed.

The order of the letters matters. The projection of B on A is not the same as the projection of A on B . Indeed, computing both, we see

> B_on_A := projBA(B,A);
A_on_B := projBA(A,B);

>

Thus, the projection operation is not commutative. Switching the order of the vectors changes which vector casts the shadow. Clearly, the lengths of the projections will not be the same if the vectors are not the same lengths.

>

Projection of B Orthogonal to A

The projection, or component, of B perpendicular to A is

=

as shown in Figure 33.1, repeated here for convenience. The green arrow, the component of B which is perpendicular to A , is the vector " B perpendicular to A ."

> p1 := arrow([0,0], A, .05,.2,.1, color=cyan):
p2 := arrow([0,0], B, .04,.2,.1, color=black):
p3 := arrow([0,0],vector([1,0]), .1,.3,.2, color=red):
p4 := arrow([0,0],vector([0,1]), .05,.3,.2, color=green):
p5 := textplot({[1,1.1,`B`],[2,.1,`A`], [.5,-.2,`B`], [.57,-.3,`A`], [-.7,.72,`B`], [-.5,.6,`_|_A`]}, font=[TIMES,BOLD,12]):
display([p||(1..5)], scaling=constrained, xtickmarks=[0,1,2], ytickmarks=[0,1]);

>

The hard work is done in getting the component of B along A . The component of B perpendicular to A is obtained by a subtraction:

> BpA := evalm(B - projBA(B,A));

>

The sum of the components should restore the vector B :

> evalm(B_on_A + BpA);
print(B);

>

Also, observe that the components are clearly perpendicular:

> dotprod(B_on_A, BpA);

>

Projections - A Matrix Formulation

We next recast in matrix form, our work on projections along, and perpendicular to, a vector.

Projection of B along A

The matrix

P =

will project any vector onto the vector A . When applied to a vector B , we can write

P B = B = A = A

Remenber, A is a vector, so A B = A . B is a scalar which commutes with A , thereby validating the middle fraction.

Example 33.1

The matrix which projects any vector B onto the vector

A = i + 2 j + 3 k

that is,

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

>

is

P = =

which we obtain in Maple with

> P := evalm(A &* transpose(A)/dotprod(A,A));

>

Applied to the vector B = i + j + k

that is, to the vector

> B := vector([1,1,1]);

>

we get

P B = B = ( i + 2 j + 2 k )

that is, we get

> B_on_A := evalm(P &* B);

>

Of course, B is also given by

A = ( i + 2 j + 2 k )

that is, by

> evalm(dotprod(B,A)/dotprod(A,A) * A);

>

Projection of B Orthogonal to A

The projection, or component of B perpendicular to A is given by

= = B P B = ( ) B = (4 i j k )

that is, by

> BpA := evalm((1-P) &* B);

>

where the Maple shortcut of using for the matrix is too convenient not to use.

This calculation yields the component of B orthogonal to A since

B =

where the matrix is computed in Maple as

> `I-P` := evalm(1-P);

>

In general, if P projects onto the vector A, then projects onto the direction orthogonal to A . Figure 33.2 shows the vectors B , A , B , and B .

> p1 := arrow([0,0,0],A, [1,0,0], .1,.3,.1,color=cyan):
p2 := arrow([0,0,0],B, [1,0,0], .1,.3,.2, color=black):
p3 := arrow([0,0,0],B_on_A, [1,1,1], .2,.5,.2, color=black):
p4 := arrow([0,0,0],BpA, [0,0,1], .1,.4,.3, color=black):
p5 := textplot3d({[1,2,1.8,`A`],[1,1,.8,`B`],[4/9,2/9,-1/9,`B`], [5/9,.7,11/9,`B`]},font=[TIMES,BOLD,12], color=black):
p6 := textplot3d([4/9,.35,-.2,`_|_`],font=[TIMES,ROMAN,6], color=black):
p7 := textplot3d({[4/9,.6,-.28,`A`],[5/9,.9,10/9,`A`]}, font=[TIMES,BOLD,8], color=black):
display3d([p||(1..7)],axes=frame, scaling=constrained, labels=[x,y,`z `], labelfont=[TIMES,ITALIC,12], orientation=[-30,65], tickmarks=[[0,1],[0,1,2],[0,1,2]]);

>

Projection Matrices

A matrix P is called a projection matrix if it is symmetric and .

(If P B is the projection of B onto A , then P ( P B ), the projection of the projection onto A is still P B , so .)

The matrix P from Example 33.1, that is the matrix

> print(P);

>

is symmetric, as we see in Maple with

> equal(P, transpose(P));

>

and a computation shows , that is,

> print(P,evalm(P^2));

>

Hence, P is a projection matrix.

>

Projection onto a Subspace

Next, we consider the concept of projecting a vector onto a subspace generated by the basis vectors A , ..., A . For example, the projection of B onto the subspace spanned by A and A would be

B = B

the component of B lying in the plane formed by A and A .

By calling the vectors A , ..., A a basis , we are stating that they span the subspace (the subspace is the collection of all possible linear combinations of the basis vectors), and that they are linearly independent. This last point is very important.

Let A be the matrix whose columns are the vectors A , ..., A so that A = [ A , ..., A ].

A recipe for a matrix P that projects vectors onto this subspace is

P = A ( A A ) A

The matrix in (33.2), namely, in

P =

can be written as A ( A A ) A because A A is a scalar. This suggests (33.4), that is, P = A ( A A ) A , is an appropriate generalization of (33.2).

Example 33.2

Form a matrix which will project vectors onto the subspace spanned by the vectors

> A1 := vector([1,1,0]);
A2 := vector([0,1,0]);

>

These two vectors are a basis for the xy -plane. To form a matrix which projects onto this plane, begin by forming the matrix A.

> A := augment(A1,A2);

>

Next, form P by the recipe

> P := evalm(A &* inverse(transpose(A) &* A) &* transpose(A));

>

Why is this result not surprising?

Project the vector

> B := vector([a,b,c]);

>

onto the subspace spanned by A and A .

> evalm(P &* B);

>

Again, why should we not be surprised by this result?

>

Example 3.3

To Project the vector

> B := vector([1,2,3,4]);

>

onto the subspace spanned by the vectors

> A1 := vector([1,0,-1,1]);
A2 := vector([0,1,1,2]);

>

form the matrix A given by

> A := augment(A1,A2);

>

Next, construct the projection matrix P according to the recipe

> P := evalm(A &* inverse(transpose(A) &* A) &* transpose(A));

>

Note that

> evalm(P^2);

>

is the same as P , that is,

> equal(P,evalm(P^2));

>

Finally, obtain the desired projection of B onto this subspace spanned by and .

> B_onA := evalm(P &* B);

>

The component of B perpendendicular to this subspace is given by ( ) B . First, compute

> `I-P` := evalm(1-P);

>

then obtain

> BpA := evalm(`I-P` &* B);

>

Notice that this vector is indeed perpendicular to both A and A .

> dotprod(BpA, A1);
dotprod(BpA, A2);

>

Final Comments

If A is an r x c matrix whose columns are independent, then

P = A ( A A ) A

projects onto the column space of A . (The column space of A is the span of the columns, that is, the set of all possible linear combinations of the columns of A .) The matrix A A will be an invertible c x c matrix and P can be constructed from this recipe.

If A has rank , then A A will be singular and not invertible. The recipe for P fails. The redundant columns in A must be deleted before the projection matrix P can be constructed.

For example,

1. The r x c matrix A is 4 x 5 so has rank at most k = 4 and the 5 x 5 matrix A A is definitely singular.

2. The r x c matrix A is 5 x 4 and has maximal rank so the 4 x 4 matrix A A is invertible.

3. The r x c matrix A is 5 x 4 and has rank so the 4 x 4 matrix A A has rank , and is singular..

>

Minimizing Property of Projection

The projection of B onto the subspace spanned by A , ..., A is the component of B in the subspace spanned by

A , ..., A . The component orthogonal to this subspace is the shortest vector from the tip of B to the subspace formed by

A , ..., A . Figure 33.3 (constructed below, in Maple, after the appropriate calculations have been executed), illustrates this for the vectors

> B := vector([-5,7,6]);

>

and

> A1 := vector([1,2,1]);
A2 := vector([-1,0,2]);

>

The matrix A = [ A , A ] is then

> A := augment(A1,A2);

>

and P = A ( A A ) A is given by

> P := evalm(A &* inverse(transpose(A) &* A) &* transpose(A));

>

The projection of B onto the span of A and A is the vector P B = B given by

> B_on_A := evalm(P &* B);

>

whereas the component of B orthogonal to the subspace spanned by A and A is B P B = B is

> BpA := evalm(B - B_on_A);

>

The vector B = B is the shortest vctor from the plane formed by A and A to the tip of B .

To see this, let

Y = A + A = i + j + k

that is,

> Y := evalm(c[1]*A1 + c[2]*A2);

>

be a general vector in the plane of A and A . By varying and the vector Y moves in the plane. Correspondingly, the vector

R = B Y = i + j + k

that is,

> R := evalm(B - Y);

>

from the tip of Y to the tip of B , varies. We then seek values of and which minimize

= || R || =

(Minimizing || R || is equivalent to, but simpler than, minimizing || R ||.)

We form the function in Maple via

> f := dotprod(R,R,orthogonal);

>

The techniques of calculus are sufficient for finding the desired minimizing values of and . They must satisfy the equations

and

that is, the equations

> q1 := diff(f,c[1]) = 0;
q2 := diff(f,c[2]) = 0;

>

and are therefore as we find in Maple by

> q := solve({q1,q2},{c[1],c[2]});

>

The minimizing Y is then

Y = 2 A + 3 A = i + 4 j + 8 k = P B

as we see by

> print(B_on_A, subs(q,op(Y)));

>

The shortest vector R is then

R = B Y = B P B = B = i + 3 j 2 k

which we can obtain in Maple by

> print(BpA, subs(q,op(R)));

>

Figure 33.3 shows B ; A and A ; B , the projection of B on the plane formed by A and A ; and the component of B orthogonal to the plane of projection.

> p1 := arrow([0,0,0],B, [4,-3,2], .2,.5,.2, color=red):
p2 := arrow([0,0,0],A1, [4,-3,2], .2,.5,.3, color=black):
p3 := arrow([0,0,0],A2, [4,-3,2], .2,.5,.3, color=black):
p4 := arrow([0,0,0],B_on_A, [4,-3,2], .2,.5,.2, color=green):
p5 := arrow(convert(B_on_A,list), BpA, [1,1,1], .4,.6,.4,color=blue):
p6 := implicitplot3d(-4*x+3*y-2*z=-.1,x=-6..2,y=-2..5,z=-1..9, color=yellow, style=patchnogrid):
display3d([p||(1..6)], axes=none, scaling=constrained, labels=[x,y,z], labelfont=[TIMES,BOLD,14], orientation=[135,160]);

>

The vector B is drawn in red, the vectors A and A are in black, and the plane spanned by A and A is in yellow.

The projection of B onto the plane spanned by A and A , namely B , is in green, and the the component of B orthogonal to the plane of projection is in blue. This blue vector is the shortest value assumed by R = B Y .

>