Lecture 12: Least Squares & Linear Regression

The idea in one sentence

When no exact answer exists, take the one whose leftover error is as small as it can be, and that answer is the projection from the last chapter read as a list of weights.

Five equations, two knobs

Here is the measurements table one last time, with all five rows asking for the same pair of coefficients $u_1$ and $u_2$ :

\begin{aligned} 5u_1 + 2u_2 &= 3 \\ 3u_1 + 4u_2 &= 8 \\ 2u_1 + 8u_2 &= 7 \\ 1u_1 + 5u_2 &= 9 \\ 4u_1 + 7u_2 &= 2 \end{aligned}

Lecture 11 showed there is no answer. Rank 2 for the matrix of coefficients, rank 3 once the target joins as a third column, so the target points somewhere the two columns cannot go. That is the normal situation for data, and a chapter that stops there has taught you nothing you can use. This one solves it anyway.

The same numbers, read two ways

Before solving, look at the failure from both sides, because the two views explain each other.

Two features span a plane, and a plane inside five-dimensional space is hard to draw. So for the pictures, use one feature and an intercept: keep the weight column and the target column. Write $w_i$ for the weight in row $i$ and $b_i$ for its target, and fit the two unknowns $C$ and $D$ in

b_i \;\approx\; C + D\,w_i .

Every argument below works for any number of columns. The two-feature answer arrives at the end of the next section.

Now the two readings. Read a row as a point. Row one says the pair $(w_1, b_1) = (5, 3)$ was measured, and a candidate line either passes near it or does not. Five rows are five dots on a page, and the unknowns $C$ and $D$ pick out one line among all lines. That is data space.

Read the same row as a constraint instead. Row one demands $C + 5D = 3$ , and the pairs $(C, D)$ that satisfy it form a straight line in a plane whose axes are the unknowns. Five rows are five lines, and a candidate line in data space is a single dot here. That is weight space.

Fig. 1

The same five rows in two rooms. On the left each row is a dot and a candidate model is a line. On the right each row is a line and a candidate model is a dot. The bold line and the bold dot are the same candidate. No dot sits on all five lines, so no line passes through all five dots.

The two rooms hold the same information turned inside out. A line over there is a point over here. Fitting the data exactly would mean finding one point that sits on all five weight-space lines at once, and the picture says there is no such point. Every pair of lines crosses somewhere different, which is another way of saying the first two rows and the third row want different coefficients.

Choosing what best means

No point sits on all five lines, so pick a rule for scoring the ones that do not. Write the model as $Ax = b$ , where $A$ holds the columns, $x$ holds the unknown coefficients, and $b$ holds the targets. Any candidate $x$ leaves a leftover $b - Ax$ , one number per row, and the score is the total of their squares:

E(x) \;=\; \norm{b - Ax}^2 \;=\; \sum_i \left( b_i - (Ax)_i \right)^2 .

Squaring is a choice, and absolute values would also measure the leftovers. Two reasons pick squares. The first is geometric. That sum is the ordinary squared distance from the point $b$ to the point $Ax$ , so making it small means finding the reachable point closest to $b$ , which is the projection Lecture 11 already built. The second is practical: squares make the problem linear, so a formula exists. Absolute values give a perfectly sensible score with no closed form.

Squaring has a cost. A leftover of 4 counts as much as sixteen leftovers of 1, so one bad outlier can drag the whole fit toward itself. Every model fit with squares carries that sensitivity.

The score has a second reading, back in weight space. Write $a_i$ for row $i$ of $A$ , so row $i$ asks for $a_i\T x = b_i$ and draws one line in the picture of unknowns. The distance from a candidate point $x$ to that line is $\lvert b_i - a_i\T x \rvert / \norm{a_i}$ , which is the leftover divided by the length of the row. Squaring and adding gives

E(x) \;=\; \sum_i \norm{a_i}^2 \, d_i(x)^2 ,

where $d_i(x)$ is that distance. So the answer settles at the point closest to all five lines at once, with each line pulling in proportion to the squared length of its row. A row of big numbers is a line that insists.

So the least-squares answer $\hat{x}$ is the coefficient list of the projection, and $p = A\hat{x}$ is the shadow of the data on what the model can reach. Nothing new has happened yet. The coefficients now have a meaning, and the leftover has a name: the residual.

