Sunday, March 7, 2010

Strang's Linear Algebra: homework


In Strang, section 4.2 we're doing projections. Example A (p. 213) asks:

Project the vector b = (3,4,4) onto the line through a = (2,2,1), and then onto the plane that also contains a* = (1,0,0).
We're supposed to find p, the projection onto a, which is equal to x̂ (x with a little hat on it) times a. It is some fraction of a. The "error" is the part of b that is perpendicular to the projection:

x̂ a + e = b
p + e = b
e = b - p

The basic equations are:

x̂ = a • b / a • a (for vectors)
x̂ = AT b / AT A (for matrices)

AT A x̂ = AT b

This comes from

e = b - x̂ a
a • e = 0
a • (b - p) = 0
a • b - a • x̂ a = 0

Furthermore

p = x̂ a
P = A (AT A)-1 AT

For the first part, with a = (2,2,1), we just have

x̂ = (2,2,1)•(3,4,4) / (2,2,1)•(2,2,1) = 18/9 = 2

p = x̂ a = (4,4,2)
e = b - p = (3,4,4) - (4,4,2) = (-1,0,2)

We can check that

e • a (or e • x̂a) = 0



Part 2

Now we also consider a* = (1,0,0) to give a plane formed from the 2 vectors.
We construct the matrix:

A = [a a*] = [ 2  1 ]
[ 2 0 ]
[ 1 0 ]

The fundamental equation is:

AT A x̂ = AT b

We could solve for , then for p = x̂ a, then for P b = p.
Instead of doing this, Strang uses the equations:

P = A (AT A)-1 AT
p* = P b
e* = b - p*

I don't want to do the arithmetic, so we'll use numpy:

>>> import numpy as np
>>> a = np.array([2,2,1])
>>> b = np.array([3,4,4])
>>> a2 = np.array([1,0,0])
>>> a.shape = (3,1)
>>> b.shape = a2.shape = a.shape
>>> A = np.hstack([a,a2])
>>> print A
[[2 1]
[2 0]
[1 0]]
>>>
>>> Atrans = A.transpose()
>>> temp = np.linalg.inv(np.dot(Atrans,A))
>>> P = np.dot(A,np.dot(temp,Atrans))
>>> P
array([[ 1. , 0. , 0. ],
[ 0. , 0.8, 0.4],
[ 0. , 0.4, 0.2]])
>>>
>>> p2 = np.dot(P,b)
>>> p2
array([[ 3. ],
[ 4.8],
[ 2.4]])
>>>
>>> e2 = b - p2
>>> e2
array([[ 0. ],
[-0.8],
[ 1.6]])
>>>

We confirm that e* (e2 in the code) is perpendicular to both a and a*.
Also, for reasons that I don't understand yet, P2 = P.
[UPDATE: The reason is simple. Suppose we do P b to get the projection of b in the plane. What happens if we do it again? Answer: nothing, we're already in the plane! So P b must equal P P b.

>>> np.dot(e2.transpose(),a)
array([[ 2.66453526e-15]])
>>> np.dot(e2.transpose(),a2)
array([[ 0.]])
>>>
>>> np.dot(P,P)
array([[ 1. , 0. , 0. ],
[ 0. , 0.8, 0.4],
[ 0. , 0.4, 0.2]])
>>>

No comments: