2.10 Least squares and regression as projection

When Ax = b has no solution, the best we can do is the point of the column space closest to b. The error of that best choice is perpendicular to every column, and that single fact gives the normal equations, the projection matrix, and the line of best fit.

Spagnolie opened Lecture 1 with three equations in two unknowns that cannot all be true at once, and said “we’ll have to decide how to make sense of this”. Real data always looks like that: more measurements than unknowns, each measurement a little off. This chapter makes sense of it. When no exact answer exists, we ask for the closest answer, and the closest answer turns out to be a projection, which Lecture 6 already taught us how to find.

Before you start

2.6: systems of linear equations and the column space. 2.5: matrix multiplication and the transpose. 2.4: inverses. 2.9: the inner product, orthogonal vectors, and projection. 2.3: null space. 2.2: linear independence. In this chapter vectors are columns of real numbers, and the inner product of $u, v$ in $\mathbb R^m$ is the dot product $u \cdot v = u^{\mathsf T} v = \sum_i u_i v_i$ .

Quick check

A

is a

3 \times 2

matrix and

x

is in

\mathbb R^2

. Where does

Ax

live?

Quick checkTwo vectors

u, v

satisfy

u \cdot v = 0

. What does that say?

Quick checkThe columns of

A

are linearly independent. What is the null space of

A

The picture

Four data points sit in the plane. A line $y = c + mt$ cannot pass through all four, so we ask for the line whose vertical misses are smallest overall. Each miss is drawn as a segment, and its square as a small square. The best line is the one for which the total area of the squares is as small as it can be. Drag the points and watch it re-fit. Then tick “try your own line” and try to beat it. You cannot.

That is the picture in the data plane. There is a second picture, in a different space, and it is the one that makes the problem solvable. Think of the four heights $y_1, \dots, y_4$ as ONE vector $b$ in $\mathbb R^4$ . The heights that ANY line could produce at those four $t$ values also form vectors in $\mathbb R^4$ , and they fill out a plane (a 2-dimensional subspace). So fitting a line means finding the point of that plane closest to $b$ . Here is the same idea one dimension down, where the “plane” is a line through the origin: the green arrow is the closest point, and the red dashed error meets the line at a right angle.

The idea, in words

Write Spagnolie’s three equations as one matrix equation:

\begin{aligned} 3x_1 + x_2 &= 1 \\ 2x_1 - x_2 &= 0 \\ 4x_1 + 7x_2 &= 2 \end{aligned} \qquad\text{that is}\qquad \underbrace{\begin{bmatrix} 3 & 1 \\ 2 & -1 \\ 4 & 7 \end{bmatrix}}_{A} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \underbrace{\begin{bmatrix} 1 \\ 0 \\ 2 \end{bmatrix}}_{b}.

The first two equations force $x_1 = 1/5$ and $x_2 = 2/5$ . The third then reads $4/5 + 14/5 = 18/5$ , not $2$ . So no $x$ works.

Chapter 2.6 told us why: $Ax = x_1 a_1 + x_2 a_2$ , where $a_1 = (3, 2, 4)$ and $a_2 = (1, -1, 7)$ are the columns of $A$ . So $Ax = b$ has a solution exactly when $b$ is in the column space of $A$ , the set of all combinations of the columns. Here the column space is a plane in $\mathbb R^3$ , and $b = (1, 0, 2)$ is not on it.

If we cannot hit $b$ , we can at least get as close as possible. Every vector $Ax$ lies in the plane, so we want the $x$ for which $Ax$ is nearest to $b$ . The leftover $b - Ax$ is the residual (the amount by which we miss). Its length is the miss distance. Making the squared length small means making the sum of the squares of its entries small. That sum of squares is where the name comes from.

Why· Why squares, and why closest means perpendicular

There are three reasons to use the squared length, and they agree.

Geometry. The closest point of a plane to $b$ is the foot of the perpendicular dropped from $b$ . Lecture 6 proved this with Pythagoras. Pythagoras speaks in squares, so squares are the natural way to measure.
Calculus. The squared length $\lVert Ax - b \rVert^2$ is a smooth function of $x$ whose gradient is $2A^{\mathsf T}(Ax - b)$ . Setting the gradient to zero gives exactly the equations we derive below by geometry. Spagnolie’s bonus note in HW4 says the same: the calculus route and the projection route give one answer, and he calls projection “by far the most elegant approach”.
Statistics (B). If each measurement is the true value plus independent Gaussian noise, the most likely $x$ is the one minimizing the sum of squares. This is why least squares is the default tool of experimental science.

Where this came from· Ceres, Gauss and Legendre

On 1 January 1801 Giuseppe Piazzi found a new body, Ceres, and followed it for about forty days before it was lost in the Sun’s glare. Astronomers had a short, noisy arc of observations and needed to say where to point their telescopes months later. Gauss, aged 24, computed an orbit from the observations (he later said he had used least squares in the calculation), and Ceres was seen again near his prediction at the end of 1801. Legendre published the method in 1805, in a work on the orbits of comets, as a way to use more measurements than unknowns. In 1809 Gauss claimed he had used it since 1795, which began a priority dispute, and in the same book he tied the method to the normal distribution of errors. The problem (more observations than unknowns, each a little off) is exactly Spagnolie’s three equations in two unknowns.

The plan for the chapter: (1) show that the closest point is the one whose error is perpendicular to the plane; (2) turn “perpendicular to every column” into the matrix equation $A^{\mathsf T}A\,x = A^{\mathsf T}b$ ; (3) find when that equation has exactly one solution; (4) package the answer as a matrix $P$ that drops any $b$ onto the plane; (5) see how computers do it without ever forming $A^{\mathsf T}A$ .

Definitions and theorems

Throughout, $A$ is an $m \times n$ matrix of real numbers with columns $a_1, \dots, a_n$ , $b$ is in $\mathbb R^m$ , and $\lVert v \rVert = \sqrt{v \cdot v}$ . Two facts about $A$ are used all the time.

$Ax = x_1 a_1 + \dots + x_n a_n$ , so the column space of $A$ is exactly $\{Ax : x \in \mathbb R^n\}$ .
The $j$ -th entry of $A^{\mathsf T} r$ is $a_j \cdot r$ : row $j$ of $A^{\mathsf T}$ is column $j$ of $A$ . So $A^{\mathsf T} r = O$ says that $r$ is perpendicular to every column of $A$ .

Definition B

Least squares solution

A vector $\hat x \in \mathbb R^n$ is a least squares solution of $Ax = b$ if

\lVert b - A\hat x \rVert \le \lVert b - Ax \rVert \quad \text{for every } x \in \mathbb R^n.

The vector $r = b - A\hat x$ is the residual. Its squared length $\lVert r \rVert^2 = r_1^2 + \dots + r_m^2$ is the sum of the squared residuals. If $Ax = b$ has an exact solution, it is a least squares solution with residual $O$ .

In one lineA vector $\hat x$ that makes $\lVert b - A\hat x \rVert$ as small as possible, even when $Ax = b$ has no solution.

Why it existsReal systems are overdetermined and slightly inconsistent. This says what answer to give anyway.

The first theorem is Lecture 6’s best approximation theorem, stated for the situation we need and proved again so the chapter stands on its own.

Theorem T Spagnolie Lecture 6; Apostol Thm 1.16

Perpendicular error means closest point

Let $W$ be a subspace of $\mathbb R^m$ , let $b \in \mathbb R^m$ , and let $p \in W$ with $(b - p) \cdot w = 0$ for every $w \in W$ . Then for every $w \in W$ ,

\lVert b - w \rVert^2 = \lVert b - p \rVert^2 + \lVert p - w \rVert^2 .

So $\lVert b - w \rVert \ge \lVert b - p \rVert$ , with equality only when $w = p$ .

In one lineIf $p$ is in a subspace $W$ and $b - p$ is perpendicular to all of $W$ , then $p$ is the closest point of $W$ to $b$ , and no other point of $W$ is as close.

Why it is trueAny other point $w$ of $W$ forms a right triangle with $b$ and $p$ , and the hypotenuse is longer than a leg.

ProofSpagnolie Lecture 6, page 2

The idea: Split the vector b - w into the part b - p, which is perpendicular to W, and the part p - w, which lies in W. Then Pythagoras applies.

b - w = (b - p) + (p - w)

Add and subtract

p

. This is the move from Lecture 6: the vector from

w

up to

b

is the vector from

w

p

followed by the vector from

p

b

p - w \in W

W

is a subspace, so it is closed under subtraction. Both

p

and

w

are in

W

(b - p) \cdot (p - w) = 0

The hypothesis:

b - p

is perpendicular to every vector of

W

, and

p - w

is one of them.

\lVert b - w \rVert^2 = \lVert b - p \rVert^2 + \lVert p - w \rVert^2

Expand

\lVert u + v \rVert^2 = (u + v)\cdot(u + v) = \lVert u \rVert^2 + 2\,u \cdot v + \lVert v \rVert^2

with

u = b - p

and

v = p - w

. The middle term is the zero from the line above. This is Pythagoras.

\lVert b - w \rVert^2 \ge \lVert b - p \rVert^2

, with equality only when

w = p

A square is never negative, and

\lVert p - w \rVert^2 = 0

only when

p - w = O

Proof workshopGive the reason for each line: the perpendicular error is the shortest

Each line of the proof is shown in order. Choose the reason that justifies it, then press Check.

$b - w = (b - p) + (p - w)$
$p - w \in W$
$(b - p) \cdot (p - w) = 0$
$\lVert b - w \rVert^2 = \lVert b - p \rVert^2 + 2\,(b - p)\cdot(p - w) + \lVert p - w \rVert^2 = \lVert b - p \rVert^2 + \lVert p - w \rVert^2$
$\lVert b - w \rVert^2 \ge \lVert b - p \rVert^2$ , with equality only if $w = p$

Why· The picture behind the proof

Here is the proof as a picture, with $W$ the line through $u$ . Think of the arrow $v$ as our $b$ and drag it anywhere. The green arrow is $p$ , and the red dashed segment is $b - p$ , at a right angle to the line. The orange point is some other point $w$ of the line; slide it with the slider. The orange dashed segment is $b - w$ , the hypotenuse of a right triangle whose legs are the red segment ( $b - p$ ) and the stretch of the line between $p$ and $w$ ( $p - w$ ). That is the Pythagoras step of the proof, and the hypotenuse is always longer than the red leg.

Quick checkWhat tells you that a point

p

of a subspace

W

is the closest point of

W

b

The theorem is about points of a subspace $W$ . Our subspace is the column space of $A$ , and the points of it are the vectors $Ax$ . Now we translate “the error is perpendicular to the column space” into a matrix equation.

Definition B

Normal equations

The normal equations for $Ax = b$ are

A^{\mathsf T}A\,x = A^{\mathsf T}b .

Here $A^{\mathsf T}A$ is $n \times n$ and $A^{\mathsf T}b$ is in $\mathbb R^n$ . (‘Normal’ is the old word for perpendicular.)

In one lineThe square system $A^{\mathsf T}A\,x = A^{\mathsf T}b$ . Its solutions are exactly the least squares solutions of $Ax = b$ .

Why it existsIt turns 'the residual is perpendicular to every column' into an ordinary linear system with as many equations as unknowns, which elimination can solve.

Theorem B Strang 4.3; Axler 6C

Least squares solutions solve the normal equations

For $\hat x \in \mathbb R^n$ and $r = b - A\hat x$ , the following are equivalent:

$\hat x$ is a least squares solution of $Ax = b$ ;
$A^{\mathsf T}r = O$ , that is, $a_j \cdot r = 0$ for every column $a_j$ ;
$A^{\mathsf T}A\,\hat x = A^{\mathsf T}b$ .

Moreover, for every $x$ , $\lVert b - Ax \rVert^2 = \lVert r \rVert^2 + \lVert A\hat x - Ax \rVert^2$ whenever (2) holds.

In one line $\hat x$ is a least squares solution of $Ax = b$ exactly when $A^{\mathsf T}A\,\hat x = A^{\mathsf T}b$ , which says the residual is perpendicular to every column of $A$ .

Why it is trueMoving $\hat x$ along a column changes the residual along that column. If the residual is not perpendicular to a column you can shrink it by moving that way, and if it is perpendicular to all of them the previous theorem says you are done.

ProofStrang 4.3

The idea: (2) and (3) are the same line of algebra. To show (2) implies (1), use the previous theorem with W the column space. To show (1) implies (2), suppose the residual is NOT perpendicular to some column and slide along that column to make it shorter, which contradicts (1).

(2) and (3) are equivalent.

Expand

A^{\mathsf T}r = A^{\mathsf T}(b - A\hat x) = A^{\mathsf T}b - A^{\mathsf T}A\hat x

. This is zero exactly when (3) holds.

Then

r \cdot (Ax) = \sum_j x_j\,(a_j \cdot r) = 0

for every

x

: the residual is perpendicular to the whole column space.

Suppose (2). Any vector in the column space is

Ax = \sum x_j a_j

, and dot products distribute over sums.

\lVert b - Ax \rVert^2 = \lVert r \rVert^2 + \lVert A\hat x - Ax \rVert^2 \ge \lVert r \rVert^2

for every

x

. This proves (2) implies (1), and the ‘moreover’ formula.

Apply the previous theorem with

W

the column space,

p = A\hat x

(which is in

W

), and

w = Ax

. Here

b - p = r

Fix

j

with

a_j \cdot r \ne 0

. For a number

t

, the residual at

\hat x + t e_j

r - t a_j

Now suppose (1) holds but (2) fails. Then some column has

a_j \cdot r \ne 0

. We will move

\hat x

along that column. The new vector

\hat x + t e_j

(with

e_j

the

j

-th standard basis vector) has residual

b - A\hat x - t a_j = r - t a_j

\lVert r - t a_j \rVert^2 = \lVert r \rVert^2 - 2t\,(a_j \cdot r) + t^2 \lVert a_j \rVert^2

Expand the square as in the previous proof. The cross term is linear in

t

and the last term is quadratic in

t

Take

t = \dfrac{a_j \cdot r}{\lVert a_j \rVert^2}

. Then

\lVert r - t a_j \rVert^2 = \lVert r \rVert^2 - \dfrac{(a_j \cdot r)^2}{\lVert a_j \rVert^2} < \lVert r \rVert^2

Choose

t

to make the right side as small as possible. Because

a_j \cdot r \ne 0

we have

a_j \ne O

, so we may divide by

\lVert a_j \rVert^2

.This is where the idea comes from: a downhill step along one direction. The best step is the bottom of a parabola in t.

So (1) forces

a_j \cdot r = 0

for every

j

, which is (2). All three are equivalent.

This contradicts (1):

\hat x + t e_j

would do strictly better than

\hat x

Proof workshopRebuild the hard direction: a least squares solution has a residual perpendicular to every column

Tap the lines in the order a proof would write them. Each line can only use what is already above it.

Why· The picture behind the proof

Go back to the best-fit line picture at the top and tick “try your own line”. The two sliders $c$ and $m$ move $\hat x$ along the two columns of $A$ : $c$ moves along the all-ones column and $m$ along the $t$ column. The sum of squares is smallest exactly when neither slider can improve it. That is the first half of the proof: if the residual were not perpendicular to some column, the slider for that column could still lower the total. At the best line, the residual (the vertical segments) is perpendicular to both columns, which is what the readout reports as “the best possible”. In the one-column picture above, the green arrow is the stopping place of the slider.

Quick checkWhich of these is NOT a way of saying that

\hat x

is a least squares solution of

Ax = b

Beyond the course: for understanding· Does a least squares solution always exist?

Yes. The orthogonal decomposition theorem of 2.9 writes $b = p + e$ with $p$ in the column space and $e$ perpendicular to it. Since $p$ is in the column space, $p = A\hat x$ for some $\hat x$ , and then $b - A\hat x = e$ is perpendicular to every column, so (2) holds. When the columns of $A$ are independent, the next theorem gives $\hat x$ by a formula and needs no other result.

Worked example Spagnolie Lecture 1

Spagnolie's three equations in two unknowns: find the best

x

for

3x_1 + x_2 = 1

2x_1 - x_2 = 0

4x_1 + 7x_2 = 2

Solution.

Write the system as $Ax = b$ . $A = \begin{bmatrix} 3 & 1 \\ 2 & -1 \\ 4 & 7 \end{bmatrix}$ and $b = (1, 0, 2)$ . We showed at the start of the chapter that no exact solution exists, so we look for the closest point of the column space instead.
Build the normal equations. $A^{\mathsf T}A = \begin{bmatrix} 29 & 29 \\ 29 & 51 \end{bmatrix}$ and $A^{\mathsf T}b = (11, 15)$ . Entry $(i, j)$ of $A^{\mathsf T}A$ is the dot product $a_i \cdot a_j$ of two columns, and entry $i$ of $A^{\mathsf T}b$ is $a_i \cdot b$ .
Solve. Subtracting the first equation $29x_1 + 29x_2 = 11$ from the second gives $22x_2 = 4$ , so $\hat x_2 = \tfrac2{11}$ and $\hat x_1 = \tfrac{11}{29} - \tfrac2{11} = \tfrac{63}{319}$ . Two equations in two unknowns. The columns are independent ( $a_2$ is not a multiple of $a_1$ ), so the solution is unique.
Find the closest point and the residual. $p = A\hat x = \tfrac1{319}(247, 68, 658)$ and $r = b - p = \tfrac1{319}(72, -68, -20)$ . The residual is the leftover, data minus the closest point.
Check that the residual is perpendicular to both columns. $a_1 \cdot r = \tfrac{216 - 136 - 80}{319} = 0$ and $a_2 \cdot r = \tfrac{72 + 68 - 140}{319} = 0$ . This is the defining property of the answer, and it is checked by a route independent of the solving.
Measure the miss. $\lVert r \rVert^2 = \tfrac{10208}{319^2} = \tfrac{32}{319} \approx 0.100$ . The first two equations alone give $(0.2, 0.4)$ , which misses the third by $\tfrac85$ and has squared miss $2.56$ . The least squares answer is about 25 times better by this measure.

The step-through below shows the same computation with every product written out.

Answer. $\hat x = (63/319,\ 2/11) \approx (0.197, 0.182)$ , with closest point $p = (247, 68, 658)/319$ and squared error $32/319$ .

Worked example

Find the point of the plane spanned by

(1, 1, 0)

and

(0, 1, 1)

that is closest to

b = (1, 2, 3)

, and the distance from

b

to the plane.

Solution.

Put the spanning vectors in the columns. $A = \begin{bmatrix} 1 & 0 \\ 1 & 1 \\ 0 & 1 \end{bmatrix}$ . The plane is the column space of $A$ , so the closest point is $A\hat x$ for the least squares solution $\hat x$ .
Form the normal equations. $A^{\mathsf T}A = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$ (the dot products $2, 1, 2$ of the columns) and $A^{\mathsf T}b = (1 + 2,\ 2 + 3) = (3, 5)$ . The residual must be perpendicular to both columns, which is $A^{\mathsf T}A\,x = A^{\mathsf T}b$ .
Solve. $2x_1 + x_2 = 3$ gives $x_2 = 3 - 2x_1$ . Then $x_1 + 2(3 - 2x_1) = 5$ , so $-3x_1 = -1$ and $x_1 = 1/3$ , $x_2 = 7/3$ . The columns are independent, so there is exactly one solution.
Compute the closest point. $p = \tfrac13(1, 1, 0) + \tfrac73(0, 1, 1) = (1/3,\ 8/3,\ 7/3)$ . $A\hat x$ is the combination of the columns with weights $\hat x$ .
Find the residual and check it. $r = b - p = (2/3, -2/3, 2/3)$ . Then $r \cdot (1, 1, 0) = 0$ and $r \cdot (0, 1, 1) = 0$ . A residual that is perpendicular to both columns confirms the arithmetic by an independent route.
Read off the distance. $\lVert r \rVert = \sqrt{3 \cdot \tfrac49} = 2/\sqrt3 \approx 1.155$ . The distance from $b$ to the plane is the length of the perpendicular error.

Answer. Closest point $(1/3, 8/3, 7/3)$ , distance $2/\sqrt3$ .

Worked example Spagnolie HW4 (Piazza @38)

Prove that if

\hat x

is a least squares solution of

Ax = b

, with

p = A\hat x

and

r = b - p

, then

\lVert b \rVert^2 = \lVert p \rVert^2 + \lVert r \rVert^2

. Then show that the squared error can be computed as

\lVert b \rVert^2 - \hat x \cdot (A^{\mathsf T}b)

Solution.

Recall what the theorem gives. The residual $r$ is perpendicular to every column of $A$ , so $r \cdot (Ax) = 0$ for every $x$ . This is condition (2) of the normal equations theorem, extended from the columns to all their combinations.
Apply it to $x = \hat x$ . $r \cdot p = r \cdot (A\hat x) = 0$ . The closest point $p$ is one of the vectors $Ax$ .
Split $b$ . $b = p + r$ . That is the definition of the residual.
Expand the squared length. $\lVert b \rVert^2 = (p + r)\cdot(p + r) = \lVert p \rVert^2 + 2\,p \cdot r + \lVert r \rVert^2$ . The dot product distributes over sums and is symmetric.
Use step 2. The middle term is $0$ , so $\lVert b \rVert^2 = \lVert p \rVert^2 + \lVert r \rVert^2$ . Pythagoras: $p$ and $r$ are the two perpendicular legs.
Rewrite $\lVert p \rVert^2$ . $\lVert p \rVert^2 = (A\hat x)\cdot(A\hat x) = \hat x^{\mathsf T}A^{\mathsf T}A\hat x = \hat x^{\mathsf T}A^{\mathsf T}b = \hat x \cdot (A^{\mathsf T}b)$ . $(A\hat x)^{\mathsf T} = \hat x^{\mathsf T}A^{\mathsf T}$ , and the normal equations replace $A^{\mathsf T}A\hat x$ by $A^{\mathsf T}b$ .
Substitute. $\lVert r \rVert^2 = \lVert b \rVert^2 - \hat x \cdot (A^{\mathsf T}b)$ . Solve step 5 for $\lVert r \rVert^2$ and use step 6. This needs only numbers we already computed when solving the normal equations.

Answer. Proved. A useful consequence: $\lVert r \rVert^2 = \lVert b \rVert^2 - \lVert p \rVert^2$ . Check with the first example: $\lVert b \rVert^2 = 5$ and $\lVert p \rVert^2 = 498597/101761 \approx 4.8997$ , and $5 - 4.8997 \approx 0.1003 = 32/319$ .

Your answerFor some data,

\lVert b \rVert^2 = 50

, and the least squares solution has

\hat x \cdot (A^{\mathsf T}b) = 42

. What is the sum of the squared residuals?

When is the answer a single vector? The key is a small lemma.

Lemma B Strang 4.2

A and A-transpose A have the same null space

For every real $m \times n$ matrix $A$ , the null space of $A^{\mathsf T}A$ equals the null space of $A$ . (We write $N(M)$ for the null space of a matrix $M$ , the set of $x$ with $Mx = O$ .)

In one line $A^{\mathsf T}Ax = O$ exactly when $Ax = O$ .

Why it is trueMultiplying by $A^{\mathsf T}$ could lose information if it sent a non-zero vector to zero, and a squared length rules that out.

ProofStrang 4.2

The idea: One inclusion is free. For the other, dot the equation with x itself, so that it becomes a squared length.

Ax = O

gives

A^{\mathsf T}Ax = A^{\mathsf T}O = O

. So

N(A) \subseteq N(A^{\mathsf T}A)

Ax = O

, multiply both sides by

A^{\mathsf T}

x^{\mathsf T}A^{\mathsf T}Ax = x^{\mathsf T}O = 0

Conversely, start with

A^{\mathsf T}Ax = O

and multiply on the left by the row vector

x^{\mathsf T}

\lVert Ax \rVert^2 = 0

(Ax)^{\mathsf T} = x^{\mathsf T}A^{\mathsf T}

, the rule for the transpose of a product. So the left side is

(Ax)^{\mathsf T}(Ax) = (Ax)\cdot(Ax)

.This is the trick of the whole chapter: a matrix times its own transpose produces a squared length.

Ax = O

. So

N(A^{\mathsf T}A) \subseteq N(A)

, and the two sets are equal.

The only vector of length zero is

O

Theorem B Strang 4.3

Independent columns give a unique least squares solution

(a) The columns of $A$ are linearly independent if and only if $A^{\mathsf T}A$ (an $n \times n$ matrix) is invertible.

(b) If they are, $Ax = b$ has exactly one least squares solution, $\hat x = (A^{\mathsf T}A)^{-1}A^{\mathsf T}b$ .

(c) In general, the fitted vector $A\hat x$ is the same for every least squares solution $\hat x$ . Different solutions differ by a vector in the null space of $A$ .

In one lineThe columns of $A$ are independent exactly when $A^{\mathsf T}A$ is invertible. Then the unique least squares solution is $\hat x = (A^{\mathsf T}A)^{-1}A^{\mathsf T}b$ .

Why it is trueIndependent columns mean each point of the column space comes from only one $x$ , so the closest point has only one name.

ProofStrang 4.3

The idea: Independence says the null space of A is trivial. The lemma moves that statement to A-transpose A. A square matrix with trivial null space is invertible, so the normal equations can be solved by multiplying by the inverse.

The columns are independent iff

N(A) = \{O\}

By definition (2.2),

a_1, \dots, a_n

are independent exactly when

x_1 a_1 + \dots + x_n a_n = O

only for

x = O

. That sum is

Ax

N(A) = \{O\}

iff

N(A^{\mathsf T}A) = \{O\}

The lemma proved above.

A square matrix

M

has

N(M) = \{O\}

iff

M

is invertible. This proves (a).

For a square matrix, trivial null space means one-to-one, which by rank–nullity (2.3) means onto, hence invertible (2.4).This is the problem in 2.3: between spaces of the same finite dimension, one-to-one and onto are the same thing.

\hat x = (A^{\mathsf T}A)^{-1}A^{\mathsf T}b

is the only least squares solution. (It exists because we wrote it down.)

Part (b). If

A^{\mathsf T}A

is invertible, the normal equations

A^{\mathsf T}Ax = A^{\mathsf T}b

have exactly one solution, found by multiplying both sides by the inverse. By the previous theorem the solutions of the normal equations are exactly the least squares solutions.

\lVert b - A\hat x' \rVert^2 = \lVert b - A\hat x \rVert^2 + \lVert A\hat x - A\hat x' \rVert^2

and both sides’ first terms are equal, so

A\hat x = A\hat x'

, that is

\hat x - \hat x' \in N(A)

Part (c). If

\hat x

and

\hat x'

are both least squares solutions, apply the 'moreover' formula of the previous theorem with

x = \hat x'

. Both have the same minimum value, so the extra term must be zero.

Quick checkThe columns of

A

are dependent and

b

is arbitrary. Which statement is true?

Remember part (c). The closest point is always unique, as a point of the plane. When the columns are dependent, several vectors $x$ give that same point, so the point has several NAMES as $Ax$ . The names are what can fail to be unique.

Worked example the picture at the top

Fit a line

y = c + mt

to the points

(-1, 0), (0, 1), (1, 1), (2, 3)

by least squares.

Solution.

Write one equation per point. $c + mt_i = y_i$ gives $c - m = 0$ , $c = 1$ , $c + m = 1$ , $c + 2m = 3$ . The line must pass through $(t_i, y_i)$ . No line passes through all four, so these four equations are inconsistent, and we look for the closest.
Put them in matrix form. $A = \begin{bmatrix} 1 & -1 \\ 1 & 0 \\ 1 & 1 \\ 1 & 2 \end{bmatrix}$ , $x = (c, m)$ , $b = (0, 1, 1, 3)$ . The first column multiplies $c$ (all ones) and the second multiplies $m$ (the $t$ values). The column space is the set of all heights some line produces.
Build the normal equations. $A^{\mathsf T}A = \begin{bmatrix} 4 & 2 \\ 2 & 6 \end{bmatrix}$ and $A^{\mathsf T}b = (5, 7)$ . The entries are $\sum 1 = 4$ , $\sum t_i = 2$ , $\sum t_i^2 = 6$ , $\sum y_i = 5$ and $\sum t_iy_i = 7$ .
Solve. $4c + 2m = 5$ and $2c + 6m = 7$ give $m = \tfrac9{10}$ and $c = \tfrac45$ . From the first, $c = (5 - 2m)/4$ . Substituting into the second: $5 + 10m = 14$ .
Find the residuals and check them. The fitted heights are $(-0.1, 0.8, 1.7, 2.6)$ and the residuals are $(0.1, 0.2, -0.7, 0.4)$ . They sum to $0$ and $\sum t_ir_i = -0.1 + 0 - 0.7 + 0.8 = 0$ . The residual must be perpendicular to both columns, the all-ones column and the $t$ column.
Measure the error. $\lVert r \rVert^2 = 0.01 + 0.04 + 0.49 + 0.16 = 0.70$ . The line $y = 0.5 + t$ has residuals $\pm0.5$ with sum of squares $1.00$ , which is worse.

The step-through below shows the same computation with every product written out.

Answer. $y = 0.8 + 0.9t$ , with sum of squared residuals $0.70$ . The same recipe fits any model that is linear in its unknowns: to fit a parabola $c + mt + qt^2$ , add a third column holding the values $t_i^2$ and change nothing else.

Quick checkWhat does the least squares line minimize?

Your answerFit a line to the points

(0, 0), (1, 1), (2, 3)

. What is the slope

m

Your answerThe constant model

y = c

is fitted to the data

2, 4, 9

. The model's single column is

(1, 1, 1)

. What is the least squares value of

c

Worked example

Fit a parabola y = c + mt + qt² to the points (0, 0), (1, 1), (2, 2), (3, 5) by least squares.

Solution.

Choose the columns. $A$ has columns $\mathbf 1$ , $t$ and $t^2$ , so the rows are $(1, 0, 0), (1, 1, 1), (1, 2, 4), (1, 3, 9)$ , and $b = (0, 1, 2, 5)$ . The model is linear in the unknowns $c, m, q$ even though it is quadratic in $t$ . Each unknown multiplies one column.
Collect the sums. $\sum 1 = 4$ , $\sum t = 6$ , $\sum t^2 = 14$ , $\sum t^3 = 36$ , $\sum t^4 = 98$ , and $\sum y = 8$ , $\sum ty = 0 + 1 + 4 + 15 = 20$ , $\sum t^2y = 0 + 1 + 8 + 45 = 54$ . Entry $(i, j)$ of $A^{\mathsf T}A$ is the dot product of columns $i$ and $j$ , which is a sum of a power of $t$ .
Write the normal equations. $\begin{bmatrix} 4 & 6 & 14 \\ 6 & 14 & 36 \\ 14 & 36 & 98 \end{bmatrix}\begin{bmatrix} c \\ m \\ q \end{bmatrix} = \begin{bmatrix} 8 \\ 20 \\ 54 \end{bmatrix}$ . $A^{\mathsf T}A\,x = A^{\mathsf T}b$ , with the sums from step 2.
Solve by elimination. $(c, m, q) = (0.1, 0.1, 0.5)$ . The columns are independent (the points have four distinct $t$ values), so there is exactly one solution.
Check in the equations. $4(0.1) + 6(0.1) + 14(0.5) = 8$ , $6(0.1) + 14(0.1) + 36(0.5) = 20$ and $14(0.1) + 36(0.1) + 98(0.5) = 54$ . Substituting back catches elimination slips.
Check the residual. The fitted heights are $0.1, 0.7, 2.3, 4.9$ , so the residual is $(-0.1, 0.3, -0.3, 0.1)$ . It is perpendicular to all three columns. With $\mathbf 1$ : $-0.1 + 0.3 - 0.3 + 0.1 = 0$ . With $t$ : $0.3 - 0.6 + 0.3 = 0$ . With $t^2$ : $0.3 - 1.2 + 0.9 = 0$ .

Answer. $y = 0.1 + 0.1t + 0.5t^2$ , with sum of squared residuals $0.01 + 0.09 + 0.09 + 0.01 = 0.20$ .

Worked example

Fit a line through the origin,

y = mt

, to the points

(1, 1), (2, 3), (4, 3)

Solution.

Write the system. There is one unknown $m$ , so $A$ is the single column $(1, 2, 4)$ and $b = (1, 3, 3)$ . The equations are $m \cdot t_i = y_i$ . A line through the origin has no constant term, so there is no column of ones.
Form the normal equation. $A^{\mathsf T}A = 1 + 4 + 16 = 21$ and $A^{\mathsf T}b = 1 + 6 + 12 = 19$ . With one column, the matrices are $1 \times 1$ numbers: $(t \cdot t)\,m = t \cdot y$ .
Solve. $m = 19/21 \approx 0.905$ . Dividing by $21$ is allowed because the column is not zero, which is what independence means for one column.
Compute the residuals. $r = b - m\,(1, 2, 4) = (2/21,\ 25/21,\ -13/21)$ . Data minus fitted values.
Check perpendicularity. $r \cdot (1, 2, 4) = (2 + 50 - 52)/21 = 0$ . An independent check that the arithmetic and the theory agree.
Size of the error. $\lVert r \rVert^2 = (4 + 625 + 169)/441 = 38/21$ . Check by the shortcut $\lVert b \rVert^2 - \hat x \cdot A^{\mathsf T}b = 19 - \tfrac{19}{21}\cdot 19 = \tfrac{399 - 361}{21} = \tfrac{38}{21}$ , which agrees.

Answer. $y = \tfrac{19}{21}t$ , with squared error $38/21$ . This is projection of $b$ onto the line through $(1, 2, 4)$ , the situation in the first picture.

Worked example

Find all least squares solutions of

Ax = b

when

A

is the

3 \times 2

matrix whose two columns are both

(1, 1, 1)

, and

b = (1, 2, 3)

Solution.

Form the normal equations. $A^{\mathsf T}A = \begin{bmatrix} 3 & 3 \\ 3 & 3 \end{bmatrix}$ and $A^{\mathsf T}b = (6, 6)$ . The theorem says least squares solutions are exactly the solutions of these equations, whether or not the columns are independent.
Notice the formula fails. $\det A^{\mathsf T}A = 9 - 9 = 0$ , so $(A^{\mathsf T}A)^{-1}$ does not exist. The two columns are equal, so they are dependent, and the invertibility theorem says $A^{\mathsf T}A$ is singular.
Solve the system anyway. Both equations say $3x_1 + 3x_2 = 6$ , that is $x_1 + x_2 = 2$ . So the solutions are $\hat x = (2 - s,\ s)$ for any number $s$ . One equation in two unknowns leaves one free variable, as in chapter 2.6.
Compare the fitted vectors. $A\hat x = (2 - s + s)(1, 1, 1) = (2, 2, 2)$ for every $s$ . Part (c) of the theorem: the closest point is unique, only its name as $Ax$ is not.
Check the residual. $r = (1, 2, 3) - (2, 2, 2) = (-1, 0, 1)$ , and $r \cdot (1, 1, 1) = 0$ . It is perpendicular to the (only) direction in the column space.

Answer. $\hat x = (2 - s, s)$ for every real $s$ ; the fitted vector is $(2, 2, 2)$ , the average of the data in every entry. Two copies of one feature cannot be told apart, so only their sum is determined. The solutions differ by multiples of $(1, -1)$ , which is in the null space of $A$ .

Worked example new (ridge regression)

Let A be any real m by n matrix, with columns possibly dependent, and let λ be a positive number. Prove that AᵀA + λI is invertible.

Solution.

Aim at the null space. Suppose $(A^{\mathsf T}A + \lambda I)x = O$ . We show $x = O$ . A square matrix with null space $\{O\}$ is invertible, by rank–nullity (2.3), as in the uniqueness proof above.
Multiply on the left by $x^{\mathsf T}$ . $x^{\mathsf T}A^{\mathsf T}Ax + \lambda x^{\mathsf T}x = 0$ . The right side stays $0$ . This is the trick of the chapter again: a matrix times its own transpose makes a squared length.
Recognise squared lengths. $\lVert Ax \rVert^2 + \lambda\lVert x \rVert^2 = 0$ . $x^{\mathsf T}A^{\mathsf T}Ax = (Ax)\cdot(Ax)$ and $x^{\mathsf T}x = x\cdot x$ .
Both terms are non-negative. Both $\lVert Ax \rVert^2 \ge 0$ and $\lambda\lVert x \rVert^2 \ge 0$ , and they add to $0$ , so each is $0$ . This is where $\lambda > 0$ is used. If $\lambda$ were $0$ , the second term would vanish for free and we could not conclude anything about $x$ .
Conclude. $\lambda\lVert x \rVert^2 = 0$ and $\lambda > 0$ give $\lVert x \rVert = 0$ , so $x = O$ . Only the zero vector has length zero.

Answer. $A^{\mathsf T}A + \lambda I$ is invertible for every $\lambda > 0$ , so ridge regression $(A^{\mathsf T}A + \lambda I)x = A^{\mathsf T}b$ always has exactly one solution, even for dependent columns. For the two equal columns above, $\begin{bmatrix}3&3\\3&3\end{bmatrix} + \lambda I$ has determinant $(3 + \lambda)^2 - 9 = \lambda(6 + \lambda) \ne 0$ .

Your answerA single column

a = (1, 2, 2)

and

b = (1, 1, 1)

are given. Solve the ridge equation

(a^{\mathsf T}a + \lambda)\,x = a^{\mathsf T}b

with

\lambda = 1

. What is

x

Worked example

Prove that the least squares line

y = c + mt

through

N

data points passes through the point

(\bar t, \bar y)

of the averages.

Solution.

Name the pieces. $A$ has two columns: the all-ones vector $\mathbf 1 = (1, \dots, 1)$ and the vector of the $t_i$ . The residual is $r_i = y_i - c - mt_i$ . The residual of the best fit is data minus fitted line, point by point.
Use perpendicularity to the first column. $\mathbf 1 \cdot r = 0$ , so $\sum_i r_i = 0$ . By the normal equations theorem, $r$ is perpendicular to every column of $A$ , and this column is $\mathbf 1$ .
Substitute the formula for $r_i$ . $0 = \sum_i (y_i - c - mt_i) = \sum_i y_i - Nc - m\sum_i t_i$ . Sums split term by term, and $\sum_i c = Nc$ because $c$ does not depend on $i$ .
Divide by $N$ . $0 = \bar y - c - m\bar t$ , that is $\bar y = c + m\bar t$ . $\sum_i y_i / N = \bar y$ and $\sum_i t_i / N = \bar t$ by the definition of an average.
Read the result. The point $(\bar t, \bar y)$ satisfies the line’s equation. That is what “the line passes through the point” means.

Answer. Proved. Check with the line fit above: $\bar t = 0.5$ , $\bar y = 1.25$ and $0.8 + 0.9 \cdot 0.5 = 1.25$ . The same proof shows that the residuals sum to zero whenever the model has a constant term.

The projection matrix

Part (b) gives the closest point $A\hat x = A(A^{\mathsf T}A)^{-1}A^{\mathsf T}b$ . Look at how it depends on $b$ : the vector $b$ is multiplied on the left by a fixed matrix that depends only on $A$ . That matrix is worth a name.

Definition B

Projection matrix

If $A$ has independent columns, the projection matrix onto its column space is

P = A\,(A^{\mathsf T}A)^{-1}A^{\mathsf T}.

It is $m \times m$ . In statistics it is called the hat matrix, because it puts a hat on $y$ : the fitted values are $\hat y = Py$ .

In one lineThe matrix $P = A(A^{\mathsf T}A)^{-1}A^{\mathsf T}$ that sends every $b$ to the closest point of the column space of $A$ .

Why it existsIt does the whole job at once: feed it any $b$ and it returns the fitted values, without solving a new system.

Theorem B Strang 4.2; Axler 6C

The projection matrix projects

Let $A$ have independent columns and $P = A(A^{\mathsf T}A)^{-1}A^{\mathsf T}$ . Then:

$Pb = A\hat x$ , the projection of $b$ onto the column space of $A$ ;
$Pv = v$ for every $v$ in the column space;
$P^2 = P$ ;
$P^{\mathsf T} = P$ ;
$(I - P)b = b - Pb$ is the residual, perpendicular to every column, so $A^{\mathsf T}(I - P) = O$ .

In one line $Pb$ is the closest point of the column space to $b$ , $P^2 = P$ and $P^{\mathsf T} = P$ .

Why it is trueProjecting a second time must change nothing, because the first projection already landed on the plane.

ProofStrang 4.2

The idea: Everything is algebra on the one expression A(A-transpose A)^{-1}A-transpose, plus the fact that the middle factor cancels against A-transpose A when it meets it.

(1)

Pb = A(A^{\mathsf T}A)^{-1}A^{\mathsf T}b = A\hat x

Substitute the formula

\hat x = (A^{\mathsf T}A)^{-1}A^{\mathsf T}b

from the previous theorem.

(2)

Pv = A(A^{\mathsf T}A)^{-1}A^{\mathsf T}Ax = A\,I\,x = Ax = v

A vector in the column space is

v = Ax

for some

x

. Then the middle factors cancel:

(A^{\mathsf T}A)^{-1}A^{\mathsf T}A = I

(3)

P^2 b = P(Pb) = Pb

for every

b

, so

P^2 = P

Pb

is in the column space by (1), so apply (2) to the vector

Pb

. This is also the geometric reason: a point already on the plane is its own shadow.

(4)

P^{\mathsf T} = (A^{\mathsf T})^{\mathsf T}\big((A^{\mathsf T}A)^{-1}\big)^{\mathsf T}A^{\mathsf T} = A\,(A^{\mathsf T}A)^{-1}A^{\mathsf T} = P

, since

A^{\mathsf T}A

is its own transpose.

Transpose reverses the order of a product:

(XYZ)^{\mathsf T} = Z^{\mathsf T}Y^{\mathsf T}X^{\mathsf T}

. The transpose of an inverse is the inverse of the transpose. And

(A^{\mathsf T}A)^{\mathsf T} = A^{\mathsf T}A

(5)

A^{\mathsf T}(I - P) = A^{\mathsf T} - A^{\mathsf T}A(A^{\mathsf T}A)^{-1}A^{\mathsf T} = A^{\mathsf T} - A^{\mathsf T} = O

Multiply out, using the cancellation

A^{\mathsf T}A(A^{\mathsf T}A)^{-1} = I

Quick check

P

is the

5 \times 5

projection matrix onto a 2-dimensional column space. Is

P

invertible?

A matrix with $P^2 = P$ is called idempotent, and one with $P^{\mathsf T} = P$ is symmetric. Properties 3 and 4 together (idempotent and symmetric) characterize orthogonal projections. A problem below proves the converse.

Why· The picture behind the proof

Look at the first picture of the chapter, the line through $u$ . The shadow $p$ is $Pv$ . Property 2 and property 3 say: if you take the shadow of a point that is already on the line (the green arrow), it stays where it is, so projecting twice is the same as once. Property 5 is the red dashed arrow: whatever the shadow leaves behind, $v - p$ , meets the line at a right angle, and $A^{\mathsf T}(I - P) = O$ says exactly that.

The projection matrix for the line-fit data. For the $A$ of the line-fit example,

P = \frac{1}{10}\begin{bmatrix} 7 & 4 & 1 & -2 \\ 4 & 3 & 2 & 1 \\ 1 & 2 & 3 & 4 \\ -2 & 1 & 4 & 7 \end{bmatrix}.

Applying it to $b = (0, 1, 1, 3)$ gives $Pb = (-0.1, 0.8, 1.7, 2.6)$ , the fitted heights, in one step. Its diagonal is $0.7, 0.3, 0.3, 0.7$ , which sums to $2$ , the number of columns. (That sum of the diagonal entries is the trace. It is no accident: the trace of $P$ always equals the number of independent columns. The proof uses $\operatorname{tr}(XY) = \operatorname{tr}(YX)$ , which holds because both equal $\sum_{i,j} X_{ij}Y_{ji}$ . Then $\operatorname{tr}P = \operatorname{tr}\big((A^{\mathsf T}A)^{-1}A^{\mathsf T}A\big) = \operatorname{tr}I_n = n$ .) The big diagonal entries at the two end points say those points have the most pull on the line, their “leverage”. Move one of them and the line follows.

Worked example

Let

P

be a projection matrix (so

P^2 = P

and

P^{\mathsf T} = P

). Prove that

I - P

is also symmetric and idempotent, that

(I - P)P = O

, and that every diagonal entry of

P

satisfies

0 \le P_{ii} \le 1

Solution.

Symmetry of $I - P$ . $(I - P)^{\mathsf T} = I^{\mathsf T} - P^{\mathsf T} = I - P$ . The transpose of a difference is the difference of transposes, and $I$ and $P$ are both symmetric.
Idempotence of $I - P$ . $(I - P)^2 = I - 2P + P^2 = I - 2P + P = I - P$ . Multiply out (matrices commute with $I$ and with themselves), then use $P^2 = P$ .
Killing what $P$ keeps. $(I - P)P = P - P^2 = P - P = O$ . $P$ lands in the plane and $I - P$ removes exactly the plane’s part, so applying $I - P$ after $P$ leaves nothing.
A formula for the diagonal. $P = P^2 = P^{\mathsf T}P$ , so $P_{ii}$ is the dot product of column $i$ of $P$ with itself: $P_{ii} = \sum_k P_{ki}^2$ . Entry $(i, i)$ of $P^{\mathsf T}P$ multiplies row $i$ of $P^{\mathsf T}$ , which is column $i$ of $P$ , by column $i$ of $P$ .
Non-negativity. A sum of squares is at least $0$ , so $P_{ii} \ge 0$ . Squares of real numbers are never negative.
The upper bound. The term $k = i$ of the sum is $P_{ii}^2$ and the other terms are at least $0$ , so $P_{ii} \ge P_{ii}^2$ . If $P_{ii} > 0$ , dividing by it gives $1 \ge P_{ii}$ . If $P_{ii} = 0$ the bound holds already. Dropping non-negative terms can only make a sum smaller.

Answer. Proved. Check with the line-fit $P$ : $P_{11} = 0.49 + 0.16 + 0.01 + 0.04 = 0.70$ , the sum of the squares of column 1, and $0 \le 0.7 \le 1$ .

Your answer

A

is a

10 \times 3

matrix with independent columns and

P

is its projection matrix. What is the trace of

P

(the sum of its diagonal entries)?

Worked example Axler 6C

A is a 5 by 2 matrix with independent columns and P is its projection matrix. Give the sizes of AᵀA, Aᵀb and P, the rank of P, and say whether P is invertible. Then prove that for any symmetric idempotent matrix P, b − Pb is perpendicular to the range of P.

Solution.

Sizes. $A^{\mathsf T}A$ is $2 \times 2$ , $A^{\mathsf T}b$ is a vector in $\mathbb R^2$ and $P = A(A^{\mathsf T}A)^{-1}A^{\mathsf T}$ is $5 \times 5$ . The product of sizes $(5 \times 2)(2 \times 2)(2 \times 5)$ is $5 \times 5$ . The normal equations always have as many unknowns as columns.
Rank. Every output $Pb$ lies in the column space of $A$ , and every vector of the column space is fixed by $P$ , so it is an output. The range of $P$ is the column space, of dimension $2$ , so $\operatorname{rank}P = 2$ . Property (2) of the projection matrix theorem says $Pv = v$ for $v$ in the column space.
Invertibility. $P$ is $5 \times 5$ of rank $2$ , so it is not invertible. An invertible $5 \times 5$ matrix has rank $5$ (2.4). The part of $b$ perpendicular to the plane is lost and cannot be recovered.
The proof: pick a vector in the range. Any vector in the range of $P$ is $Pz$ for some $z$ . That is what “range” means.
Compute the two inner products. $(Pb)\cdot(Pz) = b^{\mathsf T}P^{\mathsf T}Pz = b^{\mathsf T}P^2z = b^{\mathsf T}Pz$ , and $b\cdot(Pz) = b^{\mathsf T}Pz$ . First $P^{\mathsf T} = P$ , then $P^2 = P$ .
Subtract. $(b - Pb)\cdot(Pz) = b^{\mathsf T}Pz - b^{\mathsf T}Pz = 0$ . Both inner products from step 5 are the same number.
Conclude. By the first theorem of the chapter, with $W$ the range of $P$ and $p = Pb$ (which lies in $W$ ), $Pb$ is the closest point of $W$ to $b$ . Perpendicular error means closest point.

Answer. Sizes $2 \times 2$ , $\mathbb R^2$ and $5 \times 5$ ; rank $2$ ; not invertible. And any symmetric idempotent matrix is the orthogonal projection onto its range.

QR: how computers do it (B)

The formula $\hat x = (A^{\mathsf T}A)^{-1}A^{\mathsf T}b$ is perfect on paper and risky in a computer. Forming $A^{\mathsf T}A$ squares the “difficulty” of the problem. The condition number of a matrix measures how much a tiny error in the input can be magnified in the output. For the standard condition number (the biggest stretching the matrix does divided by the smallest), $\operatorname{cond}(A^{\mathsf T}A) = \operatorname{cond}(A)^2$ when the columns of $A$ are independent. We state this fact without proof. A concrete failure shows the damage.

Take $\varepsilon = 10^{-8}$ and $A = \begin{bmatrix} 1 & 1 \\ \varepsilon & 0 \\ 0 & \varepsilon \end{bmatrix}$ . Its columns are independent. In exact arithmetic,

A^{\mathsf T}A = \begin{bmatrix} 1 + \varepsilon^2 & 1 \\ 1 & 1 + \varepsilon^2 \end{bmatrix}.