The normal equations, twice

Two roads lead to the same equation. Take both, because each answers a question the other leaves open.

By perpendicularity. This is the argument from Lecture 11, one line long. At the closest point, the residual $e = b - A\hat{x}$ is perpendicular to every column of $A$ , and stacking those conditions is what $A\T$ does:

A\T (b - A\hat{x}) = 0 \qquad\Longrightarrow\qquad A\T A\, \hat{x} = A\T b .

That pair of equations is called the normal equations, from the old word normal meaning perpendicular. They are the whole of least squares.

By algebra. The first road shows that $\hat{x}$ is a special point. It does not show that $\hat{x}$ beats every other point, and it does not show that no second point ties with it. Expand the score instead. For any $x$ ,

E(x) = (b - Ax)\T (b - Ax) = b\T b - 2x\T A\T b + x\T A\T A x .

Let $\hat{x}$ be a solution of $A\T A \hat{x} = A\T b$ and write any other candidate as $x = \hat{x} + d$ . Substitute and collect terms. The cross terms combine into a single bracket:

E(\hat{x} + d) = E(\hat{x}) + 2\,d\T \left( A\T A \hat{x} - A\T b \right) + d\T A\T A d .

The bracket is zero, by the choice of $\hat{x}$ . And $d\T A\T A d = (Ad)\T (Ad) = \norm{Ad}^2$ . So

E(\hat{x} + d) = E(\hat{x}) + \norm{Ad}^2 \;\ge\; E(\hat{x}) .

Every other candidate scores worse by exactly $\norm{Ad}^2$ , the squared length of the change that move makes inside the column space. Equality needs $Ad = 0$ , and with independent columns that forces $d = 0$ . So $\hat{x}$ is the winner and it is the only one. Lecture 11 proved $N(A\T A) = N(A)$ , which is the same independence condition seen from the other end, and it is what lets you write $\hat{x} = (A\T A)^{-1} A\T b$ when you want the formula spelled out.

The answer, in numbers

For the line fit, $A$ has a column of ones and the weight column. Then $A\T A$ counts, sums, and sums of squares:

A\T A = \begin{bmatrix} 5 & 15 \\ 15 & 55 \end{bmatrix}, \qquad A\T b = \begin{bmatrix} 29 \\ 70 \end{bmatrix}.

The 5 is the number of rows, the 15 is $\sum w_i$ , the 55 is $\sum w_i^2$ , the 29 is $\sum b_i$ , and the 70 is $\sum w_i b_i$ . The determinant is $275 - 225 = 50$ , and solving gives $C = 10.9$ and $D = -1.7$ exactly.

Fig. 2

The line that least squares picks, with every residual drawn. The circles on the line are what the model predicts, the filled dots are what was measured, and the dashed gaps between them are what the fit could not explain.

The five fitted values are $2.4$ , $5.8$ , $7.5$ , $9.2$ , $4.1$ , and the residuals are $0.6$ , $2.2$ , $-0.5$ , $-0.2$ , $-2.1$ . Squared and added they come to $9.9$ . Try any other line and you get a bigger number, which the lab below lets you test for yourself.

Two facts about those residuals are not luck.

Fig. 3

Perpendicularity in numbers. The first row of $A\T$ is the intercept column, so its dot product with the residual is the plain sum. The second row is the weight column, so its dot product is the weighted sum. Both come out zero.

They add to zero, and their weighted total $\sum w_i e_i$ is zero as well. Both are just $A\T e = 0$ read one row at a time. The first row of $A\T$ is the row of ones, so its dot product with $e$ is the plain sum. The second row is $w$ , so its dot product is the weighted sum. Any model with an intercept column has residuals summing to zero, and that comes from perpendicularity alone. The data has no say in it.

It also explains why no better line exists. To improve the fit you would tilt or shift the line, which changes the prediction by some combination of the two columns. But the residual leans on neither column, so every such change moves you away from $b$ .

