Contact Us!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download
Views: 20
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2204
Tweet
Share
Share
Kernel: SageMath 9.1

Projections

from Strang's Introduction to Linear Algebra, 5th ed.

Section 4.2

Projection Onto a Line

Projecting a vector b\textbf{b} onto a vector a\textbf{a} will result in a vector that is in the direction of a\textbf{a}, but scaled. The scalar can be represented by x^\hat{x}, and hence the projection can be written p=x^a\mathbf{p}=\hat{x}\mathbf{a}. The error vector e\textbf{e} is the shortest distance to a\textbf{a}, and is perpendicular.

Projecting b\textbf{b} onto a\textbf{a} with error e=b−x^a\mathbf{e} = \mathbf{b} - \hat{x}\mathbf{a}

a⋅(b−x^a) or a⋅b−x^a⋅a=0a \cdot (b-\hat{x}a) \text{ or } a \cdot b - \hat{x}a \cdot a = 0x^=a⋅ba⋅a=aTbaTa\hat{x} = \frac{a \cdot b}{a \cdot a} = \frac{a^T b}{a^T a}

Note, Sage does not distinguish between column and row vectors, and interprets them as needed in operations.

b = vector([1,1,1])
a = vector([1,2,2])
xhat = (a.dot_product(b))/(a.dot_product(a))
p = xhat*a
print("xhat = ", xhat)
print("p = ", p)
xhat = 5/9 p = (5/9, 10/9, 10/9)

x^\hat{x} is determined by both b\mathbf{b} and a\mathbf{a}. What we need next is the projection matrix PP that gives p=Pb\mathbf{p} = P \mathbf{b}. The matrix PP is a transform for projecting ANY vector onto a\mathbf{a}.

p=x^a=ax^=aaTbaTa=aaTaTab=Pbp = \hat{x}a = a \hat{x} = a \frac{a^T b}{a^T a} = \frac{a a^T}{a^T a} b = P b

Note that in the numerator we have a column times a row, therefore an n×nn\times n matrix.

P = (matrix(QQ, 3, a)*matrix(QQ, 1, a))/(a*a)
P
[1/9 2/9 2/9] [2/9 4/9 4/9] [2/9 4/9 4/9]
P * b
(5/9, 10/9, 10/9)

Projection Onto a Subspace

To find the projection p\mathbf{p} of vector b\mathbf{b} onto the subspace AA, we'll proceed as before. First find x^\hat{x} in p=x^a\mathbf{p} = \hat{x}\mathbf{a}, then find p\mathbf{p}, then find PP in p=Pb\mathbf{p} = P \mathbf{b}.

As before, the error vector e=b−Ax^\mathbf{e} = \mathbf{b} - A \hat{x} is perpendicular to the subspace.

e⋅A=A⋅e=ATe=0\mathbf{e} \cdot A = A \cdot \mathbf{e} = A^T \mathbf{e} = 0AT(b−Ax^)=0A^T (\mathbf{b}-A \hat{x}) = 0

from Strang, v.5

The combination p=x1^a1+…xn^an\mathbf{p} = \hat{x_1}\mathbf{a_1} + \dots \hat{x_n}\mathbf{a_n} that is closest to b\mathbf{b}:

Find x^(n×1) where AT(b−Ax^)=0 or ATb=ATAx^\textbf{Find } \hat{x} (n\times 1) \text{ where } A^T(\mathbf{b}-A\hat{x}) = \mathbf{0} \text{ or } A^T b = A^T A \hat{x}

solve for x^\hat{x} :

(ATA)−1ATb=(ATA)−1ATAx^(A^T A)^{-1} A^T \mathbf{b} = (A^T A)^{-1} A^T A \hat{x}x^=(ATA)−1ATb\hat{x} = (A^T A)^{-1} A^T \mathbf{b}

Find p(m×1) where p=Ax^=A(ATA)−1ATb\textbf{Find } \mathbf{p} (m\times 1) \text{ where } \mathbf{p} = A\hat{x} = A(A^T A)^{-1}A^T\mathbf{b}

Find P(m×m) where P=A(ATA)−1AT and p=Pb\textbf{Find } P (m \times m) \text{ where } P = A(A^T A)^{-1}A^T \text{ and } \mathbf{p} = P \mathbf{b}

b = vector([6, 0, 0])
A = matrix(QQ, 3, [1, 0,
                  1, 1,
                  1, 2,])

A.transpose() * A
[3 3] [3 5]
A.transpose() * b
(6, 0)
# solve for xhat
xhat = (A.transpose() * A) \ (A.transpose() * b)
xhat
(5, -3)
# find projection p
A * xhat
(5, 2, -1)
from problem set 4.2

problem 5. Find the projection matrix P onto the lines a1=⟨−1,2,2⟩a_1 = \langle -1, 2, 2\rangle and a2=⟨2,2,−1⟩a_2 = \langle 2, 2, -1\rangle. Find P1P2P_1 P_2.

a1 = vector([-1, 2, 2])
P1 = (matrix(QQ, 3, a1)*matrix(QQ, 1, a1))/(a1*a1)
P1
[ 1/9 -2/9 -2/9] [-2/9 4/9 4/9] [-2/9 4/9 4/9]
a2 = vector([2, 2, -1])
P2 = (matrix(QQ, 3, a2)*matrix(QQ, 1, a2))/(a2*a2)
P2
[ 4/9 4/9 -2/9] [ 4/9 4/9 -2/9] [-2/9 -2/9 1/9]
P1 * P2
[0 0 0] [0 0 0] [0 0 0]

problem 6. Project b=⟨1,0,0⟩b = \langle 1, 0, 0\rangle onto a1a_1, a2a_2, and a3=⟨2,−1,2⟩a_3 = \langle 2, -1, 2\rangle. Find p1+p2+p3p_1 + p_2 + p_3.

b = vector([1, 0, 0])
p1 = P1 * b
p1
(1/9, -2/9, -2/9)
p2 = P2 * b
p2
(4/9, 4/9, -2/9)
a3 = vector([2, -1, 2])
P3 = (matrix(QQ, 3, a3)*matrix(QQ, 1, a3))/(a3*a3)
p3 = P3 * b
p3
(4/9, -2/9, 4/9)
p1 + p2 + p3
(1, 0, 0)

problem 7. Show that P1+P2+P3=IP_1 + P_2 + P_3 = I.

P1 + P2 + P3
[1 0 0] [0 1 0] [0 0 1]