But $\varepsilon^2 = 10^{-16}$ is below the rounding resolution of double-precision numbers, so the computer stores $1 + 10^{-16}$ as exactly $1$ . The stored matrix is $\begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}$ , which is singular. The normal equations cannot be solved, although the original problem is perfectly fine. (Numerical analysts call this the Läuchli example. Running numpy’s solve on A.T @ A reports a singular matrix here, while numpy.linalg.lstsq(A, b) and the QR method below give the right answer.)

The fix is to use the Gram–Schmidt process from 2.9 on the columns of $A$ and never form $A^{\mathsf T}A$ .

Theorem B Strang 4.4

QR from Gram–Schmidt

Let $A$ be $m \times n$ with independent columns $a_1, \dots, a_n$ . Gram–Schmidt produces orthonormal vectors $q_1, \dots, q_n$ with $\operatorname{span}(q_1, \dots, q_j) = \operatorname{span}(a_1, \dots, a_j)$ for each $j$ . Let $Q$ be the $m \times n$ matrix with columns $q_j$ . Then there is an upper triangular $n \times n$ matrix $R$ with nonzero diagonal entries such that $A = QR$ . Consequently $Q^{\mathsf T}Q = I$ , $A^{\mathsf T}A = R^{\mathsf T}R$ , the least squares solution satisfies $R\hat x = Q^{\mathsf T}b$ , and the projection matrix is $P = QQ^{\mathsf T}$ .

In one lineA matrix $A$ with independent columns factors as $A = QR$ : $Q$ has orthonormal columns and $R$ is square, upper triangular and invertible. Then $\hat x$ solves $R\hat x = Q^{\mathsf T}b$ and $P = QQ^{\mathsf T}$ .

Why it is trueGram–Schmidt rewrites each column as a combination of earlier orthonormal vectors. Writing that down column by column is a matrix product, and orthonormal columns make the normal equations collapse.

ProofStrang 4.4

The idea: Column j of A lies in the span of the first j of the q's, so it is a combination of only those. The coefficients are dot products. Stack them into R.

a_j = R_{1j}q_1 + R_{2j}q_2 + \dots + R_{jj}q_j

for some numbers

R_{ij}

Gram–Schmidt (2.9) gives

\operatorname{span}(q_1, \dots, q_j) = \operatorname{span}(a_1, \dots, a_j)

, so

a_j

is in the span of

q_1, \dots, q_j

R_{ij} = q_i \cdot a_j

for

i \le j

. Put

R_{ij} = 0

for

i > j

Dot both sides with

q_i

for

i \le j

. The

q

's are orthonormal, so

q_i \cdot q_k

1

i = k

and

0

otherwise, and every term but one vanishes.This is the Lecture 6 formula: coefficients along an orthonormal set are dot products.

A = QR

, and

R

is upper triangular.

Column

j

of the product

QR

\sum_i R_{ij}q_i

, the combination written in step 1.

Every

R_{jj} \ne 0

. A triangular matrix has determinant equal to the product of its diagonal (2.7), which is non-zero, so

R

is invertible.

R_{jj} = 0

, then

a_j

would be in

\operatorname{span}(q_1, \dots, q_{j-1}) = \operatorname{span}(a_1, \dots, a_{j-1})

, contradicting independence.

Q^{\mathsf T}Q = I

Entry

(i, k)

Q^{\mathsf T}Q

q_i \cdot q_k

The normal equations read

R^{\mathsf T}R\,\hat x = R^{\mathsf T}Q^{\mathsf T}b

A^{\mathsf T} = R^{\mathsf T}Q^{\mathsf T}

, so

A^{\mathsf T}A = R^{\mathsf T}Q^{\mathsf T}QR = R^{\mathsf T}R

and

A^{\mathsf T}b = R^{\mathsf T}Q^{\mathsf T}b

. The normal equations become

R^{\mathsf T}R\hat x = R^{\mathsf T}Q^{\mathsf T}b

R\hat x = Q^{\mathsf T}b

R^{\mathsf T}

is invertible because

R

is. Multiply both sides by

(R^{\mathsf T})^{-1}

. Since

R

is triangular, solving this is back-substitution: the last equation has one unknown, and so on upward.

P = QR\,(R^{\mathsf T}R)^{-1}R^{\mathsf T}Q^{\mathsf T} = QR\,R^{-1}(R^{\mathsf T})^{-1}R^{\mathsf T}Q^{\mathsf T} = QQ^{\mathsf T}

Substitute

A = QR

and

A^{\mathsf T}A = R^{\mathsf T}R

into the formula for

P

. Then

R R^{-1} = I

and

(R^{\mathsf T})^{-1}R^{\mathsf T} = I

The last line is Lecture 6’s formula in matrix form: $Pb = QQ^{\mathsf T}b = \sum_i (q_i \cdot b)\,q_i$ . The projection onto a subspace with an orthonormal basis is the sum of the shadows on each basis vector. In the picture: each $(q_i \cdot b)\,q_i$ is a shadow like the green arrow in the projection visualizer, and because the $q_i$ are perpendicular to each other, the shadows do not interfere and add up. Real libraries compute the factorization with Householder reflections, which are more stable than Gram–Schmidt, and numpy.linalg.lstsq uses the SVD. The idea is the same: work with orthonormal columns so the problem never gets worse.

Python labFit a line to noisy data with the normal equations

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
x = np.linspace(0, 10, 25)
y = 2.0 + 0.7 * x + rng.normal(0, 1.0, x.size)    # a true line plus noise

A = np.column_stack([np.ones_like(x), x])          # columns: 1 and x
c = np.linalg.solve(A.T @ A, A.T @ y)              # the normal equations
print("intercept, slope:", np.round(c, 4))
print("numpy's lstsq agrees:", np.round(np.linalg.lstsq(A, y, rcond=None)[0], 4))

r = y - A @ c                                      # the residual (error) vector
print("A^T r =", np.round(A.T @ r, 10))

plt.figure(figsize=(6, 3.5))
plt.scatter(x, y, s=14, label="data")
plt.plot(x, A @ c, color="C3", label="least squares line")
plt.legend(); plt.title("Line of best fit")

Try this. The last line printed is $A^{\mathsf T}r$ , the residual dotted with each column: it is zero, which is the perpendicularity the whole chapter rests on. Now fit a parabola: add the column x**2 to A (and change the plot to a smooth curve), or change the noise from 1.0 to 5.0 and watch the fitted slope drift from the true $0.7$ .

Examples and non-examples

The worked examples sit beside the ideas they practise, above. Here are three more checks.

Check with a formula from statistics. The slope should equal $\dfrac{\sum (t_i - \bar t)(y_i - \bar y)}{\sum (t_i - \bar t)^2}$ where $\bar t = 0.5$ and $\bar y = 1.25$ are the averages. The top is $\sum t_iy_i - 4\bar t \bar y = 7 - 2.5 = 4.5$ and the bottom is $\sum t_i^2 - 4\bar t^2 = 6 - 1 = 5$ , so the slope is $0.9$ . This is “covariance over variance”, and it agrees with the normal equations.

The same line by QR, with no $A^{\mathsf T}A$ .

An exact fit. If $b$ happens to be in the column space, the closest point is $b$ itself and the residual is $O$ . For example, the points $(-1, 1), (0, 0), (1, 1)$ lie exactly on $y = t^2$ . Fitting $c + mt + qt^2$ gives $(c, m, q) = (0, 0, 1)$ and zero residual. Real data with noise almost never fits a small model exactly. A model with as many columns as data points (with distinct $t$ values) always fits exactly, which is the danger of fitting too much.

Where it is used

Used inside maths· Inside maths: approximating a function, and Fourier

The argument of this chapter never used the fact that the vectors are lists of numbers. It only needs an inner product, so it works for functions too. Spagnolie’s HW4 problem asks for the best linear approximation $a + bx$ to $e^x$ on $[-1, 1]$ with $\langle f, g \rangle = \int_{-1}^{1} fg\,dx$ . The columns are now the functions $1$ and $x$ , and the normal equations have entries that are integrals: $\langle 1, 1 \rangle = 2$ , $\langle 1, x \rangle = 0$ , $\langle x, x \rangle = 2/3$ , $\langle e^x, 1 \rangle = e - 1/e$ and $\langle e^x, x \rangle = 2/e$ . Because $\langle 1, x \rangle = 0$ , the two equations separate: $2a = e - 1/e$ and $\tfrac23 b = 2/e$ . So $a = (e - 1/e)/2 \approx 1.175$ and $b = 3/e \approx 1.104$ . The course reached the same answer by projection and by calculus, and the calculus equations are exactly these normal equations.

When the columns are orthonormal (the $Q$ of the QR section), $A^{\mathsf T}A = I$ and the normal equations solve themselves: $\hat x = Q^{\mathsf T}b$ . For $\langle f, g \rangle = \int fg$ and the cosines and sines of Lecture 6, this is the Fourier series: each coefficient is one inner product, and the partial sums are the closest trigonometric polynomials to $f$ .

Used in the real world· Linear regression in ML and statistics (CS 532)

Predicting a house price from its area and number of bedrooms with a linear model is this chapter. Each of the $m$ houses is a row of $A$ : a 1 (for the intercept), the area, the bedrooms. The prices are $b$ . The weight vector $\hat x$ solves $A^{\mathsf T}A\,\hat x = A^{\mathsf T}b$ , and the predictions are $\hat y = A\hat x = Pb$ . The matrix $A^{\mathsf T}A$ has only 3 rows and columns, even for a million houses, which is the reason the normal equations were practical long before computers were large. “Polynomial regression” adds columns for $t^2, t^3, \dots$ . The model is nonlinear in $t$ and still linear in the unknowns, which is all this chapter needs, and the exact-fit example with $y = t^2$ was one. Libraries such as scikit-learn’s LinearRegression solve it with a least squares routine (not the normal equations) for the stability reason shown above.

When two features are nearly copies of each other (area in square feet and area in square metres), the columns are almost dependent, $A^{\mathsf T}A$ is nearly singular and the weights become huge and unstable. The standard cure is ridge regression: solve $(A^{\mathsf T}A + \lambda I)x = A^{\mathsf T}b$ for a small $\lambda > 0$ . A problem below proves the matrix is always invertible. Training a neural network’s last layer with a squared-error loss, and most “closed form” fitting in the first weeks of an ML course, use exactly these equations.