Watch what one move costs. Hold the slope at $-1.7$ and slide the intercept away from $10.9$ by an amount $t$ , which is the direction $d = (1, 0)$ . The algebra above says the score becomes $9.9 + \norm{Ad}^2 t^2$ , and $\norm{Ad}^2$ for this $d$ is the top left entry of $A\T A$ , which is 5. So one full unit of shift costs exactly 5.

Fig. 4

One straight walk away from the answer. Sliding the intercept by $t$ raises the score by $5t^2$ , so the bottom of this parabola sits exactly where the normal equations put it.

Here is the whole story as one surface. Score every possible pair $(C, D)$ by its total squared error, and the result is a bowl. The algebra above already told you its shape: $E(\hat{x} + d) = 9.9 + d\T A\T A d$ , so the level sets are ellipses centered on the answer.

Fig. 5

The same story in two dimensions. Each ellipse collects the lines with one fixed total squared error, and the dot at the center is the winner. The long thin shape means the score hardly changes along one direction, which is the direction a fitted model is least sure about.

The ellipses are long and thin, and that ratio matters later. It says the score barely changes along one direction and changes fast along another, so a fit can be confidently wrong about the individual coefficients while being right about their combination.

Try it. Move the data, then try to beat the fit

your intercept C = 8.0your slope D = -0.5

best line: 10.9 −1.7w
its error: 9.90
your error: 21.30

Drag any dot to change the data. Green is the line the normal equations name, warm is yours, and the blue springs are the residuals it has to live with.

Now the original two-feature question, with $w$ and $h$ as the columns and no intercept:

A\T A = \begin{bmatrix} 55 & 71 \\ 71 & 158 \end{bmatrix}, \qquad A\T b = \begin{bmatrix} 70 \\ 153 \end{bmatrix}, \qquad \det = 3649 .

Solving gives $\hat{x} = (197/3649,\; 3445/3649)$ , which is about $(0.0540,\; 0.9441)$ . The fitted values are $2.16$ , $3.94$ , $7.66$ , $4.77$ , $6.82$ , and the squared residuals add to $58.77$ . The code below computes it, checks it against np.linalg.lstsq, and confirms that the residual dots to zero against both columns. The fit is poor, and it should be: the numbers in the table were never generated by a rule. Least squares always returns an answer. Reading whether that answer means anything is your job, and the formula will not do it for you.

python

import numpy as np

w = np.array([5., 3., 2., 1., 4.])
h = np.array([2., 4., 8., 5., 7.])
b = np.array([3., 8., 7., 9., 2.])

A = np.column_stack([w, h])
print(A.T @ A)                      # [[ 55.  71.] [ 71. 158.]]
print(A.T @ b)                      # [ 70. 153.]
x_hat = np.linalg.solve(A.T @ A, A.T @ b)
print(np.round(x_hat, 4))           # [0.054  0.9441]
print(np.round(np.linalg.lstsq(A, b, rcond=None)[0], 4))   # the same numbers
e = b - A @ x_hat
print(bool(abs(A.T @ e).max() < 1e-6))     # True   e leans on neither column
print(round(float(e @ e), 4))              # 58.7745

M = np.column_stack([np.ones(5), w])      # intercept column, then one feature
print(M.T @ M)                      # [[ 5. 15.] [15. 55.]]
print(M.T @ b)                      # [29. 70.]
c, d = np.linalg.solve(M.T @ M, M.T @ b)
print(round(c, 10), round(d, 10))   # 10.9 -1.7
r = b - M @ np.array([c, d])
print(np.round(M @ np.array([c, d]), 4))   # [2.4 5.8 7.5 9.2 4.1]
print(np.round(r, 4))                      # [ 0.6  2.2 -0.5 -0.2 -2.1]
print(bool(abs(r.sum()) < 1e-6), bool(abs(w @ r) < 1e-6))    # True True
print(round(float(r @ r), 4))              # 9.9
for C, D in [(11.0, -1.7), (10.9, -1.6), (10.0, -1.5)]:
    print(C, D, round(float(((b - (C + D * w)) ** 2).sum()), 4))
# 11.0 -1.7 9.95 / 10.9 -1.6 10.45 / 10.0 -1.5 10.75: all worse than 9.9

python

import torch

w = torch.tensor([5., 3., 2., 1., 4.])
h = torch.tensor([2., 4., 8., 5., 7.])
b = torch.tensor([3., 8., 7., 9., 2.])

