## 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 = bp + e = be = 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̂ aa • e = 0a • (b - p) = 0a • b - a • x̂ a = 0`

Furthermore

 `p = x̂ aP = 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 = 2p = 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 `x̂`, then for `p = x̂ a`, then for `P b = p`.
Instead of doing this, Strang uses the equations:

 `P = A (AT A)-1 ATp* = P be* = 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))>>> Parray([[ 1. , 0. , 0. ], [ 0. , 0.8, 0.4], [ 0. , 0.4, 0.2]])>>> >>> p2 = np.dot(P,b)>>> p2array([[ 3. ], [ 4.8], [ 2.4]])>>> >>> e2 = b - p2>>> e2array([[ 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]])>>>`