Used in the real world· GPS: a position fix from more satellites than unknowns

A GPS receiver measures, for each visible satellite, the travel time of a radio signal. Multiplied by the speed of light, that is a distance (a “pseudorange”). There are four unknowns: the receiver’s position $(x, y, z)$ and its clock error. The clock error matters because light travels about 300 metres in one microsecond, and a cheap quartz clock in a phone is off by much more than a microsecond. Four unknowns need four satellites. A receiver in open sky often sees eight or more, so it has more equations than unknowns, and the measurements disagree slightly because of noise. Least squares settles it.

The distance equations are not linear, so the receiver linearizes them around a guess. A two-dimensional toy version with the same mathematics: four beacons at $(0,0), (10,0), (0,10), (10,10)$ ; the receiver is truly at $(3, 4)$ and measures distances with small errors, getting $5.10, 7.86, 6.86, 9.17$ . Guess $x_0 = (4, 4)$ . The distances from the guess to the beacons are $5.657, 7.211, 7.211, 8.485$ . Near the guess, the distance to beacon $i$ changes with the position correction $\Delta$ as $u_i \cdot \Delta$ , where $u_i$ is the unit vector pointing from beacon $i$ to the guess. That gives four equations $U\Delta = \text{measured} - \text{predicted}$ with rows $u_i$ :

\begin{bmatrix} 0.707 & 0.707 \\ -0.832 & 0.555 \\ 0.555 & -0.832 \\ -0.707 & -0.707 \end{bmatrix}\Delta = \begin{bmatrix} -0.557 \\ 0.649 \\ -0.351 \\ 0.685 \end{bmatrix}.

The normal equations $U^{\mathsf T}U\Delta = U^{\mathsf T}(\text{right side})$ give $\Delta = (-0.803, -0.082)$ , so the new guess is $(3.197, 3.918)$ . Repeating the step (with the unit vectors recomputed at the new guess) settles at $(3.186, 3.928)$ , which is about 0.2 from the true $(3, 4)$ , even though each measured distance was off by up to 0.2. (This is Gauss–Newton iteration: a least squares problem solved at every step.) A real receiver does this with a fourth column of ones for the clock error and four or more satellites.

Used in the real world· Calibrating a sensor

A temperature sensor outputs a number of ADC counts. To turn counts into degrees you measure it at known temperatures and fit a line. Suppose the readings at $0, 25, 50, 100\,^\circ\text{C}$ are $102, 205, 310, 515$ counts. We fit $T = c + m \cdot \text{counts}$ , with the all-ones column and the counts column. The normal equations give $T \approx -24.73 + 0.2420 \cdot \text{counts}$ . The residuals are $0.05, 0.12, -0.28, 0.11$ degrees, so no calibration point is off by more than a third of a degree, and any later reading is converted with the line. Using four points for two unknowns means a single bad reading shows up as a large residual instead of silently ruining the calibration.

Used in the real world· Hooke's law from noisy spring data

Hooke’s law says the force of a spring is $F = ks$ for an extension $s$ , with no constant term. (We call the extension $s$ , because $x$ is already our name for the unknowns.) Hang weights on a spring and measure the extensions $s = 0.02, 0.04, 0.06, 0.08, 0.10$ metres against forces $F = 0.21, 0.38, 0.62, 0.79, 1.02$ newtons. There is one unknown $k$ , so $A$ is the single column $s$ and the system $ks_i = F_i$ has five equations. The normal equation is $(s \cdot s)\,k = s \cdot F$ , with $s \cdot s = 0.0220$ and $s \cdot F = 0.2218$ , so

\hat k = \frac{0.2218}{0.0220} \approx 10.08\ \text{N/m}.

This is projection onto a line, the first picture in this chapter. The residuals are $0.008, -0.023, 0.015, -0.017, 0.012$ newtons. They are perpendicular to $s$ (their dot product with $s$ is zero up to rounding), and the spring constant is close to $10$ . If you also allowed an offset (fitting $F = F_0 + ks$ ) you would get $k \approx 10.15$ and $F_0 \approx -0.005$ . The offset is tiny, which supports Hooke’s law with no constant term.

Traps

Watch out· The line minimizes VERTICAL misses, not distance to the line

Least squares measures the error in $y$ only, along the vertical segments of the picture. The perpendicular distance from a point to the line is a different quantity, and minimizing it is a different problem (total least squares). The perpendicular in our theorem is in $\mathbb R^m$ , between the whole data vector $b$ and the column space, not in the data plane.

Watch out· Dependent columns: A-transpose A is singular

The formula $(A^{\mathsf T}A)^{-1}A^{\mathsf T}b$ needs independent columns. Duplicate features, a column that is the sum of others, or fewer data points than columns all break it. The fitted values are still unique (part (c)). The weights are not. In statistics this is called collinearity, and the “dummy variable trap” is one case.

Watch out· Do not invert A-transpose A on a computer

The formula is for understanding. The Läuchli example above has independent columns, and the computer still reports a singular matrix after forming $A^{\mathsf T}A$ . Use QR or a library least squares routine. Also, never build the $m \times m$ matrix $P$ for large $m$ : a million data points would need a million by a million matrix. Compute $Pb$ as $Q(Q^{\mathsf T}b)$ .

Watch out· Squaring lets one wild point steer the line

The residuals are squared, so a point ten units off counts a hundred times as much as a point one unit off. Add a single outlier $(3, 10)$ to our four points and the best line changes from $y = 0.8 + 0.9t$ to $y = 0.8 + 2.2t$ , so the slope more than doubles. When data may contain wild values, robust methods (minimizing absolute values, for example) are used instead.

Watch out· P squared equals P does not mean P is the identity

A projection matrix is NOT invertible unless the column space is everything. $P$ throws away the part of $b$ perpendicular to the plane, and that information cannot be recovered. $P^2 = P$ only says that doing it twice is the same as once.

Watch out· Computing the size of the error (the HW4 mistake)

In MATH 375 HW4, several problems asked for $\lVert f - g \rVert^2$ where $g$ is the best approximation, and the long way is to square out and integrate. The short way is Pythagoras: since $r \perp p$ , $\lVert b \rVert^2 = \lVert p \rVert^2 + \lVert r \rVert^2$ , so $\lVert r \rVert^2 = \lVert b \rVert^2 - \lVert p \rVert^2$ . Check with the line fit: $\lVert b \rVert^2 = 0 + 1 + 1 + 9 = 11$ and $\lVert p \rVert^2 = 0.01 + 0.64 + 2.89 + 6.76 = 10.30$ , and $11 - 10.30 = 0.70$ , as before.

Problems

checkStrang 4.3, adapted

Find the least squares solution of $Ax = b$ for $A = \begin{bmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2 \end{bmatrix}$ , $b = (6, 0, 0)$ . Also find the closest point and the residual, and check that the residual is perpendicular to the columns.

Solution

$A^{\mathsf T}A = \begin{bmatrix} 3 & 3 \\ 3 & 5 \end{bmatrix}$ and $A^{\mathsf T}b = (6, 0)$ . The equations $3c + 3m = 6$ and $3c + 5m = 0$ subtract to $2m = -6$ , so $m = -3$ and $c = 5$ . The closest point is $A\hat x = (5, 2, -1)$ and the residual is $r = (1, -2, 1)$ . Check: $(1, 1, 1) \cdot r = 0$ and $(0, 1, 2) \cdot r = -2 + 2 = 0$ .

check

$A$ is a $5 \times 2$ matrix with independent columns. State the sizes of $A^{\mathsf T}A$ , of $A^{\mathsf T}b$ and of $P$ . What is the rank of $P$ ?

Solution

$A^{\mathsf T}A$ is $2 \times 2$ , $A^{\mathsf T}b$ is in $\mathbb R^2$ , and $P$ is $5 \times 5$ . Every output $Pb$ is in the column space, and every vector of the column space is an output (it is fixed by $P$ ), so the range of $P$ is the column space of $A$ , which has dimension $2$ . So the rank of $P$ is $2$ .

check

Why can the formula $(A^{\mathsf T}A)^{-1}A^{\mathsf T}b$ not be used for $A = \begin{bmatrix} 1 & 2 \\ 2 & 4 \\ 3 & 6 \end{bmatrix}$ ?

Solution

The second column is twice the first, so the columns are dependent. Then $A^{\mathsf T}A = \begin{bmatrix} 14 & 28 \\ 28 & 56 \end{bmatrix}$ has determinant $14 \cdot 56 - 28 \cdot 28 = 0$ and is not invertible, as the theorem says.

use

Fit a line $y = c + mt$ to $(1, 1), (2, 3), (3, 4)$ by least squares. What is the sum of the squared residuals?

Solution

$A = \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix}$ , $b = (1, 3, 4)$ . Then $A^{\mathsf T}A = \begin{bmatrix} 3 & 6 \\ 6 & 14 \end{bmatrix}$ and $A^{\mathsf T}b = (8, 19)$ . The determinant is $42 - 36 = 6$ . Solve $3c + 6m = 8$ and $6c + 14m = 19$ : doubling the first gives $6c + 12m = 16$ , so $2m = 3$ and $m = 3/2$ ; then $c = (8 - 9)/3 = -1/3$ . The line is $y = -\tfrac13 + \tfrac32 t$ . The fitted values are $7/6, 8/3, 25/6$ , the residuals are $-1/6, 1/3, -1/6$ (they sum to $0$ and $\sum t_ir_i = -1/6 + 2/3 - 1/2 = 0$ ), and the sum of squares is $\tfrac1{36} + \tfrac19 + \tfrac1{36} = \tfrac16$ .

use

Fit a line through the origin, $y = mt$ , to $(1, 2), (2, 3), (3, 7)$ .

Solution