A = torch.stack([w, h], dim=1)
print(A.T @ A)                      # [[ 55.  71.] [ 71. 158.]]
print(A.T @ b)                      # [ 70. 153.]
x_hat = torch.linalg.solve(A.T @ A, A.T @ b)
print(torch.round(x_hat, decimals=4))                        # [0.054  0.9441]
print(torch.round(torch.linalg.lstsq(A, b).solution, decimals=4))  # the same
e = b - A @ x_hat
print(bool(abs(A.T @ e).max() < 1e-4))     # True   e leans on neither column
print(round(float(e @ e), 4))              # 58.7745

M = torch.stack([torch.ones(5), w], dim=1)  # intercept, then one feature
print(M.T @ M)                      # [[ 5. 15.] [15. 55.]]
print(M.T @ b)                      # [29. 70.]
cd = torch.linalg.solve(M.T @ M, M.T @ b)
print(round(float(cd[0]), 4), round(float(cd[1]), 4))    # 10.9 -1.7
r = b - M @ cd
print(torch.round(M @ cd, decimals=4))     # [2.4 5.8 7.5 9.2 4.1]
print(torch.round(r, decimals=4))          # [ 0.6  2.2 -0.5 -0.2 -2.1]
print(bool(abs(r.sum()) < 1e-4), bool(abs(w @ r) < 1e-4))    # True True
print(round(float(r @ r), 4))              # 9.9
for C, D in [(11.0, -1.7), (10.9, -1.6), (10.0, -1.5)]:
    print(C, D, round(float(((b - (C + D * w)) ** 2).sum()), 4))
# 11.0 -1.7 9.95 / 10.9 -1.6 10.45 / 10.0 -1.5 10.75: all worse than 9.9

python

import jax.numpy as jnp

w = jnp.array([5., 3., 2., 1., 4.])
h = jnp.array([2., 4., 8., 5., 7.])
b = jnp.array([3., 8., 7., 9., 2.])

A = jnp.column_stack([w, h])
print(A.T @ A)                      # [[ 55.  71.] [ 71. 158.]]
print(A.T @ b)                      # [ 70. 153.]
x_hat = jnp.linalg.solve(A.T @ A, A.T @ b)
print(jnp.round(x_hat, 4))          # [0.054  0.9441]
print(jnp.round(jnp.linalg.lstsq(A, b)[0], 4))    # the same numbers
e = b - A @ x_hat
print(bool(abs(A.T @ e).max() < 1e-4))     # True   e leans on neither column
print(round(float(e @ e), 4))              # 58.7745

M = jnp.column_stack([jnp.ones(5), w])    # intercept column, then one feature
print(M.T @ M)                      # [[ 5. 15.] [15. 55.]]
print(M.T @ b)                      # [29. 70.]
cd = jnp.linalg.solve(M.T @ M, M.T @ b)
print(round(float(cd[0]), 4), round(float(cd[1]), 4))    # 10.9 -1.7
r = b - M @ cd
print(jnp.round(M @ cd, 4))                # [2.4 5.8 7.5 9.2 4.1]
print(jnp.round(r, 4))                     # [ 0.6  2.2 -0.5 -0.2 -2.1]
print(bool(abs(r.sum()) < 1e-4), bool(abs(w @ r) < 1e-4))    # True True
print(round(float(r @ r), 4))              # 9.9
for C, D in [(11.0, -1.7), (10.9, -1.6), (10.0, -1.5)]:
    print(C, D, round(float(((b - (C + D * w)) ** 2).sum()), 4))
# 11.0 -1.7 9.95 / 10.9 -1.6 10.45 / 10.0 -1.5 10.75: all worse than 9.9

python

import tensorflow as tf

w = tf.constant([5., 3., 2., 1., 4.])
h = tf.constant([2., 4., 8., 5., 7.])
b = tf.constant([3., 8., 7., 9., 2.])
mv = tf.linalg.matvec
gap = lambda z: float(tf.reduce_max(tf.abs(z)))
solve = lambda G, v: tf.linalg.solve(G, v[:, None])[:, 0]

