6.4. The Complete Solution to the Normal Equations
In Chapter 6.3, we explored the idea of projecting a vector onto the column space of a matrix. If X has linearly independent columns, then the vector in colsp(X) that is closest to yβ is the vector Xwβ, where wβ is the unique solution to the normal equation,
XTXwβ=XTyβ
When X has linearly independent columns, then XTX is invertible, so we can solve for wβ uniquely.
In the case where Xβs columns are linearly dependent, we canβt invert XTX to solve for wβ. This means that
XTXw=XTyβ
has infinitely many solutions. Letβs give more thought to what these solutions actually are.
First, note that all of these solutions for wβcorrespond to the same projection,pβ=Xwβ. The βbest approximationβ of yβ in colsp(X) is always just one vector; if there are infinitely many wβ's, that just means there are infinitely many ways of describing that one best approximation. Remember, if vectors are linearly independent, then any of their linear combinations can only be expressed in one way; if they are linearly dependent, then their linear combinations can be expressed in infinitely many ways.
In other words, if X has linearly dependent columns, then there are infinitely many wβ's that satisfy the normal equation, but they all correspond to the same projection pβ=Xwβ in the figure below.
from pathlib import Path
import sys
_notes_root = next(p for p in [Path.cwd(), *Path.cwd().parents] if (p / 'myst.yml').exists())
if str(_notes_root) not in sys.path:
sys.path.insert(0, str(_notes_root))
from plane_style import (BLUE, ORANGE, PINK, GRAY, plane_figure, add_plane,
vector3, vector_label, line3, label3,
style_surface_planes, projection_figure)
import numpy as np
import plotly.graph_objects as go
y = np.array([1,3,2]) * 1.5
v1 = np.array([1,0,0.25])
v2 = np.array([1,1,0])
X = np.column_stack([v1,v2])
w = np.linalg.solve(X.T @ X, X.T @ y)
p = X @ w
fig = projection_figure(y,X,[p],notation='*')
fig.show(renderer='plotly_mimetype')
Loading...
Let me drive this point home further. Letβs suppose both w1β and w2β satisfy
XTXw=XTyβ
Then,
XTXw1ββXTXw2β=yββyβ=0
which means that
(XTX)(w1ββw2β)=0
i.e. the difference between the two vectors, w1ββw2β, is in nullsp(XTX). But, back in Chapter 5.3, we proved that XTX and X have the same null space, meaning any vector that gets sent to 0 by X also gets sent to 0 by XTX, and vice versa.
So,
X(w1ββw2β)=0
too, but that just means
Xw1β=Xw2β
meaning that even though w1β and w2β are different-looking coefficient vectors, they both still correspond to the same linear combination of Xβs columns!
Letβs see how we can apply this to an example. Let X=β£β‘β363β121β010ββ¦β€β and yβ=β£β‘β213ββ¦β€β. This is an example of a matrix with linearly dependent columns, so thereβs no unique wβ that satisfies the normal equations.
One way to find a possible vector wβ is to solve the normal equations. XTX is not invertible, so we canβt solve for wβ uniquely, but we can still try and find a solution.
Hereβs one approach: letβs just toss out the linearly dependent columns of X and solve for wβ using the remaining columns. Then, wβ for the full 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) (i.e. the set of all linear combinations of 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
Here, wβ²=(Xβ²TXβ²)β1Xβ²Tyβ=[5/2β4β]. I wonβt bore you with the calculations; you can verify them yourself.
Now, one possible wβ for the full X is β£β‘β05/2β4ββ¦β€β, which keeps the same coefficients on columns 2 and 3 as in wβ², but 0 for the column we didnβt use.
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), which is also nullsp(X). Put another way, if wsβ satisfies the normal equations, then so does wsβ+n for any nβnullsp(XTX).
So, once we have one wβ, to get the rest, just add any vector in nullsp(XTX) or nullsp(X) (since those are the same subspaces).
What is nullsp(X)? Itβs the set of vectors v such that Xv=0.
In our particular example,
X=β£β‘β363β121β010ββ¦β€β
we see that rank(X)=2, so nullsp(X) has a dimension of 3β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), and we will know that the null space is the set of scalar multiples of that vector.
Since column 1 is three times column 2, the vector v=β£β‘β1β30ββ¦β€β must be in nullsp(X).
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β'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.
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ββ=t.
These are now both the same equation; the first one is just 3 times the second. So, we can solve for w2ββ in terms of t:
w2ββ=25β6tβ
which gives us the complete solution
wβ=β£β‘βt25β6tββ4ββ¦β€β,tβR
This is the exact same line as using the null space approach! Plug in t=0 to get β£β‘β05/2β4ββ¦β€β, for example. The set of all possible wβ's is not a subspace, since it doesnβt contain the zero vector.
So far, weβve established that the vector in colsp(X) that is closest to yβ is the vector Xwβ, where wβ is the solution to the normal equations,
XTXwβ=XTyβ
If XTX is invertible, then wβ is the unique vector
wβ=(XTX)β1XTyβ
meaning that the vector in colsp(X) that is closest to yβ is
pβ=Xwβ=X(XTX)β1XTyβ
Youβll notice that the above expression also looks like a linear transformation applied to yβ, where yβ is being multiplied by the matrix
P=X(XTX)β1XT
The matrix P is called the projection matrix. In other classes, it is called the βhat matrixβ, because they might use w^ instead of wβ and y^β instead of pβ, and in that notation, y^β=Py, so P puts a βhatβ on y. (I donβt use hat notation in this class because drawing a hat on top of a vector is awkward. Doesnβt w^ look strange?)
So,
pβ=Xwβ=Pyβ
shows us that there are two ways to interpret the act of projecting yβ onto colsp(X):
The resulting vector is some optimal linear combination of X's columns.
The resulting vector is the result of applying the linear transformation P to yβ.
P=β£β‘β0.200.4β010β0.400.8ββ¦β€β contains the information we need to project yβ onto colsp(X). Each row of P tells us the right mixture of yββs components we need to construct the projection.
Notice that Pβs second row is [0β1β0β]T. This came from the fact that 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βs second column to exactly match yββs second component. Change the 154 in X to any other non-zero value and P wonβt change!
Additionally, if we consider some yβ that is already in colsp(X), then multiplying it by P doesnβt change it! For example, if we set yβ=β£β‘β31546ββ¦β€β (the sum of Xβs columns), then Pyβ=yβ.
P @ np.array([3, 154, 6])
array([ 3., 154., 6.])
Letβs work through some examples that develop our intuition for P.
No. Orthogonal matrices Q have the property that QTQ=QQT=I, meaning that
QT=Qβ1
But, as we saw, P is not invertible in general, so it canβt satisfy this property. This tells us that Pdoes not perform a rotation; projections are not rotations. Rotations can be undone but projections canβt.
Yes. Symmetric matrices A have the property that AT=A. We can show that P satisfies this property; to do so, weβll need to use the fact that (AB)T=BTAT.
Going from the second to the third line, we used the fact that XTX is symmetric, and so is its inverse. Remember that XTX is a square matrix consisting of the dot products of the columns of X with themselves.
Intuitively, this means that P2yβ is the same as Pyβ, meaning that once weβve projected yβ onto colsp(X), projecting its projection pβ again onto colsp(X) gives us back the same pβ, since pβ is already in colsp(X).
Interpret PX as a matrix made up of Px(1), Px(2), ..., Px(d) as its columns. Px(i) is the projection of x(i) onto colsp(X), but since x(i) is already in colsp(X), projecting it again onto colsp(X) gives us back the same x(i). So, PX should just be X again.
Example: Rotations, Reflections, and ProjectionsΒΆ
Suppose A is an arbitrary nΓd matrix. Describe the conditions on A that make the corresponding linear transformation f(x)=Ax a...
Rotation
Reflection
Projection
Solution
Rotation: A should be orthogonal, meaning that
ATA=AAT=I
All orthogonal matrices have a determinant of 1 or -1, since det(ATA)=det(I) implies that det(A)2=1 (remember that det(AT)=det(A)). If we want f(x)=Ax to be a rotation rather than a reflection, we also need
det(A)=1
To summarize: a rotation matrix is orthogonal and has determinant 1. Implicit here is the fact that A is square, meaning that n=d. It doesnβt make sense to talk about a rotation that also changes the dimension of the space.
Reflection: The distinction between a reflection and a rotation is subtle in higher dimensions. Weβve mostly thought of βorthogonalβ and βrotationβ as being the same thing. But, for example,
A=[01β10β]
is orthogonal, but doesnβt rotate vectors β it reflects them across the line x2β=x1β in R2. The determinant of A is -1.
In general, if A is orthogonal and has a determinant of -1, then it either involves a reflection or a reflection and a rotation.
for some unit vector u. This matrix reflects vectors across the hyperplane uTx=0 (remmeber, in R2 a line is a hyperplane; in R3 a plane is a hyperplane). If A is of that form, then it indeed satisfies the properties above: it is orthogonal and has determinant -1. (To show that its determinant is -1, we need knowledge of eigenvalues, which weβll introduce in Chapter 9.) It is furthermore symmetric, since
AT=(Iβ2uuT)T=Iβ2uuT=A
To summarize: a reflection matrix is orthogonal,symmetric, and has determinant -1. But, just because a matrix is orthogonal and has determinant -1, it doesnβt have to be a reflection matrix: it could represent a reflection followed by a rotation.
Projection: To reason about the properties of a projection matrix, letβs think about the matrix P from earlier in this section. Suppose we have a subspace, for which the columns of X are a basis (in other words, suppose Xβs columns are linearly independent, and the subspace in question is colsp(X)). Then, the matrix P that projects vectors onto that subspace is given by
P=X(XTX)β1XT
As we saw above, P is both idempotent, meaning P2=P, and symmetric, meaning PT=P.
To summarize: a projection matrix is idempotent and symmetric.
There are matrices β like [10β10β] that are idempotent but not symmetric, and you can think of them as representing non-orthogonal projections. This is not a kind of projection weβve studied yet. Letβs explore. Suppose yββRn and P projects vectors onto some subspace. The error vector is yββPyβ. As we know for orthogonal projections, yββPyβ is orthogonal to anything in the column space of P, e.g. Pyβ. But it turns out the symmetry of P is what guarantees that this happens. Letβs take the dot product of yββPyβ with Pyβ:
In order for this dot product to be zero, we need PβPTP=0nΓnβ. To make this happen, we need PT=P, since this would make PTP=P2, which would turn PβPTP into PβP2, which is guaranteed by the idempotency of P (P2=P).
What this all means is that if PT=P, then yββPyβ is guaranteed to be orthogonal to Pyβ. But if PTξ =P, then PβPTP is not necessarily zero, meaning that yββPyβ is not necessarily orthogonal to Pyβ in general.
So in short, P being idempotent makes it some kind of projection, but the symmetry of P is what makes it an orthogonal projection. But for our purposes, weβve only considered orthogonal projections, so the term βprojectionβ without any other context means βorthogonal projectionβ.
To summarize:
Transformation
Properties
Rotation
Orthogonal (ATA=I), determinant = 1
Reflection
Orthogonal (ATA=I), symmetric (AT=A) if strict reflection, determinant = -1
Our goal is to find the linear combination of Xβs columns that is closest to yβ.
This boils down to finding the vector w that minimizes β₯yββXwβ₯2.
The vector wβ that minimizes β₯yββXwβ₯2 makes the resulting error vector,
e=yββXwβ
orthogonal to the columns of X.
The wβ that makes the error vector orthogonal to the columns of X is the one that satisfies the normal equation,
XTXwβ=XTyβ
If XTX is invertible, which happens if and only if Xβs columns are linearly independent, then wβ is the unique vector
wβ=(XTX)β1XTyβ
Otherwise, there are infinitely many solutions to the normal equation. All of these infinitely many solutions correspond to the same projection, pβ=Xwβ. If wβ² is one solution (which can be found by removing the linearly dependent columns of X), then all other solutions are of the form wβ²+n, where n is any vector in nullsp(X)=nullsp(XTX).