$A$ is the single column $(1, 2, 3)$ and $b = (2, 3, 7)$ . The normal equation is $14m = 2 + 6 + 21 = 29$ , so $m = 29/14$ . The fitted values are $29/14, 29/7, 87/14$ , the residuals are $-1/14, -8/7, 11/14$ and their dot product with $(1, 2, 3)$ is $-1/14 - 16/7 + 33/14 = 0$ . The sum of squares is $(1 + 256 + 121)/196 = 378/196 = 27/14$ .

use

Find the point of the plane spanned by $u = (1, 0, 1)$ and $v = (0, 1, 1)$ that is closest to $b = (1, 1, 1)$ , and the distance.

Solution

Take $A$ with columns $u, v$ . Then $A^{\mathsf T}A = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$ and $A^{\mathsf T}b = (2, 2)$ . The equations $2x_1 + x_2 = 2$ and $x_1 + 2x_2 = 2$ give $x_1 = x_2 = 2/3$ . The closest point is $\tfrac23 u + \tfrac23 v = (2/3, 2/3, 4/3)$ , the residual is $(1/3, 1/3, -1/3)$ (perpendicular to $u$ and $v$ : $1/3 - 1/3 = 0$ for each), and the distance is $\sqrt{1/9 + 1/9 + 1/9} = 1/\sqrt3$ .

use

Fit a parabola $y = c + mt + qt^2$ to $(-1, 1), (0, 0), (1, 2), (2, 5)$ .

Solution

The columns of $A$ are $1$ , $t$ , $t^2$ at $t = -1, 0, 1, 2$ . The sums needed are $\sum 1 = 4$ , $\sum t = 2$ , $\sum t^2 = 6$ , $\sum t^3 = 8$ , $\sum t^4 = 18$ , $\sum y = 8$ , $\sum ty = -1 + 0 + 2 + 10 = 11$ , $\sum t^2y = 1 + 0 + 2 + 20 = 23$ . The normal equations are

\begin{bmatrix} 4 & 2 & 6 \\ 2 & 6 & 8 \\ 6 & 8 & 18 \end{bmatrix}\begin{bmatrix} c \\ m \\ q \end{bmatrix} = \begin{bmatrix} 8 \\ 11 \\ 23 \end{bmatrix}.

Eliminating gives $(c, m, q) = (0.3, 0.4, 1)$ . Check: $4(0.3) + 2(0.4) + 6 = 8$ , $2(0.3) + 6(0.4) + 8 = 11$ , $6(0.3) + 8(0.4) + 18 = 23$ . The parabola is $y = 0.3 + 0.4t + t^2$ , with fitted values $0.9, 0.3, 1.7, 5.1$ and residuals $0.1, -0.3, 0.3, -0.1$ .

proveHW4 shortcut

Let $\hat x$ be a least squares solution of $Ax = b$ . Prove that $\lVert A\hat x \rVert^2 = \hat x \cdot (A^{\mathsf T}b)$ , and hence that the sum of squared residuals is $\lVert b \rVert^2 - \hat x \cdot (A^{\mathsf T}b)$ .

Solution

Write $p = A\hat x$ . Then $\lVert p \rVert^2 = (A\hat x)\cdot(A\hat x) = \hat x^{\mathsf T}A^{\mathsf T}A\hat x$ . By the normal equations $A^{\mathsf T}A\hat x = A^{\mathsf T}b$ , so $\lVert p \rVert^2 = \hat x^{\mathsf T}A^{\mathsf T}b = \hat x \cdot (A^{\mathsf T}b)$ . The worked example on Pythagoras gave $\lVert r \rVert^2 = \lVert b \rVert^2 - \lVert p \rVert^2$ , so $\lVert r \rVert^2 = \lVert b \rVert^2 - \hat x \cdot (A^{\mathsf T}b)$ . Check with the line fit: $\lVert b \rVert^2 = 11$ and $\hat x \cdot A^{\mathsf T}b = 0.8 \cdot 5 + 0.9 \cdot 7 = 10.3$ , and $11 - 10.3 = 0.7$ .

proveAxler 6C

Let $P$ be an $m \times m$ matrix with $P^2 = P$ and $P^{\mathsf T} = P$ . Prove that for every $b$ , the vector $b - Pb$ is perpendicular to every vector in the range of $P$ . (So $Pb$ is the projection of $b$ onto the range of $P$ .)

Solution

A vector in the range of $P$ has the form $Pz$ . Then $(Pb)\cdot(Pz) = b^{\mathsf T}P^{\mathsf T}Pz = b^{\mathsf T}P^2z = b^{\mathsf T}Pz$ , using $P^{\mathsf T} = P$ and then $P^2 = P$ . And $b \cdot (Pz) = b^{\mathsf T}Pz$ . So $(b - Pb)\cdot(Pz) = b^{\mathsf T}Pz - b^{\mathsf T}Pz = 0$ . By the first theorem of the chapter (with $W$ the range of $P$ and $p = Pb$ , which is in $W$ ), $Pb$ is the closest point of $W$ to $b$ .

prove

You fit a plane $y = c + m_1 t + m_2 s$ to $N$ data points $(t_i, s_i, y_i)$ by least squares. Prove that the average of the fitted values equals the average of the data, and that the plane passes through the point of averages $(\bar t, \bar s, \bar y)$ .

Solution

The first column of $A$ is the all-ones vector $\mathbf 1$ , and the residual $r = y - A\hat x$ is perpendicular to it, so $\sum_i r_i = 0$ . Since $r_i = y_i - \hat y_i$ with $\hat y_i = c + m_1t_i + m_2s_i$ , this says $\sum_i \hat y_i = \sum_i y_i$ , and dividing by $N$ gives equal averages. Also $\frac1N\sum_i \hat y_i = c + m_1\bar t + m_2\bar s$ (average each term), so $\bar y = c + m_1\bar t + m_2\bar s$ , which says the plane passes through $(\bar t, \bar s, \bar y)$ .

prove

Let $A$ be any real $m \times n$ matrix (the columns may be dependent) and let $\lambda > 0$ . Prove that $A^{\mathsf T}A + \lambda I$ is invertible.

Solution

Suppose $(A^{\mathsf T}A + \lambda I)x = O$ . Multiply on the left by $x^{\mathsf T}$ : $x^{\mathsf T}A^{\mathsf T}Ax + \lambda x^{\mathsf T}x = 0$ , that is $\lVert Ax \rVert^2 + \lambda\lVert x \rVert^2 = 0$ . Both terms are non-negative and $\lambda > 0$ , so both are zero, and $\lVert x \rVert = 0$ gives $x = O$ . So the null space of this square matrix is $\{O\}$ , and by the argument in the uniqueness theorem it is invertible.

stretch

Let $P_1$ and $P_2$ be $m \times m$ matrices that are both symmetric and idempotent and have the same range $W$ . Prove $P_1 = P_2$ . (So the orthogonal projection onto a subspace is a single, well-defined matrix.)

Solution

Fix $b$ . By the prove-problem on symmetric idempotent matrices, $b - P_1b$ is perpendicular to every vector of $W$ , and $P_1b$ is in $W$ because $W$ is the range of $P_1$ . The same holds for $P_2$ . By the theorem ‘Perpendicular error means closest point’, each of $P_1b$ and $P_2b$ is the unique closest point of $W$ to $b$ . A unique closest point is one point, so $P_1b = P_2b$ . This holds for every $b$ , so $P_1 = P_2$ .

stretch

Add one more column $a$ to $A$ , giving $A' = [A \mid a]$ . Prove that the smallest possible sum of squared residuals for $A'$ is at most the smallest possible for $A$ . (More features can never raise the training error, which is why training error is a poor judge of a model.)

Solution

Let $\hat x$ be a least squares solution for $A$ with minimum squared error $E = \lVert b - A\hat x \rVert^2$ . The vector $x' = (\hat x, 0)$ gives $A'x' = A\hat x + 0 \cdot a = A\hat x$ , so it has the same squared error $E$ . The least squares solution for $A'$ minimizes over all vectors, including this one, so its minimum is at most $E$ . The inequality can be strict: adding a column enlarges the column space, so there are more points to be close to.

Summary

When $Ax = b$ has no solution, ask for the $x$ whose $Ax$ is closest to $b$ . The set of all $Ax$ is the column space, so this is projection onto the column space (2.9).
The closest point is the one whose error $r = b - A\hat x$ is perpendicular to the column space. Proof: Pythagoras on the triangle formed by $b$ , the closest point and any other point.
“Perpendicular to every column” is $A^{\mathsf T}r = O$ , which is the normal equations $A^{\mathsf T}A\,\hat x = A^{\mathsf T}b$ . They are equivalent to being a least squares solution.
$A^{\mathsf T}A$ has the same null space as $A$ (because $x^{\mathsf T}A^{\mathsf T}Ax = \lVert Ax \rVert^2$ ). So it is invertible exactly when the columns are independent, and then $\hat x = (A^{\mathsf T}A)^{-1}A^{\mathsf T}b$ is unique.
The projection matrix $P = A(A^{\mathsf T}A)^{-1}A^{\mathsf T}$ satisfies $P^2 = P$ , $P^{\mathsf T} = P$ and $\operatorname{tr}P = n$ . It sends $b$ to the fitted values.
Computers use QR, $A = QR$ : solve $R\hat x = Q^{\mathsf T}b$ and $P = QQ^{\mathsf T}$ . Forming $A^{\mathsf T}A$ squares the condition number and can destroy the problem.
The sum of squared errors of the best fit is $\lVert b \rVert^2 - \lVert p \rVert^2$ .
Regression, ridge regression, GPS fixes, sensor calibration and Hooke’s law are all this one picture: a data vector, a subspace, and a perpendicular.

Watch and read

The same ideas explained by someone else. Use these when an explanation here does not land.

2.10 Least squares and regression as projection

The picture

The idea, in words

Definitions and theorems

Least squares solution

Perpendicular error means closest point

Normal equations

Least squares solutions solve the normal equations

A and A-transpose A have the same null space

Independent columns give a unique least squares solution

The projection matrix

Projection matrix

The projection matrix projects

QR: how computers do it (B)

QR from Gram–Schmidt

Examples and non-examples

Where it is used

Traps

Problems

Summary

Watch and read

Where this sits

Comes from

Leads to

Same idea elsewhere