A = tf.stack([w, h], axis=1)
At = tf.transpose(A)
print(At @ A)                       # [[ 55.  71.] [ 71. 158.]]
print(mv(At, b))                    # [ 70. 153.]
x_hat = solve(At @ A, mv(At, b))
print(tf.round(x_hat * 1e4) / 1e4)  # [0.054  0.9441]
print(tf.round(tf.linalg.lstsq(A, b[:, None])[:, 0] * 1e4) / 1e4)  # the same
e = b - mv(A, x_hat)
print(gap(mv(At, e)) < 1e-4)        # True   e leans on neither column
print(round(float(tf.tensordot(e, e, 1)), 4))     # 58.7745

M = tf.stack([tf.ones(5), w], axis=1)     # intercept column, then one feature
Mt = tf.transpose(M)
print(Mt @ M)                       # [[ 5. 15.] [15. 55.]]
print(mv(Mt, b))                    # [29. 70.]
cd = solve(Mt @ M, mv(Mt, b))
print(round(float(cd[0]), 4), round(float(cd[1]), 4))    # 10.9 -1.7
r = b - mv(M, cd)
print(tf.round(mv(M, cd) * 1e4) / 1e4)     # [2.4 5.8 7.5 9.2 4.1]
print(tf.round(r * 1e4) / 1e4)             # [ 0.6  2.2 -0.5 -0.2 -2.1]
print(gap(tf.reduce_sum(r)) < 1e-4, gap(tf.tensordot(w, r, 1)) < 1e-4)
print(round(float(tf.tensordot(r, r, 1)), 4))     # 9.9
for C, D in [(11.0, -1.7), (10.9, -1.6), (10.0, -1.5)]:
    print(C, D, round(float(tf.reduce_sum((b - (C + D * w)) ** 2)), 4))
# 11.0 -1.7 9.95 / 10.9 -1.6 10.45 / 10.0 -1.5 10.75: all worse than 9.9

More columns, and the trap that comes with them

Nothing in the derivation cared what the columns held. They can hold anything computed from the input, and the fit is still linear in the unknowns, which is all least squares ever needed. Put $1$ , $t$ , and $t^2$ in three columns and you are fitting a parabola with the same machinery.

Take five readings of something that bends: $t = -2, -1, 0, 1, 2$ with values $7, 1, 0, 3, 9$ . With two columns the best line is $4 + 0.6t$ and the error is $56.4$ . Add the $t^2$ column and the best parabola is $2t^2 + 0.6t$ , with error $0.4$ .

Fig. 6

Five readings that bend. The dashed line is the best a straight fit can do, and the solid curve is the best parabola, found by the same normal equations with one more column.

A curve is the same projection as before, onto a column space that now happens to contain curved shapes.

Adding columns can only lower the error, since the old fit is still available with a zero coefficient on the new column. That sounds like good news until you follow it to the end. The measurements table has five rows and five distinct weights, so five columns of powers reach every possible target vector, and the error drops to exactly zero.

Fig. 7

Five points and five columns of powers. The curve hits every point exactly, so the error is zero, and it tells you nothing about a weight it has not seen.

The curve through all five points is not a discovery. It is the data written in a different alphabet. It dips below $0$ between points that never left the range $2$ to $9$ , and asked about a weight of $6$ it answers $49$ with total confidence. The error went to zero because the model ran out of things to be wrong about.

This is the training step

Nothing in this chapter mentioned learning, and yet this is what fitting a linear model is. The columns are features, the unknowns are weights, the score is the loss, and the answer is the trained model. When a library fits a linear regression it forms the normal equations, or it does something numerically better that lands in the same place.

Gradient descent lands there too. The gradient of the score is

\nabla E(x) = 2 A\T (Ax - b),

and setting it to zero is the normal equations again, rearranged. So descent walks downhill on the bowl from Figure 5 toward the point the normal equations name. It gets there slowly when the bowl is long and thin, because the step that suits one direction overshoots the other, and that zigzag has a name and a cure that arrive in Lecture 15.

Fig. 8

Gradient descent on the same bowl, started at the open circle. Every step points straight downhill, and the steepest direction points across the valley, so the walk zigzags toward the point the normal equations name in one shot.

python

import numpy as np

