Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

2.10. Projections, Revisited

We’re almost ready to return to our original motivation for studying linear algebra, which was to perform linear regression using multiple input variables. This section outlines the final piece of the puzzle.

Approximating using a Single Vector

In Chapter 2.3, we introduced the approximation problem, which asked:

Among all vectors of the form kx⃗k \color{#3d81f6} \vec x, which one is closest to y⃗{\color{orange}\vec y}?

We now know the answer is the vector p⃗\color{#004d40} \vec p, where

p⃗=(y⃗⋅x⃗x⃗⋅x⃗)x⃗{\color{#004d40} \vec p} = \left( \frac{{\color{orange}\vec y} \cdot \color{#3d81f6} \vec x}{{\color{#3d81f6}\vec x} \cdot {\color{#3d81f6}\vec x}} \right) \color{#3d81f6} \vec x

p⃗\color{#004d40} \vec p is called the orthogonal projection of y⃗\color{orange} \vec y onto x⃗\color{#3d81f6} \vec x.

Note that I’ve used y⃗\color{orange} \vec y and x⃗\color{#3d81f6} \vec x here rather than u⃗\color{orange} \vec u and v⃗\color{#3d81f6} \vec v, just to make the notation more consistent with the notation we’ll use as we move back into the world of machine learning.

Loading...
Loading...
Image produced in Jupyter

As we’ve studied, the resulting error vector,

e⃗=y⃗−p⃗{\color{#d81a60} \vec e} = {\color{orange} \vec y} - {\color{#004d40} \vec p}

is orthogonal to x⃗\color{#3d81f6} \vec x.

In our original look at the approximation problem, we were approximating y⃗\color{orange} \vec y using a scalar multiple of just a single vector, x⃗\color{#3d81f6} \vec x. The set of all scalar multiples of x⃗\color{#3d81f6} \vec x, denoted by span({x⃗})\text{span}(\{ {\color{#3d81f6} \vec x}\}), is a line in Rn\mathbb{R}^n.

Key idea: instead of projecting onto the subspace spanned by just a single vector, how might we project onto the subspace spanned by multiple vectors?

Approximating using Multiple Vectors

Equipped with our understanding of linear independence, spans, subspaces, and column spaces, we’re ready to tackle a more advanced version of the approximation problem.

All three statements at the bottom of the box above are asking the exact same question; I’ve presented all three forms so that you see more clearly how the ideas of spans, column spaces, and matrix-vector multiplication fit together. I will tend to refer to the latter two versions of the problem the most. In what follows, suppose X\color{#3d81f6} X is an n×dn \times d matrix whose columns x⃗(1){\color{#3d81f6} \vec x^{(1)}}, x⃗(2){\color{#3d81f6} \vec x^{(2)}}, ..., x⃗(d){\color{#3d81f6} \vec x^{(d)}} are the building blocks we want to approximate y⃗\color{orange} \vec y with.

First, let’s get the trivial case out of the way. If y⃗∈colsp(X){\color{orange} \vec y} \in \text{colsp}({\color{#3d81f6} X}), then the vector in colsp(X)\text{colsp}({\color{#3d81f6} X}) that is closest to y⃗\color{orange} \vec y is just y⃗\color{orange} \vec y itself. If that’s the case, there exists some w⃗\vec w such that y⃗=Xw⃗{\color{orange} \vec y} = {\color{#3d81f6} X} \vec w exactly. This w⃗\vec w is unique only if X{\color{#3d81f6} X}'s columns are linearly independent; otherwise, there will be infinitely many good w⃗\vec w’s.

But, that’s not the case I’m really interested in. I care more about when y⃗\color{orange} \vec y is not in colsp(X)\text{colsp}({\color{#3d81f6} X}). (Remember, this is the case we’re interested in when we’re doing linear regression: usually, it’s not possible to make our predictions 100% correct, and we’ll have to settle for some error.)

Then what?

In general, colsp(X)\text{colsp}({\color{#3d81f6} X}) is an rr-dimensional subspace of Rn\mathbb{R}^n, where r=rank(X)r = \text{rank}({\color{#3d81f6} X}). In the diagram below, I’ve used a plane to represent colsp(X)\text{colsp}({\color{#3d81f6} X}); just remember that X{\color{#3d81f6} X} may have more than 3 rows or columns.

Loading...

Remember that colsp(X)\text{colsp}({\color{#3d81f6} X}) is the set of linear combinations of X{\color{#3d81f6} X}'s columns, so it’s the set of all vectors that can be written as Xw⃗{\color{#3d81f6} X} \vec w, where w⃗∈Rd\vec w \in \mathbb{R}^d.

Let’s consider two possible vectors of the form Xw⃗{\color{#3d81f6} X} \vec w, and look at their corresponding error vectors, e⃗=y⃗−Xw⃗{\color{#d81a60} \vec e} = {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w. I won’t draw the columns of X{\color{#3d81f6} X}, since those would clutter up the picture.

Loading...

Our problem boils down to finding the w⃗\vec w that minimizes the norm of the error vector. Since it’s a bit easier to work with squared norms (remember that ∥x⃗∥2=x⃗⋅x⃗\lVert \vec x \rVert^2 = \vec x \cdot \vec x), we’ll minimize the squared norm of the error vector instead; this is an equivalent problem, since the norm is non-negative to begin with.

∥e⃗∥2=∥y⃗−Xw⃗∥2⏟which w⃗ minimizes this?\underbrace{\lVert {\color{#d81a60} \vec e} \rVert^2 = \lVert {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w \rVert^2}_{\text{which $\vec w$ minimizes this?}}

Think of ∥y⃗−Xw⃗∥2\lVert {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w \rVert^2 as a function of w⃗\vec w only; X\color{#3d81f6} X and y⃗\color{orange} \vec y should be thought of as fixed. This is a least squares problem: we’re looking for the w⃗\vec w that minimizes the sum of squared errors between y⃗\color{orange} \vec y and Xw⃗{\color{#3d81f6} X} \vec w.

There are two ways we’ll minimize this function of w⃗\vec w:

  1. Using a geometric argument, as we did in the single vector case.

  2. Using calculus. This is more involved than before, since the input variable is a vector, not a scalar, but it can be done, as we’ll see in Chapter 4.

Let’s focus on the geometric argument. What does our intuition tell us? Extending the single vector case, we expect the vector in colsp(X)\text{colsp}({\color{#3d81f6} X}) that is closest to y⃗\color{orange} \vec y to be the orthogonal projection of y⃗\color{orange} \vec y onto colsp(X)\text{colsp}({\color{#3d81f6} X}): that is, its error should be orthogonal to colsp(X)\text{colsp}({\color{#3d81f6} X}).

We could see this intuitively in the visual above. w⃗o\vec w_\text{o} was chosen to make e⃗o\color{#d81a60}\vec e_\text{o} orthogonal to colsp(X)\text{colsp}({\color{#3d81f6} X}), meaning that e⃗o\color{#d81a60}\vec e_\text{o} is orthogonal to every vector in colsp(X)\text{colsp}({\color{#3d81f6} X}). (The subscript “o” stands for “orthogonal”.) w⃗′\vec w' was some other arbitrary vector, leading e⃗′\color{#d81a60}\vec e' to not be orthogonal to colsp(X)\text{colsp}({\color{#3d81f6} X}). Clearly e⃗o\color{#d81a60}\vec e_\text{o} is shorter than e⃗′\color{#d81a60}\vec e'.

To prove that the optimal choice of w⃗\vec w comes from making the error vector orthogonal to colsp(X)\text{colsp}({\color{#3d81f6} X}), you could use the same argument as in the single vector case: if you consider two vectors, w⃗o\vec w_\text{o} with an orthogonal error vector e⃗o\color{#d81a60}\vec e_\text{o}, and w⃗′\vec w' with an error vector e⃗′\color{#d81a60}\vec e' that is not orthogonal to colsp(X)\text{colsp}({\color{#3d81f6} X}), then we can draw a right triangle with e⃗′\color{#d81a60}\vec e' as the hypotenuse and e⃗o\color{#d81a60}\vec e_\text{o} as one of the legs, making

∥e⃗′∥2>∥e⃗o∥2\lVert {\color{#d81a60}\vec e'} \rVert^2 > \lVert {\color{#d81a60}\vec e_\text{o}} \rVert^2

This is such an important idea that I want to redraw the picture above with just the orthogonal projection. Note that I’ve replaced w⃗o\vec w_\text{o} with w⃗∗\vec w^* to indicate that this is the optimal choice of w⃗\vec w.

Image produced in Jupyter

A Proof that the Orthogonal Error Vector is the Shortest

In office hours, a student asked for more justification that the shortest possible error vector is one that’s orthogonal to the column space of X\color{#3d81f6} X, especially because it’s hard to visualize what orthogonality looks like in higher dimensions. Remember, all it means for two vectors to be orthogonal is that their dot product is 0.

Given that, let’s assume only that

  • w⃗o\vec w_\text{o} is chosen so that e⃗o=y⃗−Xw⃗o{\color{#d81a60} \vec e_\text{o}} = {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w_\text{o} is orthogonal to colsp(X)\text{colsp}({\color{#3d81f6} X}), and

  • w⃗′\vec w' is any other choice of w⃗\vec w, with a corresponding error vector e⃗′=y⃗−Xw⃗′{\color{#d81a60} \vec e'} = {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w'.

Just with these facts alone, we can show that e⃗o\color{#d81a60} \vec e_\text{o} is the shortest possible error vector. To do so, let’s start by considering the (squared) magnitude of e⃗′\color{#d81a60} \vec e':

∥e⃗′∥2=∥y⃗−Xw⃗′∥2=∥y⃗−Xw⃗′+Xw⃗o−Xw⃗o∥2⏟seems arbitrary, but it’s a legal operation that brings w⃗oback in=∥y⃗−Xw⃗o⏟‘‘a"+X(w⃗o−w⃗′)⏟‘‘b"∥2=∥y⃗−Xw⃗o∥2+∥X(w⃗o−w⃗′)∥2+2(y⃗−Xw⃗o)⋅X(w⃗o−w⃗′)⏟∥a⃗+b⃗∥2=∥a⃗∥2+∥b⃗∥2+2a⃗⋅b⃗=∥e⃗o∥2+∥X(w⃗o−w⃗′)∥2+2e⃗o⋅X(w⃗o−w⃗′)⏟0, because e⃗o is orthogonal to the columns of X=∥e⃗o∥2+∥X(w⃗o−w⃗′)∥2≥∥e⃗o∥2\begin{align*} \lVert {\color{#d81a60} \vec e'} \rVert^2 &= \lVert {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w' \rVert^2 \\ &= \underbrace{\lVert {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w' + {\color{#3d81f6} X} \vec w_\text{o} - {\color{#3d81f6} X} \vec w_\text{o} \rVert^2}_{\text{seems arbitrary, but it's a legal operation that brings } \vec w_\text{o} \text{back in}} \\ &= \lVert \underbrace{{\color{orange} \vec y} - {\color{#3d81f6} X} \vec w_\text{o}}_{``a"} + \underbrace{{\color{#3d81f6} X} (\vec w_\text{o} - \vec w')}_{``b"} \rVert^2 \\ &= \underbrace{\lVert {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w_\text{o} \rVert^2 + \lVert {\color{#3d81f6} X} (\vec w_\text{o} - \vec w') \rVert^2 + 2 ({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w_\text{o}) \cdot {\color{#3d81f6} X} (\vec w_\text{o} - \vec w')}_{\lVert \vec a + \vec b \rVert^2 = \lVert \vec a \rVert^2 + \lVert \vec b \rVert^2 + 2 \vec a \cdot \vec b} \\ &= \lVert {\color{#d81a60} \vec e_\text{o}} \rVert^2 + \lVert {\color{#3d81f6} X} (\vec w_\text{o} - \vec w') \rVert^2 + \underbrace{2 {\color{#d81a60} \vec e_\text{o}} \cdot {\color{#3d81f6} X} (\vec w_\text{o} - \vec w')}_{\text{0, because } {\color{#d81a60} \vec e_\text{o}} \text{ is orthogonal to the columns of } \color{#3d81f6} X} \\ &= \lVert {\color{#d81a60} \vec e_\text{o}} \rVert^2 + \lVert {\color{#3d81f6} X} (\vec w_\text{o} - \vec w') \rVert^2 \\ \geq \lVert {\color{#d81a60} \vec e_\text{o}} \rVert^2 \end{align*}

So, no matter what choice of w⃗′\vec w' we make, the magnitude of e⃗′\color{#d81a60} \vec e' can’t be smaller than the magnitude of e⃗o\color{#d81a60} \vec e_\text{o}. This means that the error vector that is orthogonal to the column space of X\color{#3d81f6} X is the shortest possible error vector.

This is really just the same proof as in Chapter 2.3, where we argued that e⃗o{\color{#d81a60} \vec e_\text{o}}, X(w⃗o−w⃗′){\color{#3d81f6} X}(\vec w_\text{o} - \vec w'), and e⃗′\color{#d81a60} \vec e' form a right triangle, where e⃗′\color{#d81a60} \vec e' is the hypotenuse.

The Normal Equation

We’ve come to the conclusion that in order to find the w⃗\vec w that minimizes

∥e⃗∥2=∥y⃗−Xw⃗∥2\lVert {\color{#d81a60} \vec e} \rVert^2 = \lVert {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w \rVert^2

we need to find the w⃗\vec w that makes the error vector e⃗=y⃗−Xw⃗{\color{#d81a60} \vec e} = {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w orthogonal to colsp(X)\text{colsp}({\color{#3d81f6} X}). colsp(X)\text{colsp}({\color{#3d81f6} X}) is the set of all linear combinations of X{\color{#3d81f6} X}'s columns. So, if we can find an e⃗\color{#d81a60} \vec e that is orthogonal to every column of X{\color{#3d81f6} X}, then it must be orthogonal to any of their linear combinations, too.

So, we’re looking for a e⃗=y⃗−Xw⃗{\color{#d81a60} \vec e} = {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w that satisfies

x⃗(1)⋅(y⃗−Xw⃗)=0x⃗(2)⋅(y⃗−Xw⃗)=0⋮x⃗(d)⋅(y⃗−Xw⃗)=0\begin{align*} {\color{#3d81f6} \vec x^{(1)}} \cdot ({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) &= 0 \\ {\color{#3d81f6} \vec x^{(2)}} \cdot ({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) &= 0 \\ &\vdots \\ {\color{#3d81f6} \vec x^{(d)}} \cdot ({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) &= 0 \end{align*}

As you might have guessed, there’s an easier way to write these dd equations simultaneously. Above, we’re taking the dot product of y⃗−Xw⃗{\color{orange} \vec y} - {\color{#3d81f6} X} \vec w with each of the columns of X\color{#3d81f6} X. We’ve learned that Av⃗A \vec v contains the dot products of v⃗\vec v with the rows of AA. So how do we get the dot products of y⃗−Xw⃗{\color{orange} \vec y} - {\color{#3d81f6} X} \vec w with the columns of X\color{#3d81f6} X? Transpose it!

XT(y⃗−Xw⃗)=[—x⃗(1)T——x⃗(2)T—⋮—x⃗(d)T—](y⃗−Xw⃗)=[x⃗(1)⋅(y⃗−Xw⃗)x⃗(2)⋅(y⃗−Xw⃗)⋮x⃗(d)⋅(y⃗−Xw⃗)]{\color{#3d81f6} X^T}({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) = {\color{#3d81f6}\begin{bmatrix} {\color{#3d81f6} —} & {\color{#3d81f6} \vec{x}^{(1)}}^T & {\color{#3d81f6} —} \\ {\color{#3d81f6} —} & {\color{#3d81f6} \vec{x}^{(2)}}^T & {\color{#3d81f6} —} \\ & \vdots & \\ {\color{#3d81f6} —} & {\color{#3d81f6} \vec{x}^{(d)}}^T & {\color{#3d81f6} —} \end{bmatrix}} ({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) = \begin{bmatrix} {\color{#3d81f6}\vec x^{(1)}} \cdot ({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) \\ {\color{#3d81f6}\vec x^{(2)}} \cdot ({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) \\ \vdots \\ {\color{#3d81f6}\vec x^{(d)}} \cdot ({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) \end{bmatrix}

So, if we want y⃗−Xw⃗{\color{orange} \vec y} - {\color{#3d81f6} X} \vec w to be orthogonal to each of the columns of X\color{#3d81f6} X, then we need XT(y⃗−Xw⃗)=0⃗{\color{#3d81f6} X^T}({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) = \vec 0 (note that this is the vector 0⃗∈Rd\vec 0 \in \mathbb{R}^d, not the scalar 0.) Another way of saying this is that we need the error vector to be in the left null space of X\color{#3d81f6} X, i.e. e⃗∈null(XT){\color{#d81a60} \vec e} \in \text{null}({\color{#3d81f6} X^T}).

XTe⃗=0XT(y⃗−Xw⃗)=0⃗XTy⃗−XTXw⃗=0⃗XTXw⃗=XTy⃗\begin{align*} {\color{#3d81f6}X^T}{\color{#d81a60} \vec e} &= 0 \\{\color{#3d81f6} X^T}({\color{orange} \vec y} - {\color{#3d81f6} X} \vec w) &= \vec 0 \\ {\color{#3d81f6} X^T} {\color{orange} \vec y} - {\color{#3d81f6} X^T} {\color{#3d81f6} X} \vec w &= \vec 0 \\ {\color{#3d81f6} X^T} {\color{#3d81f6} X} \vec w &= {\color{#3d81f6} X^T} {\color{orange} \vec y} \end{align*}

The final equation above is called the normal equation. “Normal” means “orthogonal”. Sometimes it’s called the normal equations to reference the fact that it’s a system of dd equations and dd unknowns, where the unknowns are components of w⃗\vec w (w1,w2,…,wdw_1, w_2, \ldots, w_d.)

Note that XTXw⃗=XTy⃗{\color{#3d81f6} X^T} {\color{#3d81f6} X} \vec w = {\color{#3d81f6} X^T} {\color{orange} \vec y} looks a lot like Xw⃗=y⃗{\color{#3d81f6} X} \vec w = \color{orange} \vec y, with added factors of XT\color{#3d81f6} X^T on the left. Remember that if y⃗\color{orange} \vec y is in colsp(X)\text{colsp}({\color{#3d81f6} X}), then Xw⃗=y⃗{\color{#3d81f6} X} \vec w = \color{orange} \vec y has a solution, but that’s usualy not the case, hence why we’re attempting to approximate y⃗\color{orange} \vec y with a linear combination of X\color{#3d81f6} X’s columns.

Is there a unique vector w⃗\vec w that satisfies the normal equation? That depends on whether XTX\color{#3d81f6} X^TX is invertible. XTX\color{#3d81f6} X^TX is a d×dd \times d matrix with the same rank as the n×dn \times d matrix X\color{#3d81f6} X, as we proved in Chapter 2.8.

rank(XTX)=rank(X)\text{rank}({\color{#3d81f6} X^TX}) = \text{rank}({\color{#3d81f6} X})

So, XTX\color{#3d81f6} X^TX is invertible if and only if rank(X)=d\text{rank}({\color{#3d81f6} X}) = d, meaning all of X\color{#3d81f6} X’s columns are linearly independent. In that case, the best choice of w⃗\vec w is the unique vector

w⃗∗=(XTX)−1XTy⃗\boxed{\vec w^* = ({\color{#3d81f6} X^TX})^{-1} \color{#3d81f6}X^T \color{orange} \vec y}

w⃗∗\vec w^* has a star on it, denoting that it is the best choice of w⃗\vec w. I don’t ask you to memorize much (you get to bring a notes sheet into your exams, after all), but this equation is perhaps the most important of the semester! It might even look familiar: back in the single vector case in Chapter 2.3, the optimal coefficient on x⃗\vec x was x⃗⋅y⃗x⃗⋅x⃗=x⃗Ty⃗x⃗Tx⃗\frac{\vec x \cdot \vec y}{\vec x \cdot \vec x} = \frac{\vec x^T \vec y}{\vec x^T \vec x}, which looks similar to the one above. The difference is that here, XTX\color{#3d81f6} X^T X is a matrix, not a scalar. (But, if X\color{#3d81f6} X is just a matrix with a single column, then XTX\color{#3d81f6} X^T X is just the dot product of X\color{#3d81f6} X with itself, which is a scalar, and the boxed formula above reduces to the formula from Chapter 2.3.)

What if XTX\color{#3d81f6} X^TX isn’t invertible? Then, there are infinitely many w⃗∗\vec w^*'s that satisfy the normal equation,

XTXw⃗=XTy⃗{\color{#3d81f6} X^T} {\color{#3d81f6} X} \vec w = {\color{#3d81f6} X^T} {\color{orange} \vec y}

It’s not immediately obvious what it means for there to be infinitely many solutions to the normal equation; I’ve dedicated a whole subsection to it below to give this idea the consideration it deserves.

In the examples that follow, we’ll look at how to find all of the solutions to this equation when there are infinitely many.

First, let’s start with a straightforward example. Let

X=[1021110−1],y⃗=[1045]{\color{#3d81f6} X = \begin{bmatrix} 1 & 0 \\ 2 & 1 \\ 1 & 1 \\ 0 & -1 \end{bmatrix}}, \quad {\color{orange} \vec y = \begin{bmatrix} 1 \\ 0 \\ 4 \\ 5 \end{bmatrix}}

The vector in colsp(X)\text{colsp}({\color{#3d81f6} X}) that is closest to y⃗\color{orange} \vec y is the vector Xw⃗∗{\color{#3d81f6} X} \vec w^*, where w⃗∗\vec w^* is the solution to the normal equations,

XTXw⃗∗=XTy⃗{\color{#3d81f6} X^T X} \vec w^* = {\color{#3d81f6} X^T} \color{orange} \vec y

The first step is to compute XTX{\color{#3d81f6} X^T X}, which is a 2×22 \times 2 matrix of dot products of the columns of X{\color{#3d81f6} X}.

XTX=[6333]{\color{#3d81f6} X^T X = \begin{bmatrix} 6 & 3 \\ 3 & 3 \end{bmatrix}}

XTX{\color{#3d81f6} X^T X} is invertible, so we can solve for w⃗∗\vec w^* uniquely. Remember that in practice, we’d ask Python to solve np.linalg.solve(X.T @ X, X.T @ y), but here XTX\color{#3d81f6} X^T X is small enough that we can invert it by hand.

XTX=[6333]  ⟹  (XTX)−1=19[3−3−36]=[1/3−1/3−1/32/3]{\color{#3d81f6} X^T X = \begin{bmatrix} 6 & 3 \\ 3 & 3 \end{bmatrix}} \implies ({\color{#3d81f6} X^T X})^{-1} = \frac{1}{9} \begin{bmatrix} 3 & -3 \\ -3 & 6 \end{bmatrix} = \begin{bmatrix} 1/3 & -1/3 \\ -1/3 & 2/3 \end{bmatrix}

Then,

w⃗∗=(XTX)−1XTy⃗=[1/3−1/3−1/32/3]⏟(XTX)−1[1210011−1]⏟XT[1045]=[2−8/3]\vec w^* = ({\color{#3d81f6} X^T X})^{-1} {\color{#3d81f6} X^T} {\color{orange} \vec y} = \underbrace{\begin{bmatrix} 1/3 & -1/3 \\ -1/3 & 2/3 \end{bmatrix}}_{({\color{#3d81f6} X^T X})^{-1}} \underbrace{{\color{#3d81f6}\begin{bmatrix} 1 & 2 & 1 & 0 \\ 0 & 1 & 1 & -1 \end{bmatrix}}}_{\color{#3d81f6}X^T} {\color{orange}\begin{bmatrix} 1 \\ 0 \\ 4 \\ 5 \end{bmatrix}} = \begin{bmatrix} 2 \\ -8/3 \end{bmatrix}

The magic is in the interpretation of the numbers in w⃗∗\vec w^*, 2 and −83-\frac{8}{3}. These are the coefficients of the columns of X{\color{#3d81f6} X} in the linear combination that is closest to y⃗\color{orange} \vec y. Meaning,

Xw⃗∗=2[1210]⏟x⃗(1)−73[011−1]⏟x⃗(2)=[25/3−1/37/3]{\color{#3d81f6} X} \vec w^* = 2 \underbrace{\color{#3d81f6}\begin{bmatrix} 1 \\ 2 \\ 1 \\ 0 \end{bmatrix}}_{\color{#3d81f6}\vec x^{(1)}} - \frac{7}{3} \underbrace{\color{#3d81f6}\begin{bmatrix} 0 \\ 1 \\ 1 \\ -1 \end{bmatrix}}_{\color{#3d81f6}\vec x^{(2)}} = \begin{bmatrix} 2 \\ 5/3 \\ -1/3 \\ 7/3 \end{bmatrix}

is the vector in colsp(X)\text{colsp}({\color{#3d81f6} X}) that is closest to y⃗\color{orange} \vec y. This vector is the orthogonal projection of y⃗\color{orange} \vec y onto colsp(X)\text{colsp}({\color{#3d81f6} X}).

Examples

The first example above is the most concrete. The examples that follow will build our understanding of how orthogonal projections really work.

Example: Point and Plane

Find the point on the plane 6x−3y+2z=06x - 3y + 2z = 0 that is closest to the point (1,1,1)(1, 1, 1).

Solution

At first, this doesn’t seem like a projection problem, but it is. The plane

6x−3y+2z=06x - 3y + 2z = 0

is a 2-dimensional subspace of R3\mathbb{R}^3, meaning it can be described as the span of two non-collinear vectors. So, all we need to do is find some matrix XX with those two vectors as columns, and then use the formula we derived above to project y⃗=[111]\vec y = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} onto colsp(X)\text{colsp}(X). The fact that we said “point” instead of “vector” doesn’t change the problem: in settings like these, points and vectors are equivalent.

Two vectors that lie in the plane (but point in different directions) are [120]\begin{bmatrix} 1 \\ 2 \\ 0 \end{bmatrix} and [10−3]\begin{bmatrix} 1 \\ 0 \\ -3 \end{bmatrix}. There’s nothing special about these two vectors, other than that they have relatively small integer components. Let’s use them as columns of XX:

X=[11200−3]X = \begin{bmatrix} 1 & 1 \\ 2 & 0 \\ 0 & -3 \end{bmatrix}

Now, to formulate the normal equations, XTXw⃗∗=XTy⃗X^TX \vec w^* = X^T \vec y, we need to compute XTXX^TX and XTy⃗X^T \vec y.

XTX=[12010−3][11200−3]=[51110]X^TX = \begin{bmatrix} 1 & 2 & 0 \\ 1 & 0 & -3 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 2 & 0 \\ 0 & -3 \end{bmatrix} = \begin{bmatrix} 5 & 1 \\ 1 & 10 \end{bmatrix}
XTy⃗=[12010−3][111]=[3−2]X^T \vec y = \begin{bmatrix} 1 & 2 & 0 \\ 1 & 0 & -3 \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} = \begin{bmatrix} 3 \\ -2 \end{bmatrix}

So, the normal equations are

[51110][w0∗w1∗]⏟w⃗∗=[3−2]\begin{bmatrix} 5 & 1 \\ 1 & 10 \end{bmatrix} \underbrace{\begin{bmatrix} w_0^* \\ w_1^* \end{bmatrix}}_{\vec w^*} = \begin{bmatrix} 3 \\ -2 \end{bmatrix}

We could use w⃗∗=(XTX)−1XTy⃗\vec w^* = (X^TX)^{-1}X^T \vec y, but often it’s easier to solve the system directly.

5w0∗+w1∗=3w0∗+10w1∗=−2\begin{align*} 5w_0^* + w_1^* &= 3 \\ w_0^* + 10w_1^* &= -2 \end{align*}

Solving this system gives us w0∗=3249w_0^* = \frac{32}{49} and w1∗=−1349w_1^* = -\frac{13}{49}. (Sorry, I know the numbers aren’t pretty in this example. But that’s what happens in the real world.)

Now that we’ve solved for w⃗∗=[3249−1349]\vec w^* = \begin{bmatrix} \frac{32}{49} \\ -\frac{13}{49} \end{bmatrix}, the projection of y⃗\vec y onto colsp(X)\text{colsp}(X) – which is the point on the plane that’s closest to y⃗\vec y – is

Xw⃗∗=[11200−3][3249−1349]=[3249−134964493949]=[194964493949]\begin{align*} X\vec w^* &= \begin{bmatrix} 1 & 1 \\ 2 & 0 \\ 0 & -3 \end{bmatrix} \begin{bmatrix} \frac{32}{49} \\ -\frac{13}{49} \end{bmatrix} \\ &= \begin{bmatrix} \frac{32}{49} - \frac{13}{49} \\ \frac{64}{49} \\ \frac{39}{49} \end{bmatrix} \\ &= \begin{bmatrix} \frac{19}{49} \\ \frac{64}{49} \\ \frac{39}{49} \end{bmatrix} \end{align*}

So, after all that work, we can conclude that the point on the plane 6x−3y+2z=06x - 3y + 2z = 0 that is closest to (1,1,1)(1, 1, 1) is (1949,6449,3949)(\frac{19}{49}, \frac{64}{49}, \frac{39}{49}).

Example: What if y⃗∈colsp(X){\color{orange} \vec y} \in \text{colsp}({\color{#3d81f6} X})?

Find the orthogonal projection of y⃗=[132−1]\color{orange} \vec y = \begin{bmatrix} 1 \\ 3 \\ 2 \\ -1 \end{bmatrix} onto colsp(X)\text{colsp}({\color{#3d81f6} X}), where X=[1021110−1]{\color{#3d81f6} X = \begin{bmatrix} 1 & 0 \\ 2 & 1 \\ 1 & 1 \\ 0 & -1 \end{bmatrix}}.

Solution

First, notice that y⃗\color{orange} \vec y is just the sum of the two columns of colsp(X)\text{colsp}({\color{#3d81f6} X}). So intuitively, because y⃗{\color{orange} \vec y} is already in the column space of X{\color{#3d81f6} X}, the projection is just y⃗{\color{orange} \vec y} itself.

But let’s make sure the math works out that way.

First, we solve the normal equations for w⃗∗\vec w^*.

XTX=[1210011−1][1021110−1]=[1+4+12+12+11+1+1]=[6333]XTy⃗=[1210011−1][132−1]=[96]\begin{align*} {\color{#3d81f6} X^TX} &= {\color{#3d81f6}\begin{bmatrix} 1 & 2 & 1 & 0 \\ 0 & 1 & 1 & -1 \end{bmatrix}\begin{bmatrix} 1 & 0 \\ 2 & 1 \\ 1 & 1 \\ 0 & -1 \end{bmatrix}} \\&= {\color{#3d81f6}\begin{bmatrix} 1 + 4 + 1 & 2 + 1 \\ 2 + 1 & 1 + 1 + 1 \end{bmatrix}} \\&= {\color{#3d81f6}\begin{bmatrix} 6 & 3 \\ 3 & 3 \end{bmatrix}} \\ \\ {\color{#3d81f6}X^T}{\color{orange}\vec y} &= {\color{#3d81f6}\begin{bmatrix} 1 & 2 & 1 & 0 \\ 0 & 1 & 1 & -1 \end{bmatrix}}{\color{orange}\begin{bmatrix} 1 \\ 3 \\ 2 \\ -1 \end{bmatrix}} \\&= \begin{bmatrix} 9 \\ 6 \end{bmatrix} \end{align*}

So, w⃗∗\vec w^* is the solution to

[6333]⏟XTXw⃗∗=[96]⏟XTy⃗\underbrace{{\color{#3d81f6} \begin{bmatrix} 6 & 3 \\ 3 & 3 \end{bmatrix}}}_{\color{#3d81f6}X^TX}\vec w^* = \underbrace{\begin{bmatrix} 9 \\ 6 \end{bmatrix}}_{\color{#3d81f6}X^T\color{orange}\vec y}

and here, it’s pretty clear that w⃗∗=[11]\vec w^* = \begin{bmatrix} 1 \\ 1 \end{bmatrix}. (Again, we could go through the hassle of inverting XTXX^TX to get the same answer, but there’s no need to.)

This means that the projection of y⃗\color{orange} \vec y onto colsp(X)\text{colsp}({\color{#3d81f6} X}) is

Xw⃗∗=[1021110−1][11]=[132−1]=y⃗{\color{#3d81f6}X}\vec w^* = \begin{bmatrix} 1 & 0 \\ 2 & 1 \\ 1 & 1 \\ 0 & -1 \end{bmatrix}\begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 3 \\ 2 \\ -1 \end{bmatrix} = \color{orange}\vec y

as we predicted from the start.

Example: Orthogonality with the Columns of XX

Let X=[12320212]\color{#3d81f6} X = \begin{bmatrix} 1 & 2 \\ 3 & 2 \\ 0 & 2 \\ 1 & 2 \end{bmatrix}, Z=[100−10130]\color{#004d40} Z = \begin{bmatrix} 1 & 0 \\ 0 & -1 \\ 0 & 1 \\ 3 & 0 \end{bmatrix}, and y⃗\color{orange} \vec y be any vector in R4\mathbb{R}^4.

Let p⃗X\vec p_X and p⃗Z\vec p_Z be the orthogonal projections of y⃗\color{orange} \vec y onto colsp(X)\text{colsp}({\color{#3d81f6} X}) and colsp(Z)\text{colsp}({\color{#004d40} Z}), respectively.

Explain why it is guaranteed that the components of the vector

e⃗X=y⃗−p⃗X{\color{#d81a60} \vec e_X} = {\color{orange} \vec y} - \vec p_X

sum to zero, but the components of the vector e⃗Z=y⃗−p⃗Z{\color{#d81a60} \vec e_Z} = {\color{orange} \vec y} - \vec p_Z do not necessarily.

Solution

The error vectors e⃗X{\color{#d81a60} \vec e_X} and e⃗Z{\color{#d81a60} \vec e_Z} are orthogonal to the column spaces of their respective matrices. In other words,

XTe⃗X=0⃗,ZTe⃗Z=0⃗{\color{#3d81f6} X^T}{\color{#d81a60} \vec e_X}=\vec 0, \qquad {\color{#004d40} Z^T}{\color{#d81a60} \vec e_Z}=\vec 0

Suppose e⃗X\color{#d81a60} \vec e_X’s components are e1,e2,e3,e4e_1, e_2, e_3, e_4. In the case of X\color{#3d81f6} X, we have that

[13012222]⏟XT[e1e2e3e4]⏟e⃗X=[e1+3e2+e42e1+2e2+2e3+2e4]=[00]\underbrace{\begin{bmatrix} 1 & 3 & 0 & 1 \\ 2 & 2 & 2 & 2 \end{bmatrix}}_{\color{#3d81f6} X^T} \underbrace{\begin{bmatrix} e_1 \\ e_2 \\ e_3 \\ e_4 \end{bmatrix}}_{\color{#d81a60} \vec e_X} = \begin{bmatrix} e_1 + 3e_2 + e_4 \\ 2e_1 + 2e_2 + 2e_3 + 2e_4 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}

The second component of the above equation implies that 2e1+2e2+2e3+2e4=02e_1 + 2e_2 + 2e_3 + 2e_4 = 0, so e1+e2+e3+e4=0e_1 + e_2 + e_3 + e_4 = 0, meaning the components of e⃗X\color{#d81a60} \vec e_X sum to 0.

In the case of Z\color{#004d40} Z, we don’t have that same assurance.

[10030−110]⏟ZT[e1e2e3e4]⏟e⃗Z=[e1+3e4−e2+e3]=[00]\underbrace{\begin{bmatrix} 1 & 0 & 0 & 3 \\ 0 & -1 & 1 & 0\end{bmatrix}}_{\color{#004d40} Z^T} \underbrace{\begin{bmatrix} e_1 \\ e_2 \\ e_3 \\ e_4 \end{bmatrix}}_{\color{#d81a60} \vec e_Z} = \begin{bmatrix} e_1 + 3e_4 \\ -e_2 + e_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}

The condition ZTe⃗Z=[00]{\color{#004d40} Z^T}{\color{#d81a60} \vec e_Z} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} gives us some relationships involving e1,e2,e3,e4e_1, e_2, e_3, e_4, but not enough to guarantee that the components sum to 0.

Key takeaway: Remember that XTe⃗X=[00]\color{#3d81f6} X^T \color{#d81a60} \vec e_X = \begin{bmatrix} 0 \\ 0 \end{bmatrix} implies that e⃗X\color{#d81a60} \vec e_X is orthogonal to any linear combination of the columns of X\color{#3d81f6} X. So, if you can create a column of all 1’s using a linear combination of X\color{#3d81f6} X’s columns, then the components of e⃗X\color{#d81a60} \vec e_X will sum to 0, no matter which vector y⃗\color{orange} \vec y you choose to project onto colsp(X)\text{colsp}({\color{#3d81f6} X}).

This is one of the most important consequences of orthogonal projections, especially as it relates to linear regression.

Example: XX with Orthogonal Columns

Let X=[3/13−4/54/133/512/130]\color{#3d81f6} X = \begin{bmatrix} 3/13 & -4/5 \\ 4/13 & 3/5 \\ 12/13 & 0 \end{bmatrix} and y⃗=[101]\color{orange} \vec y = \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix}.

Find the value of w⃗∗\vec w^* that minimizes ∥y⃗−Xw⃗∥2\lVert {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w \rVert^2. What about X\color{#3d81f6} X makes this easier than in other examples?

Solution

Notice that both of X\color{#3d81f6} X’s columns are orthogonal to each other, and are unit vectors, i.e. they are orthonormal. This means that

XTX=[3/134/135/13−4/53/50][3/13−4/54/133/55/130]=[1001]\color{#3d81f6} X^T \color{#3d81f6} X = \begin{bmatrix} 3/13 & 4/13 & 5/13 \\ -4/5 & 3/5 & 0 \end{bmatrix} \begin{bmatrix} 3/13 & -4/5 \\ 4/13 & 3/5 \\ 5/13 & 0 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}

is the identity matrix, so

w⃗∗=(XTX)−1⏟IXTy⃗=XTy⃗=[15/13−4/5]\vec w^* = \underbrace{({\color{#3d81f6} X^T \color{#3d81f6} X})^{-1}}_{I} {\color{#3d81f6} X^T} {\color{orange} \vec y} = {\color{#3d81f6} X^T} {\color{orange} \vec y} = \begin{bmatrix} 15/13 \\ -4/5 \end{bmatrix}

Usually, our data doesn’t come to us with orthogonal columns. But, using the Gram-Schmidt process introduced to you in Homework 7, you can convert a set of linearly independent vectors into an orthogonal set with the same span, allowing you to leverage that QTQ=IQ^TQ = I and simplify the projection process.

Key takeaway: Orthonormal vectors are very easy to work with!

Example: Why Can’t We Separate?

Why can’t we separate

w⃗∗=(XTX)−1XTy⃗\vec w^* = ({\color{#3d81f6} X^T X})^{-1} {\color{#3d81f6} X^T} \color{orange} \vec y

into

w⃗∗=X−1(XT)−1XT⏟Iy⃗=X−1y⃗\vec w^* = {\color{#3d81f6} X^{-1}} \underbrace{({\color{#3d81f6} X^T})^{-1} {\color{#3d81f6} X^T}}_{I} {\color{orange} \vec y} = {\color{#3d81f6} X^{-1}} \color{orange} \vec y

in general?

Solution

In general, X\color{#3d81f6} X is not a square matrix, so it can’t be invertible!

If X\color{#3d81f6} X is square and invertible, then the steps above are valid. But, if X\color{#3d81f6} X is a square and invertible n×nn \times n matrix, then colsp(X)=Rn\text{colsp}({\color{#3d81f6} X}) = \mathbb{R}^n, meaning that any y⃗\color{orange} \vec y in Rn\mathbb{R}^n can be written as a linear combination of X\color{#3d81f6} X’s columns, meaning that y⃗\color{orange} \vec y is already in colsp(X)\text{colsp}({\color{#3d81f6} X}), and there is no projection error (or need to do a projection in the first place).

What if XX’s Columns are Linearly Dependent?

In the case where X\color{#3d81f6} X’s columns are linearly dependent, we can’t invert XTX\color{#3d81f6} X^T \color{#3d81f6} X to solve for w⃗∗\vec w^*. This means that

XTXw⃗=XTy⃗{\color{#3d81f6} X^T} {\color{#3d81f6} X} \vec w = {\color{#3d81f6} X^T} {\color{orange} \vec y}

has infinitely many solutions. Let’s give more thought to what these solutions actually are.

First, note that all of these solutions for w⃗∗\vec w^* correspond to the same projection, p⃗=Xw⃗∗\vec p = X \vec w^*. The “best approximation” of y⃗\color{orange} \vec y in colsp(X)\text{colsp}({\color{#3d81f6}X}) is always just one vector; if there are infinitely many w⃗∗\vec w^*'s, that just means there are infinitely many ways of describing that one best approximation. This is because if vectors are linearly dependent, then their linear combinations aren’t unique, but if they are linearly independent, their linear combinations are unique.

Let me drive this point home further. Let’s suppose both w⃗1\vec w_1 and w⃗2\vec w_2 satisfy

XTXw⃗=XTy⃗{\color{#3d81f6} X^T} {\color{#3d81f6} X} \vec w = {\color{#3d81f6} X^T} {\color{orange} \vec y}

Then,

XTXw⃗1−XTXw⃗2=y⃗−y⃗=0⃗X^TX \vec w_1 - X^TX \vec w_2 = \vec y - \vec y = \vec 0

which means that

(XTX)(w⃗1−w⃗2)=0⃗(X^TX)(\vec w_1 - \vec w_2) = \vec 0

i.e. the difference between the two vectors, w⃗1−w⃗2\vec w_1 - \vec w_2, is in nullsp(XTX)\text{nullsp}(X^TX). But, back in Chapter 2.8, we proved that XTXX^TX and XX have the same null space, meaning any vector that gets sent to 0⃗\vec 0 by XX also gets sent to 0⃗\vec 0 by XTXX^TX, and vice versa.

So,

X(w⃗1−w⃗2)=0⃗X(\vec w_1 - \vec w_2) = \vec 0

too, but that just means

Xw⃗1=Xw⃗2X \vec w_1 = X \vec w_2

meaning that even though w⃗1\vec w_1 and w⃗2\vec w_2 are different-looking coefficient vectors, they both still correspond to the same linear combination of XX’s columns!

Let’s see how we can apply this to an example. Let X=[310621310]{\color{#3d81f6} X = \begin{bmatrix} 3 & 1 & 0 \\ 6 & 2 & 1 \\ 3 & 1 & 0 \end{bmatrix}} and y⃗=[213]{\color{orange} \vec y = \begin{bmatrix} 2 \\ 1 \\ 3 \end{bmatrix}}. This is an example of a matrix with linearly dependent columns, so there’s no unique w⃗∗\vec w^* that satisfies the normal equations.

Finding One Solution

One way to find a possible vector w⃗∗\vec w^* is to solve the normal equations. XTX\color{#3d81f6} X^T \color{#3d81f6} X is not invertible, so we can’t solve for w⃗∗\vec w^* uniquely, but we can try and find a parameterized solution.

I’m going to take a different approach here: instead, let’s just toss out the linearly dependent columns of X\color{#3d81f6} X and solve for w⃗∗\vec w^* using the remaining columns. Then, w⃗∗\vec w^* for the full X\color{#3d81f6} X can use the same coefficients for the linearly independent columns, but 0s for the dependent ones. Removing the linearly dependent columns does not change colsp(X)\text{colsp}({\color{#3d81f6} X}) (i.e. the set of all linear combinations of X\color{#3d81f6} X’s columns), so the projection is the same.

The easy solution is to keep columns 2 and 3, since their numbers are smallest. So, for now, let’s say

X′=[102110],y⃗=[213]X' = \begin{bmatrix} 1 & 0 \\ 2 & 1 \\ 1 & 0 \end{bmatrix}, \quad \color{orange} \vec y = \begin{bmatrix} 2 \\ 1 \\ 3 \end{bmatrix}

Here, w⃗′=(X′TX′)−1X′Ty⃗=[5/2−4]\vec w' = (X'^T X')^{-1} X'^T {\color{orange} \vec y} = \begin{bmatrix} 5/2 \\ -4 \end{bmatrix}. I won’t bore you with the calculations; you can verify them yourself.

Now, one possible w⃗∗\vec w^* for the full X\color{#3d81f6} X is [05/2−4]\begin{bmatrix} 0 \\ 5/2 \\ -4\end{bmatrix}, which keeps the same coefficients on columns 2 and 3 as in w⃗′\vec w', but 0 for the column we didn’t use.

Finding All Solutions

As I mentioned above, if there are infinitely many solutions to the normal equation, then the difference between any two solutions is in nullsp(XTX)\text{nullsp}({\color{#3d81f6} X^TX}), which is also nullsp(X)\text{nullsp}({\color{#3d81f6} X}). Put another way, if w⃗s\vec w_s satisfies the normal equations, then so does w⃗s+n⃗\vec w_s + \vec n for any n⃗∈nullsp(XTX)\vec n \in \text{nullsp}({\color{#3d81f6} X^TX}).

XTXw⃗s=XTy⃗XTX(w⃗s+n⃗)=XTXw⃗s+XTXn⃗⏟0⃗, by definition of null space=XTy⃗+0⃗\begin{align*} {\color{#3d81f6} X^T X} \vec w_s &= {\color{#3d81f6} X^T} \color{orange} \vec y \\ {\color{#3d81f6} X^T X} (\vec w_s + \vec n) &= {\color{#3d81f6} X^T X} \vec w_s + \underbrace{{\color{#3d81f6} X^T X} \vec n}_{\vec 0, \text{ by definition of null space}} &= {\color{#3d81f6} X^T} {\color{orange} \vec y} + \vec 0 \\ \end{align*}

So, once we have one w⃗∗\vec w^*, to get the rest, just add any vector in nullsp(XTX)\text{nullsp}({\color{#3d81f6} X^TX}) or nullsp(X)\text{nullsp}({\color{#3d81f6} X}) (since those are the same spaces).

What is nullsp(X)\text{nullsp}({\color{#3d81f6} X})? It’s the set of vectors v⃗\vec v such that Xv⃗=0⃗{\color{#3d81f6} X} \vec v = \vec 0.

In our particular example,

X=[310621310]{\color{#3d81f6} X = \begin{bmatrix} 3 & 1 & 0 \\ 6 & 2 & 1 \\ 3 & 1 & 0 \end{bmatrix}}

we see that rank(X)=2\text{rank}({\color{#3d81f6} X}) = 2, so nullsp(X)\text{nullsp}({\color{#3d81f6} X}) has a dimension of 3−2=13 - 2 = 1 (by the rank-nullity theorem), so it’s going to be the span of a single vector. All we need to do now is find one vector in nullsp(X)\text{nullsp}({\color{#3d81f6} X}), and we will know that the null space is the set of scalar multiples of that vector.

X[1−30]=[310621310][1−30]=[000]{\color{#3d81f6} X} \begin{bmatrix} 1 \\ -3 \\ 0 \end{bmatrix} = {\color{#3d81f6}\begin{bmatrix} 3 & 1 & 0 \\ 6 & 2 & 1 \\ 3 & 1 & 0 \end{bmatrix}}\begin{bmatrix} 1 \\ -3 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}

So, since nullsp(X)=nullsp(XTX)=span({[1−30]})\text{nullsp}({\color{#3d81f6} X}) = \text{nullsp}({\color{#3d81f6} X^TX}) = \text{span}\left(\left\{\begin{bmatrix} 1 \\ -3 \\ 0 \end{bmatrix}\right\}\right), we know that the set of all possible w⃗∗\vec w^*'s is

[05/2−4]+t[1−30],t∈R⏟there are infinitely many solutions to the normal equations, but they’re all of this form\underbrace{\begin{bmatrix} 0 \\ 5/2 \\ -4 \end{bmatrix} + t\begin{bmatrix} 1 \\ -3 \\ 0 \end{bmatrix}, t \in \mathbb{R}}_\text{there are infinitely many solutions to the normal equations, but they're all of this form}

This is not a subspace, since it doesn’t contain the zero vector.


There’s another way to arrive at this set of possible w⃗∗\vec w^*'s: we can solve the normal equations directly. I wouldn’t recommend this second approach since it’s much longer, but I’ll add it here for completeness.

XTX=[363121010][310621310]=[541861862621]XTy⃗=[363121010][213]=[2171]\begin{align*} {\color{#3d81f6}X^TX}&={\color{#3d81f6}\begin{bmatrix} 3 & 6 & 3 \\ 1 & 2 & 1 \\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} 3 & 1 & 0 \\ 6 & 2 & 1 \\ 3 & 1 & 0 \end{bmatrix}}={\color{#3d81f6}\begin{bmatrix}54 & 18 & 6 \\ 18 & 6 & 2 \\ 6 & 2 & 1 \end{bmatrix}} \\ \\ {\color{#3d81f6}X^T}{\color{orange}\vec y}&={\color{#3d81f6}\begin{bmatrix} 3 & 6 & 3 \\ 1 & 2 & 1 \\ 0 & 1 & 0 \end{bmatrix}}{\color{orange}\begin{bmatrix} 2 \\ 1 \\ 3 \end{bmatrix}}=\begin{bmatrix} 21 \\ 7 \\ 1 \end{bmatrix} \end{align*}

Then, the normal equations XTXw⃗∗=XTy⃗{\color{#3d81f6}X^TX}\vec w^*={\color{#3d81f6}X^T}\color{orange}\vec y give us

54w1∗+18w2∗+6w3∗=2118w1∗+6w2∗+2w3∗=76w1∗+2w2∗+w3∗=1\begin{align*} 54w^*_1 + 18w^*_2 + 6w^*_3 &= 21 \\18w^*_1 + 6w^*_2 + 2w^*_3&=7 \\6w^*_1 + 2w^*_2 + w^*_3 &= 1 \end{align*}

The first and second equations are just scalar multiples of each other, so we can disregard one of them, and solve for a form where we can use one unknown as a parameter for the other two. To illustrate, let’s pick w1∗=tw_1^* = t.

18t+6w2∗+2w3∗=76t+2w2∗+w3∗=1\begin{align*} 18t + 6w^*_2 + 2w^*_3&=7 \tag{2} \\6t + 2w^*_2 + w^*_3 &= 1 \tag{3} \end{align*}

(2)−3⋅(3)(2) - 3 \cdot (3) gives us w3∗=−4w_3^* = -4. Plugging this into both equations gives us

18t+6w2∗−8=7  ⟹  18t+6w2∗=156t+2w2∗−4=1  ⟹  6t+2w2∗=5\begin{align*} 18t + 6w_2^* - 8 &= 7 \implies 18t + 6w_2^* &= 15 \\6t + 2w_2^* - 4 &= 1 \implies 6t + 2w_2^* &= 5 \end{align*}

These are now both the same equation; the first one is just 3 times the second. So, we can solve for w2∗w_2^* in terms of tt:

w2∗=5−6t2w_2^* = \frac{5-6t}{2}

which gives us the complete solution

w⃗∗=[t5−6t2−4],t∈R\vec w^* = \begin{bmatrix} t \\ \frac{5-6t}{2} \\ -4 \end{bmatrix}, t \in \mathbb{R}

This is the exact same line as using the null space approach! Plug in t=0t = 0 to get [05/2−4]\begin{bmatrix} 0 \\ 5/2 \\ -4 \end{bmatrix}, for example. Again, this is not a subspace, since it doesn’t contain the zero vector.

The Projection Matrix

So far, we’ve established that the vector in colsp(X)\text{colsp}({\color{#3d81f6} X}) that is closest to y⃗\color{orange} \vec y is the vector Xw⃗∗{\color{#3d81f6} X} \vec w^*, where w⃗∗\vec w^* is the solution to the normal equations,

XTXw⃗∗=XTy⃗{\color{#3d81f6} X^T X} \vec w^* = {\color{#3d81f6} X^T} \color{orange} \vec y

If XTXX^TX is invertible, then w⃗∗\vec w^* is the unique vector

w⃗∗=(XTX)−1XTy⃗\vec w^* = ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T \color{orange} \vec y

meaning that the vector in colsp(X)\text{colsp}({\color{#3d81f6} X}) that is closest to y⃗\color{orange} \vec y is

p⃗=Xw⃗∗=X(XTX)−1XTy⃗\vec p = {\color{#3d81f6} X} \vec w^* = {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T \color{orange} \vec y

You’ll notice that the above expression also looks like a linear transformation applied to y⃗\color{orange} \vec y, where y⃗\color{orange} \vec y is being multiplied by the matrix

P=X(XTX)−1XTP = {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T

The matrix PP is called the projection matrix. In other classes, it is called the “hat matrix”, because they might use w^\hat{\mathbf{w}} instead of w⃗∗\vec w^* and y^\hat{\mathbf{y}} instead of p⃗\vec p, and in that notation, y^=Py\hat{\mathbf{y}} = P \mathbf{y}, so PP puts a “hat” on y\mathbf{y}. (I don’t use hat notation in this class because drawing a hat on top of a vector is awkward. Doesn’t w⃗^\hat{\vec w} look strange?)

So,

p⃗=Xw⃗∗=Py⃗\vec p = {\color{#3d81f6} X} \vec w^* = P \color{orange} \vec y

shows us that there are two ways to interpret the act of projecting y⃗\color{orange} \vec y onto colsp(X)\text{colsp}({\color{#3d81f6} X}):

  1. The resulting vector is some optimal linear combination of X{\color{#3d81f6} X}'s columns.

  2. The resulting vector is the result of applying the linear transformation PP to y⃗\color{orange} \vec y.

Let’s work out an example. Suppose

X=[30015460],y⃗=[123]{\color{#3d81f6} X = \begin{bmatrix} {\color{#3d81f6} 3} & {\color{#3d81f6} 0} \\ {\color{#3d81f6} 0} & {\color{#3d81f6} 154} \\ {\color{#3d81f6} 6} & {\color{#3d81f6} 0} \end{bmatrix}}, \quad {\color{orange} \vec y = \begin{bmatrix} {\color{orange} 1} \\ {\color{orange} 2} \\ {\color{orange} 3} \end{bmatrix}}

X\color{#3d81f6} X’s columns are linearly independent, so XTX\color{#3d81f6} X^T X is invertible, and

P=X(XTX)−1XTP = {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T

is well-defined.

X = np.array([[3, 0], 
              [0, 154], 
              [6, 0]])

P = X @ np.linalg.inv(X.T @ X) @ X.T
P
array([[0.2, 0. , 0.4], [0. , 1. , 0. ], [0.4, 0. , 0.8]])
P @ np.array([1, 2, 3])
array([1.4, 2. , 2.8])

P=[0.200.40100.400.8]P = \begin{bmatrix} 0.2 & 0 & 0.4 \\ 0 & 1 & 0 \\ 0.4 & 0 & 0.8 \end{bmatrix} contains the information we need to project y⃗\color{orange} \vec y onto colsp(X)\text{colsp}({\color{#3d81f6} X}). Each row of PP tells us the right mixture of y⃗\color{orange} \vec y’s components we need to construct the projection.

Notice that PP’s second row is [010]T\begin{bmatrix} 0 & 1 & 0 \end{bmatrix}^T. This came from the fact that X\color{#3d81f6} X’s first column had a second component of 0 while its second column had a non-zero second component but zeros in the other two components, meaning that we can scale X\color{#3d81f6} X’s second column to exactly match y⃗\color{orange} \vec y’s second component. Change the 154 in X\color{#3d81f6} X to any other non-zero value and PP won’t change!

Additionally, if we consider some y⃗\color{orange} \vec y that is already in colsp(X)\text{colsp}({\color{#3d81f6} X}), then multiplying it by PP doesn’t change it! For example, if we set y⃗=[31546]\color{orange} \vec y = \begin{bmatrix} 3 \\ 154 \\ 6 \end{bmatrix} (the sum of X\color{#3d81f6} X’s columns), then Py⃗=y⃗P {\color{orange} \vec y} = {\color{orange} \vec y}.

P @ np.array([3, 154, 6])
array([ 3., 154., 6.])

Let’s work through some examples that develop our intuition for PP.

Example: Is PP invertible?

Suppose P=X(XTX)−1XTP = {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T exists, meaning XTX\color{#3d81f6} X^T X is invertible. Is PP invertible? If so, what is its inverse?

Solution

Before we do any calculations, intuitively the answer should be no. Once we’ve projected onto colsp(X)\text{colsp}({\color{#3d81f6} X}), we’ve lost information, since we went from an arbitrary vector in Rn\mathbb{R}^n to a vector in a smaller subspace, so it shouldn’t be possible to reverse the projection. Put another way, two different vectors in Rn\mathbb{R}^n might have the same “shadow” onto colsp(X)\text{colsp}({\color{#3d81f6} X}).

Even in the most recent example, P=[0.200.40100.400.8]P = \begin{bmatrix} 0.2 & 0 & 0.4 \\ 0 & 1 & 0 \\ 0.4 & 0 & 0.8 \end{bmatrix} is not invertible, since column 3 is a multiple of column 1.

Let’s think about this a bit more from the perspective of ranks. It turns out that rank(P)=rank(X)\text{rank}(P) = \text{rank}({\color{#3d81f6}X}); I’ve provided a proof of this at the bottom of the solutions box, but you might want to attempt it on your own for practice.

Remember that X\color{#3d81f6} X is an n×dn \times d matrix, meaning rank(X)≤min⁡(n,d)\text{rank}({\color{#3d81f6}X}) \leq \min(n, d). X\color{#3d81f6} X doesn’t need to have a rank of nn for XTX\color{#3d81f6} X^TX to be invertible; it just needs to have a rank of dd.

Since PP is an n×nn \times n matrix, in general it won’t be the case that rank(P)=n\text{rank}(P) = n. To illustrate, in the example above where P=[0.200.40100.400.8]P = \begin{bmatrix} 0.2 & 0 & 0.4 \\ 0 & 1 & 0 \\ 0.4 & 0 & 0.8 \end{bmatrix}, XX was a 3×23 \times 2 matrix with rank 2.

The only case in which rank(P)=n\text{rank}(P) = n is when rank(X)=n\text{rank}({\color{#3d81f6}X}) = n, which also only happens when X\color{#3d81f6} X is an n×nn \times n square matrix that is also invertible. In such a case, colsp(X)=Rn\text{colsp}({\color{#3d81f6}X}) = \mathbb{R}^n, and so we probably wouldn’t set out to project onto colsp(X)\text{colsp}({\color{#3d81f6} X}) in the first place, since any vector in Rn\mathbb{R}^n is already a linear combination of X{\color{#3d81f6} X}'s columns.


Extra: Proof that rank(P)=rank(X)\text{rank}(P) = \text{rank}({\color{#3d81f6} X})

We can show that this is the case by showing that PP and X\color{#3d81f6} X both have the same column spaces. This proof also helps explain why the normal equation, XTXw⃗=XTy⃗{\color{#3d81f6} X^TX} \vec w = {\color{#3d81f6}X^T} \color{orange} \vec y, always has at least one solution for w⃗\vec w, even when XTX\color{#3d81f6} X^TX isn’t invertible.

  1. Show colsp(P)⊆colsp(X)\text{colsp}(P) \subseteq \text{colsp}({\color{#3d81f6} X}), i.e. show that any vector in the column space of PP is also in the column space of X\color{#3d81f6} X.

    If v⃗∈colsp(P)\vec v \in \text{colsp}(P), then v⃗\vec v can be written as a linear combination of PP’s columns. Say v⃗=Pu⃗\vec v = P \vec u for some u⃗∈Rn\vec u \in \mathbb{R}^n. Then,

    v⃗=Pu⃗=X((XTX)−1XTu⃗)⏟some vector\vec v = P \vec u = {\color{#3d81f6} X} \underbrace{\left( ({\color{#3d81f6} X^T X})^{-1} {\color{#3d81f6} X^T} \vec u \right)}_{\text{some vector}}

    Here, we see that if v⃗=Pu⃗\vec v = P \vec u, then v⃗\vec v is also in X{\color{#3d81f6} X}'s column space. So,

    v⃗ is a lin. comb. of P’s columns  ⟹  v⃗ is a lin. comb. of X’s columns\vec v \text{ is a lin. comb. of $P$'s columns} \implies \vec v \text{ is a lin. comb. of ${\color{#3d81f6} X}$'s columns}

  2. Show colsp(X)⊆colsp(P)\text{colsp}({\color{#3d81f6} X}) \subseteq \text{colsp}(P), i.e. show that any vector in the column space of X\color{#3d81f6} X is also in the column space of PP.

    This direction is a bit more involved. Let’s start by considering some vector v⃗∈colsp(X)\vec v \in \text{colsp}({\color{#3d81f6} X}), meaning

    v⃗=Xu⃗\vec v = {\color{#3d81f6} X} \vec u

    for some u⃗∈Rd\vec u \in \mathbb{R}^d. What happens if we multiply both sides by PP?

    v⃗=Xu⃗Pv⃗=PXu⃗Pv⃗=X(XTX)−1XT⏟PXu⃗Pv⃗=X(XTX)−1XTX⏟Iu⃗Pv⃗=Xu⃗=v⃗⏟original definition of v⃗\begin{align*} \vec v &= {\color{#3d81f6} X} \vec u \\ P \vec v &= P {\color{#3d81f6} X} \vec u \\ P \vec v &= \underbrace{{\color{#3d81f6} X}({\color{#3d81f6} X^T X})^{-1} {\color{#3d81f6} X^T}}_P {\color{#3d81f6}X} \vec u \\ P \vec v &= {\color{#3d81f6} X}\underbrace{({\color{#3d81f6} X^T X})^{-1} {\color{#3d81f6} X^T} {\color{#3d81f6}X}}_I \vec u \\ P \vec v &= \underbrace{{\color{#3d81f6} X} \vec u = \vec v}_{\text{original definition of $\vec v$}} \end{align*}

    If v⃗∈colsp(X)\vec v \in \text{colsp}({\color{#3d81f6} X}), then v⃗\vec v can be written as a linear combination of X{\color{#3d81f6} X}'s columns. Say v⃗=Xw⃗\vec v = {\color{#3d81f6} X} \vec w for some w⃗∈Rd\vec w \in \mathbb{R}^d. Then,

    v⃗=Xw⃗=Pw⃗\vec v = {\color{#3d81f6} X} \vec w = P \vec w

    So, if v⃗=Xu⃗\vec v = {\color{#3d81f6} X} \vec u, then v⃗=Pv⃗\vec v = P \vec v, meaning that v⃗\vec v is also in PP’s column space if it’s in X{\color{#3d81f6} X}'s column space. Intuitively, v⃗=Pv⃗\vec v = P \vec v means that if v⃗\vec v is already in the span of X\color{#3d81f6} X’s columns, then projecting it onto colsp(X)\text{colsp}({\color{#3d81f6} X}) doesn’t change it.

Now that we’ve shown that colsp(P)⊆colsp(X)\text{colsp}(P) \subseteq \text{colsp}({\color{#3d81f6} X}) and colsp(X)⊆colsp(P)\text{colsp}({\color{#3d81f6} X}) \subseteq \text{colsp}(P), we can conclude that colsp(P)=colsp(X)\text{colsp}(P) = \text{colsp}({\color{#3d81f6} X}). If two sets are subsets of each other, then they must be equal.

Example: Is PP orthogonal?

Is PP orthogonal?

Solution

No. Orthogonal matrices QQ have the property that QTQ=QQT=IQ^TQ = QQ^T = I, meaning that

QT=Q−1Q^T = Q^{-1}

But, as we saw, PP is not invertible in general, so it can’t satisfy this property. This tells us that PP does not perform a rotation; projections are not rotations. Rotations can be undone but projections can’t.

Example: Is PP symmetric?

Is PP symmetric?

Solution

Yes. Symmetric matrices AA have the property that AT=AA^T = A. We can show that PP satisfies this property; to do so, we’ll need to use the fact that (AB)T=BTAT(AB)^T = B^T A^T.

PT=(X(XTX)−1XT)T=(XT)T((XTX)−1)TXT=X(XTX)−1XT=P\begin{align*} P^T &= \left( {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} {\color{#3d81f6} X^T} \right)^T \\ &= \left( {\color{#3d81f6} X^T} \right)^T \left( ({\color{#3d81f6} X^T X})^{-1} \right)^T {\color{#3d81f6} X^T} \\ &= {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T \\ &= P \end{align*}

Going from the second to the third line, we used the fact that XTX\color{#3d81f6} X^TX is symmetric, and so is its inverse. Remember that XTX\color{#3d81f6} X^T X is a square matrix consisting of the dot products of the columns of X\color{#3d81f6} X with themselves.

Example: Is PP idempotent?

Recall, an idemponent matrix AA satisfies A2=AA^2 = A. Is PP idempotent?

Solution

Yes.

P2=P⋅P=(X(XTX)−1XT)X(XTX)−1XT=X((XTX)−1XTX)⏟I(XTX)−1XT=X(XTX)−1XT=P\begin{align*} P^2 &= P \cdot P \\ &= \left({\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} {\color{#3d81f6} X^T} \right) {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T \\ &= {\color{#3d81f6} X} \underbrace{\left( ({\color{#3d81f6} X^T X})^{-1} {\color{#3d81f6} X^T} {\color{#3d81f6} X} \right)}_{I} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T \\ &= {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T \\ &= P \end{align*}

Intuitively, this means that P2y⃗P^2 {\color{orange} \vec y} is the same as Py⃗P {\color{orange} \vec y}, meaning that once we’ve projected y⃗\vec y onto colsp(X)\text{colsp}({\color{#3d81f6} X}), projecting its projection p⃗\vec p again onto colsp(X)\text{colsp}({\color{#3d81f6} X}) gives us back the same p⃗\vec p, since p⃗\vec p is already in colsp(X)\text{colsp}({\color{#3d81f6} X}).

Example: What is PXPX, and why?

What is PXPX? What does the result mean?

Solution
PX=X(XTX)−1XTX=X(XTX)−1XTX=XI=X\begin{align*} P{\color{#3d81f6} X} &= {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T {\color{#3d81f6} X} \\ &= {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T {\color{#3d81f6} X} \\ &= {\color{#3d81f6} X} I \\ &= {\color{#3d81f6} X} \end{align*}

Interpret PXP {\color{#3d81f6} X} as a matrix made up of Px⃗(1)P \color{#3d81f6} \vec x^{(1)}, Px⃗(2)P \color{#3d81f6} \vec x^{(2)}, ..., Px⃗(d)P \color{#3d81f6} \vec x^{(d)} as its columns. Px⃗(i)P \color{#3d81f6} \vec x^{(i)} is the projection of x⃗(i)\color{#3d81f6} \vec x^{(i)} onto colsp(X)\text{colsp}({\color{#3d81f6} X}), but since x⃗(i)\color{#3d81f6} \vec x^{(i)} is already in colsp(X)\text{colsp}({\color{#3d81f6} X}), projecting it again onto colsp(X)\text{colsp}({\color{#3d81f6} X}) gives us back the same x⃗(i)\color{#3d81f6} \vec x^{(i)}. So, PXP {\color{#3d81f6} X} should just be X\color{#3d81f6} X again.

Example: Rotations, Reflections, and Projections

Suppose AA is an arbitrary n×dn \times d matrix. Describe the conditions on AA that make the corresponding linear transformation f(x⃗)=Ax⃗f(\vec x) = A \vec x a...

  1. Rotation

  2. Reflection

  3. Projection

Summary

Let’s take a step back and walk through our logic from this section once more, since it’s that important.

Suppose X{\color{#3d81f6} X} is an n×dn \times d matrix and y⃗\color{orange} \vec y is some vector in Rn\mathbb{R}^n.

Orthogonal Projections

  1. Our goal is to find the linear combination of X\color{#3d81f6} X’s columns that is closest to y⃗\color{orange} \vec y.

  2. This boils down to finding the vector w⃗\vec w that minimizes ∥y⃗−Xw⃗∥2\lVert {\color{orange} \vec y} - {\color{#3d81f6}X } \vec w \rVert^2.

  3. The vector w⃗∗\vec w^* that minimizes ∥y⃗−Xw⃗∥2\lVert {\color{orange} \vec y} - {\color{#3d81f6}X } \vec w \rVert^2 makes the resulting error vector, e⃗=y⃗−Xw⃗∗{\color{#d81a60} \vec e} = {\color{orange} \vec y} - {\color{#3d81f6} X} \vec w^*, orthogonal to the columns of X\color{#3d81f6} X.

  4. The w⃗∗\vec w^* that makes the error vector orthogonal to the columns of X\color{#3d81f6} X is the one that satisfies the normal equation,

    XTXw⃗∗=XTy⃗{\color{#3d81f6} X^T} {\color{#3d81f6} X} \vec w^* = {\color{#3d81f6} X^T} {\color{orange} \vec y}

  5. If XTX{\color{#3d81f6} X^T} {\color{#3d81f6} X} is invertible, which happens if and only if X\color{#3d81f6}X’s columns are linearly independent, then w⃗∗\vec w^* is the unique vector

    w⃗∗=(XTX)−1XTy⃗\vec w^* = ({\color{#3d81f6} X^T} {\color{#3d81f6} X})^{-1} {\color{#3d81f6} X^T} \color{orange} \vec y

    otherwise, there are infinitely many solutions to the normal equation. All of these infinitely many solutions correspond to the same projection, p⃗=Xw⃗∗\vec p = {\color{#3d81f6} X} \vec w^*. If w⃗′\vec w' is one solution (which can be found by removing the linearly dependent columns of X\color{#3d81f6} X), then all other solutions are of the form w⃗′+n⃗\vec w' + \vec n, where n⃗\vec n is any vector in nullsp(X)=nullsp(XTX)\text{nullsp}({\color{#3d81f6} X}) = \text{nullsp}({\color{#3d81f6} X^TX}).

The Projection Matrix

Assuming X\color{#3d81f6} X has linearly independent columns, the projection matrix is

P=X(XTX)−1XTP = {\color{#3d81f6} X} ({\color{#3d81f6} X^T X})^{-1} \color{#3d81f6} X^T

PP is defined such that Py⃗P \color{orange} \vec y is the vector in colsp(X)\text{colsp}({\color{#3d81f6} X}) that is closest to y⃗\color{orange} \vec y. PP is symmetric and idemponent, but not invertible nor orthogonal.

We’re now finally ready to head back to the land of machine learning.

EECS 245 Course Notes
2.9. Inverses
EECS 245 Course Notes
Chapter 3: Regression using Linear Algebra