t = np.array([-2., -1., 0., 1., 2.])
y = np.array([7., 1., 0., 3., 9.])
line = np.column_stack([np.ones(5), t])
para = np.column_stack([np.ones(5), t, t ** 2])
c_line = np.linalg.solve(line.T @ line, line.T @ y)
c_para = np.linalg.solve(para.T @ para, para.T @ y)
sse = lambda X, c: round(float(((y - X @ c) ** 2).sum()), 2)
print(np.round(c_line, 4), sse(line, c_line))
print(np.round(c_para, 4), sse(para, c_para))
# [4.  0.6] 56.4        a line cannot bend
# [0.  0.6 2. ] 0.4     one more column, and it can

w = np.array([5., 3., 2., 1., 4.])
b = np.array([3., 8., 7., 9., 2.])
for k in (2, 3, 5):                     # k columns of powers of w
    V = np.vander(w, k)
    c = np.linalg.lstsq(V, b, rcond=None)[0]
    print(k, round(float(((b - V @ c) ** 2).sum()), 2))
# 2 9.9 / 3 9.83 / 5 0.0   five columns, five rows, nothing left to get wrong
print(f"{float(np.linalg.cond(np.vander(w, 5))):.1e}")     # 2.6e+04

A = np.column_stack([np.ones(5), w])   # gradient descent finds the same point
x = np.zeros(2)
for step in range(3000):
    x = x - 0.008 * (2 * A.T @ (A @ x - b))
print(np.round(x, 4))                   # [10.9 -1.7]

python

import torch

t = torch.tensor([-2., -1., 0., 1., 2.])
y = torch.tensor([7., 1., 0., 3., 9.])
line = torch.stack([torch.ones(5), t], dim=1)
para = torch.stack([torch.ones(5), t, t ** 2], dim=1)
c_line = torch.linalg.solve(line.T @ line, line.T @ y)
c_para = torch.linalg.solve(para.T @ para, para.T @ y)
sse = lambda X, c: round(float(((y - X @ c) ** 2).sum()), 2)
print(torch.round(c_line, decimals=4), sse(line, c_line))
print(torch.round(c_para, decimals=4), sse(para, c_para))
# [4.  0.6] 56.4        a line cannot bend
# [0.  0.6 2. ] 0.4     one more column, and it can

w = torch.tensor([5., 3., 2., 1., 4.])
b = torch.tensor([3., 8., 7., 9., 2.])
for k in (2, 3, 5):                     # k columns of powers of w
    V = torch.stack([w ** j for j in range(k - 1, -1, -1)], dim=1)
    c = torch.linalg.lstsq(V, b, driver='gelsd').solution
    print(k, round(float(((b - V @ c) ** 2).sum()), 2))
# 2 9.9 / 3 9.83 / 5 0.0   five columns, five rows, nothing left to get wrong
V5 = torch.stack([w ** j for j in range(4, -1, -1)], dim=1)
print(f"{float(torch.linalg.cond(V5)):.1e}")               # 2.6e+04

A = torch.stack([torch.ones(5), w], dim=1)  # gradient descent, same point
x = torch.zeros(2)
for step in range(3000):
    x = x - 0.008 * (2 * A.T @ (A @ x - b))
print(torch.round(x, decimals=4))       # [10.9 -1.7]

python

import jax.numpy as jnp

t = jnp.array([-2., -1., 0., 1., 2.])
y = jnp.array([7., 1., 0., 3., 9.])
line = jnp.column_stack([jnp.ones(5), t])
para = jnp.column_stack([jnp.ones(5), t, t ** 2])
c_line = jnp.linalg.solve(line.T @ line, line.T @ y)
c_para = jnp.linalg.solve(para.T @ para, para.T @ y)
sse = lambda X, c: round(float(((y - X @ c) ** 2).sum()), 2)
print(jnp.round(c_line, 4), sse(line, c_line))
print(jnp.round(c_para, 4), sse(para, c_para))
# [4.  0.6] 56.4        a line cannot bend
# [0.  0.6 2. ] 0.4     one more column, and it can

w = jnp.array([5., 3., 2., 1., 4.])
b = jnp.array([3., 8., 7., 9., 2.])
for k in (2, 3, 5):                     # k columns of powers of w
    V = jnp.vander(w, k)
    c = jnp.linalg.lstsq(V, b)[0]
    print(k, round(float(((b - V @ c) ** 2).sum()), 2))
# 2 9.9 / 3 9.83 / 5 0.0   five columns, five rows, nothing left to get wrong
print(f"{float(jnp.linalg.cond(jnp.vander(w, 5))):.1e}")   # 2.6e+04

A = jnp.column_stack([jnp.ones(5), w])  # gradient descent, same point
x = jnp.zeros(2)
for step in range(3000):
    x = x - 0.008 * (2 * A.T @ (A @ x - b))
print(jnp.round(x, 4))                  # [10.9 -1.7]

python

import tensorflow as tf

mv = tf.linalg.matvec
solve = lambda G, v: tf.linalg.solve(G, v[:, None])[:, 0]
fit = lambda X, v: solve(tf.transpose(X) @ X, mv(tf.transpose(X), v))
sse = lambda X, c, v: round(float(tf.reduce_sum((v - mv(X, c)) ** 2)), 2)

t = tf.constant([-2., -1., 0., 1., 2.])
y = tf.constant([7., 1., 0., 3., 9.])
line = tf.stack([tf.ones(5), t], axis=1)
para = tf.stack([tf.ones(5), t, t ** 2], axis=1)
c_line, c_para = fit(line, y), fit(para, y)
print(tf.round(c_line * 1e4) / 1e4, sse(line, c_line, y))
print(tf.round(c_para * 1e4) / 1e4, sse(para, c_para, y))
# [4.  0.6] 56.4        a line cannot bend
# [0.  0.6 2. ] 0.4     one more column, and it can

w = tf.constant([5., 3., 2., 1., 4.])
b = tf.constant([3., 8., 7., 9., 2.])
for k in (2, 3, 5):                     # k columns of powers of w
    V = tf.stack([w ** float(j) for j in range(k - 1, -1, -1)], axis=1)
    c = tf.linalg.lstsq(V, b[:, None], fast=False)[:, 0]
    print(k, sse(V, c, b))
# 2 9.9 / 3 9.83 / 5 0.0   five columns, five rows, nothing left to get wrong
V5 = tf.stack([w ** float(j) for j in range(4, -1, -1)], axis=1)
s = tf.linalg.svd(V5, compute_uv=False)
print(f"{float(s[0] / s[-1]):.1e}")     # 2.6e+04

A = tf.stack([tf.ones(5), w], axis=1)  # gradient descent finds the same point
x = tf.zeros(2)
for step in range(3000):
    x = x - 0.008 * (2 * mv(tf.transpose(A), mv(A, x) - b))
print(tf.round(x * 1e4) / 1e4)          # [10.9 -1.7]

Where this is going

Least squares is the first honest answer this course has given to a question with no answer. It came from one geometric fact, that the shortest distance to a subspace is measured along a perpendicular, and it produced a formula you can compute.

Two threads run onward. The first is numerical. Forming $A\T A$ is the wrong way to compute $\hat{x}$ , and Lecture 13 shows the right way by building orthonormal columns first, which turns the normal equations into back substitution and stops the condition number from being squared. The second is what happens when the columns are dependent, so $A\T A$ has no inverse and the best answer is a whole flat set of them. Choosing sensibly among that set is what the pseudoinverse does, and Part VI cashes it in with the singular value decomposition.

Check yourself

Fit $b = C + Dt$ to the three points $(0, 6)$ , $(1, 0)$ , $(2, 0)$ . Build $A\T A$ and $A\T b$ by hand, solve, and check the residuals sum to zero.
In the same problem, compute $p = A\hat{x}$ and $e = b - p$ , and verify $e$ is perpendicular to both columns of $A$ .
Why does an intercept column force the residuals to add to zero? Answer using $A\T e = 0$ and one sentence.
Take the line fit in Figure 2 and change the target of the third row from 7 to 8. Predict, before computing, whether $D$ goes up or down, then check with the lab.
A model has 5 data points and 5 independent columns. What is the least-squares error, what is the projection matrix $P$ , and why is neither fact worth celebrating?
Redraw Figure 1 from memory, both panels and their labels. The drawing is the test.