Part II · Chapter 8

Linear regression

By the end of this page you can derive the least-squares line from scratch — twice, by completing the square and by calculus — write down the normal equations $\theta^\star=(X^\top X)^{-1}X^\top y$ and say exactly when they break.

Essentials

This chapter fits a line (then a hyperplane) to real-valued labels by minimising squared error, and solves it exactly — in closed form.

Everything below is the explanation and the practice. If you are short on time, read this box, then go to the problems.

Why anybody fits a line

The least glamorous possible situation: $\X = \R$, $\Y = \R$, and $N$ training pairs

$$\D = \{(x^{(1)},y^{(1)}),\ (x^{(2)},y^{(2)}),\ \dots,\ (x^{(N)},y^{(N)})\}.$$

Notation, first time: the superscript in $x^{(i)}$ is an index, not a power — $x^{(3)}$ is the third training input, not $x$ cubed. The parentheses exist precisely so it can't be read as an exponent, and later $x^{(i)}_j$ will mean "the $j$-th coordinate of the $i$-th sample".

The professor's running example is selling a used car: eight Pontiac Vibe listings, with mileage $x$ and asking price $y$.

mileage $x$ (miles)25,00034,00045,00070,00093,000110,000125,000160,000
price $y$ (USD)23,00028,00022,00022,00012,5005,5005,0004,000

Plot it and the eye draws a line through the cloud. This chapter is the machine version of that reflex.

Figure 1 · The data, the best line, and the residuals it leaves behind
What to notice: the line does not pass through a single data point. Every point contributes a residual — the vertical gap $y^{(i)}-h(x^{(i)})$, drawn as a segment. Fitting a line means choosing the line that makes these segments collectively small. Vertical, not perpendicular: we are predicting $y$ from $x$, so only error in $y$ counts.

Two reasons to care. First, a line is the simplest family that outputs a number and actually responds to the input: two knobs, a closed-form optimum, and over a limited range almost everything looks linear. Second, the last layer of a neural network is a linear model: freeze a pretrained network's features $z^{(i)} = \phi(x^{(i)}) \in \R^d$ and retrain only the final weights, and the problem you are left with is $\argmin_{w,b}\frac1N\sum_i (y^{(i)} - (z^{(i)\top}w + b))^2$ — exactly this chapter's problem with $z$ in place of $x$. That is what "linear probing" means.

The hypothesis class: two numbers

A learning problem is fixed by a hypothesis class $\H$ and a loss (see the ML formulation chapter). Here:

Definition

$\H$ is the set of all functions $\R \to \R$ that can be written as a line, each pinned down by two real parameters $w \in \R$ and $b \in \R$:

$$h_{w,b}(x) \;=\; x\cdot w + b.$$

$w$ is the slope — units of $y$ per unit of $x$. $b$ is the intercept or bias — the value of $h$ at $x=0$. (Algebra class wrote $y = mx+c$; same line.)

Because a hypothesis is determined by $(w,b)$, searching over $\H$ is searching over $\R^2$ — the move every parametric model in this course makes.

Common mistake

"$h(x)=wx+b$ is a linear function." In the linear algebra sense it is not: linearity requires $f(\alpha x + \beta v)=\alpha f(x)+\beta f(v)$, which forces $f(0)=0$, and $h(0)=b \neq 0$. A line-plus-constant is affine.

The twist that matters more: which variable? $h(x)=wx+b$ is affine in $x$, but with the data held fixed, $f(w,b)=wx+b$ is genuinely linear in $(w,b)$. In ML the parameters are the variables, which is why $f(x)=w_2x^2+w_1x+b$ — a quadratic in $x$ — still counts as a "linear model". Hold that thought until feature engineering.

ERM with squared loss

We need one number saying how bad a candidate line is. For each point the error is $\text{err}^{(i)} = y^{(i)} - \hat y^{(i)}$ with $\hat y^{(i)} = wx^{(i)}+b$. That is $N$ numbers; the pre-reading tries three ways to collapse them to one:

So the loss is $\ell(y,\hat y) = (y-\hat y)^2$ — squared loss — and ERM reads:

Definition $$\argmin_{h \in \H}\ \frac1N\sum_{i=1}^{N}\Big(y^{(i)} - h\big(x^{(i)}\big)\Big)^2$$

Over all lines $h$ in our class, find the one minimising the average over the $N$ training points of the squared gap between $y^{(i)}$ and $h(x^{(i)})$. The $\frac1N$ makes it a mean, so the number does not grow just because you collected more data.

Why squaring is the natural choice for real-valued labels: it removes the sign without removing smoothness, so the calculus of Chapter 9 applies and a closed form exists; it punishes big misses disproportionately (double the error, quadruple the contribution); and its unconstrained minimiser is the mean — $\argmin_c \frac1N\sum_i(y^{(i)}-c)^2 = \bar y$, whereas absolute error's answer is the median. We often report $\text{RMSE}=\sqrt{\text{MSE}}$ because it lives in the units of $y$: for our eight cars the best line has $\text{RMSE} \approx \$3{,}156$, i.e. "typically off by about three thousand dollars".

Common mistake

Squared loss' sensitivity to big errors is a feature and a bug. One mistyped listing — a $\$250{,}000$ Vibe — drags the fitted line noticeably, because that single point contributes an enormous squared term. Least squares is not robust to outliers; absolute error is, and robustness is exactly what we traded for smoothness and a closed form.

You fit a line and the eight residuals are $+400, -400, +600, -600, +100, -100, +9000, -9000$. Which statement is true?
Why: The errors sum to exactly $0$ — and the line is terrible, since two cars are off by $\$9{,}000$. That is precisely the "Attempt 1 fails" argument: signed errors cancel, so the sum says nothing. MSE does not cancel: $\frac18(400^2+400^2+600^2+600^2+100^2+100^2+9000^2+9000^2) = \frac18(162{,}020{,}000) \approx 2.03\times10^7$, i.e. RMSE $\approx \$4{,}500$.

The derivation, line by line (with $b=0$)

The notes simplify first: set $b=0$, pinning every candidate line to the origin so only the slope is free. It is a scaffold — with one unknown every step is visible — and we restore $b$ in the next section. ERM becomes a one-variable minimisation:

$$\argmin_{w \in \R}\ \frac1N\sum_{i=1}^{N}\Big(y^{(i)} - x^{(i)}w\Big)^2.$$

Step 1 — expand the square

For each $i$, $(y^{(i)} - x^{(i)}w)^2 = y^{(i)2} + x^{(i)2}w^2 - 2\,y^{(i)}x^{(i)}w$, which is $(p-q)^2 = p^2+q^2-2pq$. Sum over $i$ and pull $w$ out of each sum:

$$= \argmin_{w\in\R}\ \frac1N\left[\ \sum_{i=1}^N y^{(i)2} \;+\; w^2\sum_{i=1}^N x^{(i)2} \;-\; 2w\sum_{i=1}^N y^{(i)}x^{(i)}\ \right].$$

Everything in the bracket is a number computed from the data except $w$. Name the three data constants:

$$S \;=\; \sum_{i=1}^N x^{(i)2}, \qquad P \;=\; \sum_{i=1}^N y^{(i)}x^{(i)}, \qquad Q \;=\; \sum_{i=1}^N y^{(i)2}.$$

So we minimise $J(w)=\frac1N\left[Q + Sw^2 - 2Pw\right]$: a parabola in $w$, opening upwards whenever $S>0$, hence with exactly one minimum.

Step 2 — completing the square

Proof

Take the two $w$-dependent terms and factor $S$ out:

$$Sw^2 - 2Pw \;=\; S\left(w^2 - \frac{2P}{S}w\right).$$

Inside the bracket, add and subtract the number that completes the square. Since $w^2 - 2\frac{P}{S}w$ is the start of $\left(w-\frac{P}{S}\right)^2 = w^2 - 2\frac{P}{S}w + \frac{P^2}{S^2}$, we add and subtract $\frac{P^2}{S^2}$:

$$= S\left[\left(w-\frac{P}{S}\right)^2 - \frac{P^2}{S^2}\right] \;=\; S\left(w-\frac{P}{S}\right)^2 - \frac{P^2}{S}.$$

Now push the $S$ inside the square, which means multiplying the bracket by $\sqrt S$ (legal because $S \geq 0$, being a sum of squares):

$$S\left(w-\frac{P}{S}\right)^2 \;=\; \left(\sqrt S\left(w - \frac{P}{S}\right)\right)^2 \;=\; \left(w\sqrt{S} \;-\; \frac{P}{\sqrt S}\right)^2,$$

using $\sqrt S \cdot \frac{P}{S} = \frac{P}{\sqrt S}$. Writing the sums back out, that inner expression is

$$w\sqrt{\sum_{i=1}^N x^{(i)2}} \;-\; \frac{\sum_{i=1}^N y^{(i)}x^{(i)}}{\sqrt{\sum_{i=1}^N x^{(i)2}}},$$

which is precisely the form in the lecture notes. Putting the whole objective back together:

$$J(w) \;=\; \frac1N\left[\left(w\sqrt S - \frac{P}{\sqrt S}\right)^2 \;+\; \underbrace{Q - \frac{P^2}{S}}_{\text{no } w \text{ in here}}\right].$$

Step 3 — why we may throw the constant away

The notes drop $Q-\frac{P^2}{S}$ without comment. It is legal, and it is a move you will make a hundred more times this semester.

Key idea

For any function $f$, any constant $c$, and any constant $\lambda > 0$:

$$\argmin_w\ \big[\lambda f(w) + c\big] \;=\; \argmin_w f(w).$$

Proof. $\lambda f(w)+c \leq \lambda f(w')+c \iff \lambda f(w) \leq \lambda f(w') \iff f(w) \leq f(w')$, the last step because dividing by a positive number preserves the inequality. The two functions order every pair of candidates identically, so they have the same minimisers. $\square$

Adding $c$ slides the graph vertically; multiplying by $\lambda>0$ stretches it. Neither moves the location of the bottom. Two caveats: the $\argmin$ is unchanged but the $\min$ value is not, and $\lambda<0$ would turn the minimum into a maximum. This is also why the $\frac1N$ never affects the answer.

Step 4 — read off the minimiser

After dropping the constant and the $\frac1N$ we are minimising a single squared quantity:

$$\argmin_{w\in\R}\ \left(w\sqrt S - \frac{P}{\sqrt S}\right)^2.$$

A square is always $\ge 0$, so the smallest value possible is $0$ — and it is attained, at the $w$ making the inside vanish. Because $0$ is achievable, that $w$ is the minimiser. (This "lower bound plus achievability" argument is what makes completing the square a proof rather than a trick.)

Set the inside to zero and multiply through by $\sqrt S$:

$$w\sqrt S - \frac{P}{\sqrt S} = 0 \quad\Longleftrightarrow\quad wS - P = 0 \quad\Longleftrightarrow\quad w\sum_{i=1}^N x^{(i)2} - \sum_{i=1}^N y^{(i)}x^{(i)} = 0.$$
Key idea $$\boxed{\;w^\star \;=\; \frac{\sum_{i=1}^N y^{(i)}x^{(i)}}{\sum_{i=1}^N x^{(i)2}}\;}$$

The best-fitting line through the origin is $y = w^\star x$. Numerator: how much $x$ and $y$ move together. Denominator: how much $x$ moves at all. The minimum loss achieved is $\frac1N\left(Q - \frac{P^2}{S}\right)$ — non-zero in general, since the points don't lie on a line.

Step 5 — the calculus cross-check

Proof

Differentiate $J(w) = \frac1N\left[Q + Sw^2 - 2Pw\right]$ term by term. $Q$ is constant so contributes $0$; $\frac{d}{dw}w^2 = 2w$; $\frac{d}{dw}w = 1$:

$$\frac{dJ}{dw} \;=\; \frac1N\left[2Sw - 2P\right].$$

Set it to zero — at an interior minimum of a differentiable function the derivative vanishes:

$$2Sw - 2P = 0 \quad\Longrightarrow\quad w^\star = \frac{P}{S} = \frac{\sum_i y^{(i)}x^{(i)}}{\sum_i x^{(i)2}}.$$

The same answer. And the second derivative is $\frac{d^2J}{dw^2} = \frac{2S}{N} > 0$ whenever $S>0$, confirming it is a minimum and not a maximum or an inflection. $\square$

Intuition

Why show it twice? Completing the square proves a global minimum for free — the objective is literally a square plus a constant — and generalises to matrices, which is what we do below. Setting the derivative to zero only finds a stationary point, so it needs the second-derivative test, but it generalises to messy functions.

Guided problemFit a through-the-origin line to four points, by hand

Training data: $(1,1),\ (2,3),\ (3,4),\ (4,8)$ — that is $x^{(1)}=1,y^{(1)}=1$, and so on, $N=4$. Assume $b=0$. Find $w^\star$, then the MSE it achieves.

Step 1 — What do I actually need to compute?

Only two numbers: $P=\sum_i y^{(i)}x^{(i)}$ and $S=\sum_i x^{(i)2}$. Then $w^\star = P/S$. Nothing else about the data matters — which is itself informative: the fit depends on the data only through these two summaries.

Step 2 — Build the table
$i$$x^{(i)}$$y^{(i)}$$x^{(i)}y^{(i)}$$x^{(i)2}$
111$1\cdot1=1$$1$
223$2\cdot3=6$$4$
334$3\cdot4=12$$9$
448$4\cdot8=32$$16$
totals$P=51$$S=30$
Step 3 — Divide

$$w^\star = \frac{P}{S} = \frac{51}{30} = 1.7.$$

The fitted line is $h(x) = 1.7x$.

Step 4 — What loss does it achieve?

Predictions: $h(1)=1.7$, $h(2)=3.4$, $h(3)=5.1$, $h(4)=6.8$. Residuals $y-\hat y$: $1-1.7=-0.7$; $3-3.4=-0.4$; $4-5.1=-1.1$; $8-6.8=+1.2$.

$$\text{MSE} = \tfrac14\left(0.49 + 0.16 + 1.21 + 1.44\right) = \tfrac14(3.30) = 0.825.$$

Cross-check with the formula $\frac1N(Q - P^2/S)$: $Q = 1+9+16+64 = 90$, so $\frac14\left(90 - \frac{51^2}{30}\right) = \frac14\left(90 - \frac{2601}{30}\right) = \frac14(90-86.7) = \frac14(3.3)=0.825$. ✓

Answer & the takeaway

$w^\star = 1.7$, MSE $=0.825$, RMSE $\approx 0.908$.

Takeaway: notice the residuals do not sum to zero here ($-0.7-0.4-1.1+1.2 = -1.0$). That is a symptom of forcing $b=0$: without an intercept the line cannot slide up or down to balance itself. Fix that next.

Corner case: when is $\sum_i x^{(i)2} = 0$?

The notes ask this and leave it hanging ("how would the plot of training points look?"). Here is the answer.

The algebra. $S = \sum_i x^{(i)2}$ is a sum of non-negative terms, and such a sum is zero if and only if every term is zero — one positive term could not be cancelled, since none is negative. Since $x^{(i)2}=0 \iff x^{(i)}=0$,

$$\sum_{i=1}^N x^{(i)2} = 0 \quad\Longleftrightarrow\quad x^{(i)} = 0 \ \text{ for every } i = 1,\dots,N.$$

The picture. All $N$ points sit stacked on the vertical axis — a column of dots at $x=0$ at assorted heights, with no horizontal spread at all.

Figure 2 · The degenerate scatter: every $x^{(i)} = 0$
x y 0 all N points sit here, at x = 0 every slope w scores identically
What to notice: all four candidate lines pass through the origin and therefore all predict $\hat y = w\cdot 0 = 0$ at the only input value that appears in the data. The slope is invisible to the loss, so the loss cannot choose one.

Why the slope is genuinely undefined. Plug $x^{(i)}=0$ into the objective:

$$J(w) = \frac1N\sum_{i=1}^N\left(y^{(i)} - 0\cdot w\right)^2 = \frac1N\sum_{i=1}^N y^{(i)2} = \frac{Q}{N}.$$

$J$ does not contain $w$: it is flat, not a parabola. Every $w \in \R$ achieves the same loss $Q/N$, so $\argmin_w J(w) = \R$ — not a point, but the whole real line. Every line through the origin fits these points equally well, or rather equally badly: they all predict $0$ at $x=0$, and the data never probes any other $x$. (They all fit perfectly in the sub-case $y^{(i)}=0$ for all $i$, when all the data sits at $(0,0)$.)

Key idea

The formula "behaves weirdly" ($0/0$, or $P/0$) for a real reason, not a numerical accident: the problem itself has no unique answer. Division by zero in a closed form is nearly always the algebra reporting that the optimisation problem is degenerate. Never patch it; diagnose it.

In one dimension $S$ is exactly $X^\top X$ for the $N\times1$ design matrix built below, so $S=0$ means $X$ is rank deficient. With $d$ features the same disease appears whenever the columns of $X$ are linearly dependent and $X^\top X$ becomes singular — same symptom (infinitely many tied optima), same cure. See the section on $X^\top X$ not being invertible.

Putting the intercept back: the general scalar solution

Real fits need $b$ — a Pontiac Vibe with zero miles is not worth zero dollars. The two-parameter solution reuses everything above.

$$J(w,b) \;=\; \frac1N\sum_{i=1}^N\left(y^{(i)} - w x^{(i)} - b\right)^2.$$
Proof

Take the partial derivative with respect to $b$ first. Using the chain rule, $\frac{\partial}{\partial b}(y_i - wx_i - b)^2 = 2(y_i - wx_i - b)\cdot(-1)$:

$$\frac{\partial J}{\partial b} = -\frac2N\sum_{i=1}^N\left(y_i - wx_i - b\right) = 0 \;\Longrightarrow\; \sum_i y_i - w\sum_i x_i - Nb = 0.$$

Divide by $N$ and write $\bar x = \frac1N\sum_i x_i$, $\bar y = \frac1N\sum_i y_i$ for the sample means:

$$\boxed{\,b = \bar y - w\bar x\,}$$

Notice this holds for whatever $w$ turns out to be. It also says the residuals sum to zero at the optimum — the equation above is literally "$\sum_i \text{residual}_i = 0$". That's the intercept doing its job: sliding the line vertically until the errors balance.

Now substitute that $b$ back in. The $i$-th residual becomes

$$y_i - wx_i - (\bar y - w\bar x) \;=\; (y_i - \bar y) - w(x_i - \bar x).$$

So

$$J(w) = \frac1N\sum_{i=1}^N\Big((y_i-\bar y) - w(x_i - \bar x)\Big)^2,$$

which is exactly the $b=0$ problem we already solved, run on the centred data $\tilde x_i = x_i - \bar x$, $\tilde y_i = y_i - \bar y$. Apply $w^\star = \sum \tilde y_i\tilde x_i / \sum \tilde x_i^2$:

$$\boxed{\;w^\star = \frac{\sum_{i=1}^N (x_i - \bar x)(y_i - \bar y)}{\sum_{i=1}^N (x_i-\bar x)^2},\qquad b^\star = \bar y - w^\star \bar x.\;}\qquad\square$$
Intuition

Adding an intercept is the same as moving the origin to the centre of mass of the data and fitting a through-the-origin line there — so the $b=0$ case was the general case in disguise. The numerator is $N$ times the sample covariance of $x$ and $y$ and the denominator $N$ times the sample variance of $x$, so $w^\star = \widehat{\mathrm{Cov}}(x,y)/\widehat{\mathrm{Var}}(x)$: how much they move together, per unit of how much $x$ moves.

The line always passes through $(\bar x, \bar y)$

$$h(\bar x) = w^\star \bar x + b^\star = w^\star\bar x + (\bar y - w^\star\bar x) = \bar y.$$

The least-squares line is pinned to the centroid of the data and pivots about it as the slope changes — so if you know $\bar x,\bar y$ and the slope you know the whole line, and you can sanity-check any fit by asking whether it passes through the average point.

The matching corner case: $\sum_i(x_i-\bar x)^2 = 0$ exactly when every $x_i$ equals $\bar x$ — all points stacked in a vertical column, now anywhere rather than at the axis. Same disease, same diagnosis: the data never varies $x$, so it can say nothing about the slope.

Example

The eight cars, with mileage in units of 10,000 miles and price in thousands of dollars, give $\bar x = 8.275$, $\bar y = 15.25$, and

$$w^\star = -1.8906\ \text{(thousand dollars per 10{,}000 miles)}, \qquad b^\star = 30.895\ \text{(thousand dollars)}.$$

In raw units: each extra mile costs about 18.9 cents, and a zero-mileage Vibe is priced at about $\$30{,}895$. Best achievable RMSE: $\$3{,}156$ — beating the pre-reading's eyeballed "good" line, which managed $\$3{,}620$. Predicted price at 80,000 miles: $30{,}895 - 0.18906\times80{,}000 \approx \$15{,}770$.

Feel it: drag the line, watch the loss

Two figures, one set of controls: the line and its residuals in data space, then the same $(w,b)$ as a dot on the loss surface $J(w,b)$ that gradient descent will roll down.

Figure 3 · Interactive line fitter — eight cars, mileage in 10k miles, price in thousands of dollars
What to notice: you can make MSE small in many mediocre ways, but only one $(w,b)$ makes it smallest. Try dragging $w$ while watching the residuals flip from "all too high on the left" to "all too high on the right" — that sign pattern is the derivative telling you which way to move.
Figure 4 · The same $(w,b)$ on the loss surface $J(w,b)$
Driven by the sliders above. Dark bands are low loss. What to notice: the valley is a long, tilted ellipse, not a circular bowl — $w$ and $b$ are coupled, because raising the intercept can be partly undone by lowering the slope. The cross marks $(w^\star,b^\star)=(-1.908,\,31.042)$, the exact answer from the closed form. That elongated, tilted valley is precisely the geometry that makes plain gradient descent zig-zag, which is the story of the next chapter.
On the loss surface above, you sit at some $(w,b)$ and increase $b$ by a little. The loss goes up. What can you conclude?
Why: a partial derivative tells you about one direction only. $\partial J/\partial b > 0$ says the slice of the surface at your current $w$ is rising in $b$. Because the valley is tilted, the joint optimum can still have a larger $b$ than yours — you'd just have to change $w$ at the same time. This coupling is the reason we solve for $w$ and $b$ together (or run gradient descent on both coordinates at once) rather than one at a time.

Many features: the design matrix

Real inputs have several features — mileage, age, engine size, previous owners. So let $\X = \R^d$, keep $\Y = \R$, and let a training input be a vector $x^{(i)} \in \R^d$ whose $j$-th coordinate $x^{(i)}_j$ is its $j$-th feature. The hypothesis class becomes all hyperplanes: $h(x) = x^\top w + b$ with $w \in \R^d$, $b\in\R$, where $x^\top w = \sum_{j=1}^d x_jw_j$ is a weighted sum of features and each $w_j$ says "how many dollars per unit of feature $j$".

The bias trick: absorbing $b$

Carrying $b$ separately is a nuisance in the algebra, so fix the first coordinate of every $x$ to be $1$ and stack $b$ on top of $w$:

$$x = \begin{bmatrix}1\\ x_1 \\ \vdots \\ x_{d-1}\end{bmatrix}, \qquad \theta = \begin{bmatrix} b \\ w \end{bmatrix}, \qquad\text{so}\qquad x^\top\theta = 1\cdot b + \sum_j x_j w_j = x^\top w + b.$$

Now $h_\theta(x) = x^\top\theta$ with no separate bias term — the intercept is the weight on a constant feature. From here on $d$ counts that constant feature and $\theta \in \R^d$ is the single parameter vector we solve for.

Definition

The design matrix $X \in \R^{N\times d}$ stacks the training inputs as rows, and the label vector $y \in \R^N$ stacks the labels:

$$X = \begin{bmatrix} \leftarrow\ x^{(1)\top}\ \rightarrow \\ \leftarrow\ x^{(2)\top}\ \rightarrow \\ \vdots \\ \leftarrow\ x^{(N)\top}\ \rightarrow\end{bmatrix} \in \R^{N\times d}, \qquad y = \begin{bmatrix}y^{(1)}\\ y^{(2)} \\ \vdots \\ y^{(N)}\end{bmatrix} \in \R^{N}.$$

Rows are samples, columns are features. Row $i$ is one car; column $j$ is one feature measured across all cars. The first column is all ones.

The single most useful consequence: $X\theta \in \R^N$ is the vector of all $N$ predictions at once, because row $i$ of $X\theta$ is $x^{(i)\top}\theta = h_\theta(x^{(i)})$. So $y - X\theta$ is the vector of all $N$ residuals, and

$$\|y - X\theta\|_2^2 = \sum_{i=1}^N\left(y^{(i)} - x^{(i)\top}\theta\right)^2$$

(recall $\|v\|_2^2 = \sum_k v_k^2$, the squared Euclidean length). Therefore ERM in matrix form is simply

$$\theta^\star = \argmin_{\theta \in \R^d}\ \frac1N\|y - X\theta\|_2^2.$$
Intuition

Nothing was added, only repackaged: $N$ separate squared errors became one squared length of one residual vector in $\R^N$. Least squares is: make the residual vector short.

Answering the notes: "can you simply set $\theta = X^{-1}y$?"

The notes pose this parenthetically and answer "No." There are four distinct reasons.

1. $X$ is usually not square, so $X^{-1}$ is not a hard object to compute — it is not an object at all. Inverses are defined only for square matrices, and $X$ is $N\times d$ with $N$ normally much larger than $d$.

2. Even when $N=d$, $X$ is usually still singular. A square matrix is invertible only if its rows (equivalently columns) are linearly independent; two identical cars give two identical rows, and two features that are the same measurement in different units give two dependent columns.

3. The system is overdetermined. "$\theta = X^{-1}y$" asks for $X\theta = y$ exactly — a line through every data point, i.e. $N=8$ constraints on $d=2$ unknowns. Generic data admits no solution: look at Figure 1, where no straight line touches all eight dots.

Key idea

As $\theta$ ranges over $\R^d$, the prediction vector $X\theta$ traces out exactly the column space of $X$ — a $d$-dimensional subspace of $\R^N$. When $d < N$ that is a thin slice, and thanks to noise $y$ almost surely does not lie in it. So "solve $X\theta = y$" is the wrong question; the right one is which point of the column space is closest to $y$? The answer is the orthogonal projection of $y$ onto $\mathrm{col}(X)$.

Figure 5 · Least squares is a projection onto the column space
col(X) = { Xθ : θ ∈ ℝᵈ }, a d-dimensional slice of ℝᴺ 0 y (the labels) not in the plane — noise lifts it off ŷ = Xθ* the closest point in the plane residual y − Xθ* ⊥ to every column of X some other Xθ longer
What to notice: the shortest segment from a point to a plane is the perpendicular one. Writing "perpendicular" algebraically — every column of $X$ is orthogonal to the residual, i.e. $X^\top(y - X\theta^\star) = 0$ — rearranges instantly to $X^\top X\theta^\star = X^\top y$. That is where the normal equations get their name: normal means perpendicular.

4. And even if you could invert $X$, you shouldn't want to. With $N=d$ and $X$ invertible, $\theta = X^{-1}y$ would pass exactly through every training point, measurement errors included: training loss $0$ and a model that has memorised the noise. That least squares refuses to interpolate when $N > d$ is a feature.

The normal equations, derived in full

Drop the $\frac1N$ (positive constant) and expand the squared norm using $\|v\|_2^2 = v^\top v$:

$$\argmin_\theta\ (y - X\theta)^\top(y - X\theta).$$

Step 1 — multiply it out

Using $(a-b)^\top(c-d) = a^\top c - a^\top d - b^\top c + b^\top d$ and $(X\theta)^\top = \theta^\top X^\top$:

$$(y-X\theta)^\top(y-X\theta) \;=\; y^\top y \;-\; \theta^\top X^\top y \;-\; y^\top X\theta \;+\; \theta^\top X^\top X\theta.$$

Step 2 — the two middle terms are the same number

Proof

Check the shapes: $\theta^\top$ is $1\times d$, $X^\top$ is $d\times N$, $y$ is $N\times1$, so $\theta^\top X^\top y$ is $1\times1$ — a scalar. Likewise $y^\top X\theta$. A $1\times1$ matrix equals its own transpose, and transposing the first term gives the second:

$$\left(\theta^\top X^\top y\right)^\top = y^\top \left(X^\top\right)^\top \left(\theta^\top\right)^\top = y^\top X\theta.$$

So $\theta^\top X^\top y = y^\top X\theta$, and the two terms combine into $-2\theta^\top X^\top y$. $\square$

Also $y^\top y$ contains no $\theta$ — drop it by the same argmin rule. We are left with the notes' expression:

$$\argmin_\theta\ -2\theta^\top X^\top y + \theta^\top X^\top X\theta.$$

Step 3 — name the pieces

Definition $$A = X^\top X \in \R^{d\times d}, \qquad b = X^\top y \in \R^{d}.$$

(Careful: this $b$ is the notes' reuse of the letter — a vector, not the scalar intercept, which now lives inside $\theta$.)

$A$ is the matrix version of $S = \sum_i x^{(i)2}$ — its $(j,k)$ entry is $\sum_{i=1}^N x^{(i)}_j x^{(i)}_k$ — and $b$ is the matrix version of $P = \sum_i y^{(i)}x^{(i)}$, its $j$-th entry being $\sum_{i=1}^N x^{(i)}_j y^{(i)}$. The answer will be $A^{-1}b$: "$b$ divided by $A$".

Two facts about $A$ we will need. $A$ is symmetric: $A^\top = (X^\top X)^\top = X^\top (X^\top)^\top = X^\top X = A$. $A$ is positive semidefinite: for any $v\in\R^d$, $v^\top A v = (Xv)^\top(Xv) = \|Xv\|_2^2 \ge 0$. If $A$ is additionally invertible then $\|Xv\|^2 > 0$ for all $v\neq0$, i.e. $A$ is positive definite. Assume that for now.

Step 4 — complete the square, in matrices

Proof

Claim. For symmetric invertible $A$,

$$-2\theta^\top b + \theta^\top A\theta \;=\; (\theta - A^{-1}b)^\top A (\theta - A^{-1}b) \;-\; b^\top A^{-1}b.$$

Verify by expanding the right-hand side. Write $m = A^{-1}b$. Then

$$(\theta-m)^\top A(\theta-m) = \theta^\top A\theta - \theta^\top A m - m^\top A\theta + m^\top A m.$$

Handle the four pieces:

Substituting:

$$(\theta-m)^\top A(\theta-m) = \theta^\top A\theta - \theta^\top b - \theta^\top b + b^\top A^{-1}b = \theta^\top A\theta - 2\theta^\top b + b^\top A^{-1}b.$$

Subtract $b^\top A^{-1}b$ from both sides and the claim follows. $\square$

Key idea

Symmetry of $A$ is doing real work, not cosmetics. Without $A^\top = A$ the two cross terms $\theta^\top Am$ and $m^\top A\theta$ would not collapse into $2\theta^\top b$, and there would be no clean completed square. This is exactly why the derivation is stated for $A = X^\top X$: that product is symmetric no matter what $X$ is.

Step 5 — read off the minimiser

Substituting the identity and dropping the constant $-b^\top A^{-1}b$:

$$\argmin_\theta\ (\theta - A^{-1}b)^\top A(\theta - A^{-1}b).$$

Since $A = X^\top X$, that quadratic form is

$$(\theta - A^{-1}b)^\top X^\top X(\theta - A^{-1}b) = \big\|X\left(\theta - (X^\top X)^{-1}X^\top y\right)\big\|_2^2,$$

the notes' final line. A squared norm is $\ge 0$ and hits $0$ when its argument does, at $\theta = A^{-1}b$ — achievable, so it is the minimum, and with $A$ positive definite it is the unique one.

Key idea $$\boxed{\;\theta^\star \;=\; A^{-1}b \;=\; \left(X^\top X\right)^{-1}X^\top y\;}$$

The normal equations. Compare the scalar answer $w^\star = P/S$: same formula, with $S \to X^\top X$, $P \to X^\top y$, and "divide" $\to$ "multiply by the inverse". In practice nobody forms the inverse; you solve the linear system $(X^\top X)\theta = X^\top y$, which is faster and numerically far better behaved.

Step 6 — the calculus route, for cross-check

The optimization chapter proves $\nabla_\theta(\theta^\top v) = v$ and $\nabla_\theta(\theta^\top A\theta) = (A + A^\top)\theta$. Applied to $J(\theta) = y^\top y - 2\theta^\top X^\top y + \theta^\top X^\top X\theta$, using symmetry of $X^\top X$ at the last step:

$$\nabla_\theta J(\theta) \;=\; 0 \;-\; 2X^\top y \;+\; (X^\top X + (X^\top X)^\top)\theta \;=\; -2X^\top y + 2X^\top X\theta,$$ $$\nabla_\theta J = 0 \;\Longrightarrow\; X^\top X\theta = X^\top y \;\Longrightarrow\; \theta^\star = (X^\top X)^{-1}X^\top y.$$

Same answer, third time. And $X^\top X\theta = X^\top y$ rearranges to $X^\top(y - X\theta) = 0$ — the geometric statement from Figure 5 that the residual is orthogonal to every column of $X$. Algebra, calculus and geometry land on the same equation.

Guided problemNormal equations by hand, with a 2×2 inverse

Three points: $(x,y) = (0,1),\ (1,2),\ (2,5)$. Fit $h(x) = wx + b$ using $\theta^\star = (X^\top X)^{-1}X^\top y$ with the bias absorbed. Do every arithmetic step.

Reminder: $\begin{bmatrix}p&q\\r&s\end{bmatrix}^{-1} = \dfrac{1}{ps-qr}\begin{bmatrix}s&-q\\-r&p\end{bmatrix}$.

Step 1 — Build $X$ and $y$

Absorbing the bias means each input becomes $\begin{bmatrix}1 \\ x\end{bmatrix}$ and $\theta = \begin{bmatrix}b\\w\end{bmatrix}$. Rows of $X$ are the inputs:

$$X = \begin{bmatrix}1 & 0\\ 1 & 1\\ 1 & 2\end{bmatrix} \in \R^{3\times2}, \qquad y = \begin{bmatrix}1\\2\\5\end{bmatrix}\in\R^3.$$

Sanity check: $N=3$ rows, $d=2$ columns. The first column is all ones — that's the constant feature.

Step 2 — Compute $A = X^\top X$

$X^\top = \begin{bmatrix}1&1&1\\0&1&2\end{bmatrix}$ is $2\times3$, so $A$ is $2\times2$. Entry $(j,k)$ is the dot product of column $j$ and column $k$ of $X$:

  • $A_{11} = 1^2+1^2+1^2 = 3$ (this is always $N$ for the ones-column)
  • $A_{12} = A_{21} = 1\cdot0 + 1\cdot1 + 1\cdot2 = 3$ (this is $\sum_i x_i$)
  • $A_{22} = 0^2+1^2+2^2 = 5$ (this is $\sum_i x_i^2$)
$$A = \begin{bmatrix}3 & 3\\ 3 & 5\end{bmatrix}.$$

Symmetric, as promised.

Step 3 — Compute $b = X^\top y$
  • First entry: $1\cdot1 + 1\cdot2 + 1\cdot5 = 8$ (this is $\sum_i y_i$)
  • Second entry: $0\cdot1 + 1\cdot2 + 2\cdot5 = 12$ (this is $\sum_i x_iy_i$)
$$b = \begin{bmatrix}8\\12\end{bmatrix}.$$
Step 4 — Invert $A$

$\det A = 3\cdot5 - 3\cdot3 = 15 - 9 = 6 \neq 0$, so the inverse exists.

$$A^{-1} = \frac16\begin{bmatrix}5 & -3\\ -3 & 3\end{bmatrix}.$$

Check: $A^{-1}A = \frac16\begin{bmatrix}5\cdot3-3\cdot3 & 5\cdot3-3\cdot5\\ -3\cdot3+3\cdot3 & -3\cdot3+3\cdot5\end{bmatrix} = \frac16\begin{bmatrix}6&0\\0&6\end{bmatrix} = I$. ✓

Step 5 — Multiply
$$\theta^\star = A^{-1}b = \frac16\begin{bmatrix}5&-3\\-3&3\end{bmatrix}\begin{bmatrix}8\\12\end{bmatrix} = \frac16\begin{bmatrix}40-36\\ -24+36\end{bmatrix} = \frac16\begin{bmatrix}4\\12\end{bmatrix} = \begin{bmatrix}2/3\\ 2\end{bmatrix}.$$

So $b^\star = 2/3$ and $w^\star = 2$: the fitted line is $h(x) = 2x + \tfrac23$.

Answer & the takeaway

$h(x) = 2x + \tfrac23$. Cross-check with the scalar with-intercept formulas: $\bar x = 1$, $\bar y = 8/3$.

$$w^\star = \frac{(-1)(1-\tfrac83) + 0 + (1)(5-\tfrac83)}{(-1)^2+0^2+1^2} = \frac{\tfrac53 + \tfrac73}{2} = \frac{4}{2} = 2,\qquad b^\star = \tfrac83 - 2\cdot1 = \tfrac23.\ \checkmark$$

Residuals: $1 - \tfrac23 = \tfrac13$; $2 - \tfrac83 = -\tfrac23$; $5-\tfrac{14}{3} = \tfrac13$. They sum to $\tfrac13-\tfrac23+\tfrac13 = 0$, exactly as the $\partial J/\partial b = 0$ condition predicted.

Takeaway: every entry of $A$ and $b$ was a familiar scalar sum in disguise. The matrix machinery is bookkeeping for sums you already know how to compute — its value is that it keeps working when $d = 500$.

Corner case: when $X^\top X$ is not invertible

The notes flag this and say "you don't have to worry about it". You should at least be able to recognise it.

When. $\det(X^\top X) = 0$ exactly when $\mathrm{rank}(X) < d$, i.e. the columns of $X$ are linearly dependent. Two routes there:

Example

Record each car's weight twice, in thousands of kilograms and thousands of pounds, so column 2 is $2.2\times$ column 1:

$$X = \begin{bmatrix}1 & 2.2\\ 2 & 4.4\\ 3 & 6.6\end{bmatrix},\qquad X^\top X = \begin{bmatrix}14 & 30.8\\ 30.8 & 67.76\end{bmatrix}, \qquad \det = 14\times67.76 - 30.8^2 = 0.$$

Singular — and the fit is visibly ambiguous: if $\theta = (\theta_1,\theta_2)$ fits, so does $(\theta_1 + 2.2c,\ \theta_2 - c)$ for any $c$. "Kilograms matter and pounds don't" and the reverse are the same model.

What it means. If $v \neq 0$ satisfies $Xv = 0$ (a null-space direction, which exists precisely when the columns are dependent), then for any optimal $\theta^\star$,

$$X(\theta^\star + cv) = X\theta^\star + cXv = X\theta^\star \quad\text{for every } c \in \R.$$

Identical predictions, identical loss: infinitely many optima, forming a whole line (or plane) of tied solutions. Crucially the projection $\hat y = X\theta^\star$ is still unique — Figure 5 still has one closest point in the plane. What is not unique is the coordinates $\theta$ naming that point.

What the pseudo-inverse picks. The Moore–Penrose pseudo-inverse $X^+$ is defined for every matrix, and $\theta = X^+y$ always returns a least-squares solution — specifically the tied one with the smallest $\|\theta\|_2$. When $X^\top X$ is invertible, $X^+ = (X^\top X)^{-1}X^\top$, so nothing is lost by always using it (numpy.linalg.lstsq does; numpy.linalg.inv blows up). beyond syllabus The other standard fix is ridge regression: minimise $\|y-X\theta\|^2 + \lambda\|\theta\|^2$, giving $(X^\top X + \lambda I)^{-1}X^\top y$, invertible for any data because $\lambda I$ shifts every eigenvalue up by $\lambda$.

You have $N=50$ documents and $d=5{,}000$ word-count features. What happens if you try $\theta = (X^\top X)^{-1}X^\top y$?
Why: $X^\top X$ is always square ($d\times d$) — squareness is never the problem. Rank is. $\mathrm{rank}(X^\top X) = \mathrm{rank}(X) \le \min(N,d) = 50$, far below $5{,}000$, so the determinant is $0$. With $d \gg N$ you can fit the 50 labels perfectly in infinitely many ways; the training loss is $0$ for a whole $4{,}950$-dimensional family of $\theta$. You need the pseudo-inverse, or regularisation, or fewer features — and you should be very suspicious of the resulting test error.

What it costs, and why the next chapter exists

The closed form is wonderful, and for large $d$ unusable. Count the operations:

StepWhat it doesCost
$X^\top X$$d^2$ entries, each a dot product of two length-$N$ columns$O(Nd^2)$
$X^\top y$$d$ entries, each a length-$N$ dot product$O(Nd)$
Solve $(X^\top X)\theta = X^\top y$Cholesky factorisation of a $d\times d$ system$O(d^3)$
Totalplus $O(d^2)$ memory just to hold $X^\top X$$\mathbf{O(Nd^2 + d^3)}$

At $d = 10$, $N = 10^6$ that is $10^8$ operations — instant, always use the closed form. At $d = 10^5$ (a modest last layer, or bag-of-words) it is $d^3 = 10^{15}$ operations and $\approx$ 80 GB just for $X^\top X$: hopeless, and the cubic term does not care how much data you have.

Key idea

Gradient descent computes $\nabla_\theta J = -2X^\top y + 2X^\top X\theta$, but it never forms $X^\top X$. It evaluates $X^\top(X\theta)$ as two matrix–vector products: $O(Nd)$ per iteration, $O(d)$ memory. Run it for $T$ iterations and you pay $O(TNd)$ — linear in $d$, not cubic. For $d = 10^5$ that is the difference between "runs in a minute" and "does not run".

And linear regression is the easy case: it is one of the very few models with any closed form at all. Logistic regression has none. Neural networks have none. Gradient descent is the general-purpose answer, and it is the subject of the next chapter.

Guided problemA fit from summary statistics — should I sell the Vibe?

Someone hands you only these numbers for the eight cars (mileage in miles, price in dollars): $\bar x = 82{,}750$, $\bar y = 15{,}250$, $\sum_i (x_i-\bar x)(y_i-\bar y) = -2.1268\times10^9$, $\sum_i (x_i-\bar x)^2 = 1.1249\times10^{10}$.

(a) Find the fitted line. (b) Predict the price of a Vibe with 80,000 miles. (c) The car actually has 140,000 miles; what does the model say, and what should you notice?

Step 1 — Which formula applies?

We want an intercept (a zero-mileage car is not free), so use the centred formulas, not $P/S$:

$$w^\star = \frac{\sum_i(x_i-\bar x)(y_i-\bar y)}{\sum_i(x_i-\bar x)^2},\qquad b^\star = \bar y - w^\star\bar x.$$

Note we never needed the raw data — those four summaries are sufficient.

Step 2 — Slope and intercept
$$w^\star = \frac{-2.1268\times10^9}{1.1249\times10^{10}} = -0.18906\ \text{dollars per mile}.$$ $$b^\star = 15{,}250 - (-0.18906)(82{,}750) = 15{,}250 + 15{,}645 = 30{,}895\ \text{dollars}.$$

So $h(x) = 30{,}895 - 0.18906\,x$. Sign check: the slope is negative, which is what "more miles, less money" should look like. If you got a positive slope, you swapped the sign somewhere.

Step 3 — Predict at 80,000 miles
$$h(80{,}000) = 30{,}895 - 0.18906\times80{,}000 = 30{,}895 - 15{,}125 = 15{,}770.$$

About $\$15{,}770$. Sanity check with the centroid property: $80{,}000$ is just below $\bar x = 82{,}750$, and indeed $15{,}770$ is just above $\bar y = 15{,}250$. The line passes through $(\bar x,\bar y)$, so predictions near the mean mileage should land near the mean price.

Step 4 — Predict at 140,000 miles, and be suspicious
$$h(140{,}000) = 30{,}895 - 0.18906\times140{,}000 = 30{,}895 - 26{,}468 = 4{,}427.$$

About $\$4{,}427$. Two things to notice.

(i) Extrapolation. $140{,}000$ is inside the training range ($25{,}000$–$160{,}000$), so this is interpolation and reasonably trustworthy. But the same formula at $x = 200{,}000$ gives $30{,}895 - 37{,}812 = -\$6{,}917$: a negative price. The model has no idea that prices are bounded below by zero. A line is linear everywhere, including where reality is not.

(ii) RMSE is $\$3{,}156$. Against a prediction of $\$4{,}427$ that is a huge relative uncertainty. A point estimate without its error bar is a half-answer.

Answer & the takeaway

(a) $h(x) = 30{,}895 - 0.18906x$. (b) $\approx\$15{,}770$. (c) $\approx\$4{,}427$, with RMSE $\approx\$3{,}156$, and the model goes negative beyond about $163{,}400$ miles.

Takeaway: the fit depends on the data only through four summary numbers, the fitted line always passes through $(\bar x,\bar y)$, and the model's confident-looking arithmetic says nothing about whether a line was the right hypothesis class. Checking that is what train/test splits are for.

Guided problemDiagnose three broken designs

For each, say whether the closed-form solution is unique, and why.

(A) $b=0$ model, data $(0,4),(0,-1),(0,7)$.    (B) with-intercept model, data $(5,3),(5,9),(5,1)$.    (C) $d=3$ features (constant, mileage, mileage-in-kilometres), $N=100$ cars.

(A) — do the arithmetic before reasoning

$S=\sum x_i^2 = 0+0+0 = 0$ and $P = \sum x_iy_i = 0$. The formula reads $0/0$: undefined. And indeed $J(w) = \frac13(16+1+49) = 22$ for every $w$ — flat. Infinitely many optima; the slope is unidentifiable because the data never varies $x$. This is exactly the notes' corner case.

(B) — same disease, shifted

$\bar x = 5$, so every $(x_i - \bar x) = 0$ and the denominator $\sum(x_i-\bar x)^2 = 0$. Undefined again. Picture: three points stacked vertically at $x=5$. Any line through $(5,\bar y) = (5,\tfrac{13}{3})$ fits identically, whatever its slope. The intercept did not save us: the problem is zero horizontal spread, not the position of the origin.

(C) — rank deficiency with plenty of data

Column 3 is $1.609\times$ column 2, so $\mathrm{rank}(X) = 2 < 3 = d$ even though $N=100 \gg d$. $X^\top X$ is $3\times3$ with determinant $0$. Not unique: shifting weight between the miles column and the kilometres column leaves every prediction untouched. Lots of data does not fix redundant features. Drop one column, or use the pseudo-inverse.

Answer & the takeaway

None of the three has a unique solution. (A) and (B) are the scalar corner case — zero variation in $x$. (C) is its multivariate form — linearly dependent columns.

Takeaway: "the denominator is zero" and "$X^\top X$ is singular" are the same statement at two different values of $d$. In both cases the fitted predictions may still be perfectly well defined; it is the parameters that are ambiguous. If your job is prediction you can shrug; if your job is to interpret the coefficients ("how much is a mile worth?"), you cannot.

Which of these is not a valid reason that $\theta^\star = (X^\top X)^{-1}X^\top y$ rather than $X^{-1}y$?
Why: $X^\top X$ is not always invertible — that's the whole point of the pseudo-inverse discussion. It is invertible exactly when $X$ has full column rank, which fails when $N < d$ or when features are redundant. The other two statements are correct and are the real reasons.

Recap

One feature: $w^\star=\frac{\sum_i y_ix_i}{\sum_i x_i^2}$ through the origin, or the centred formulas with an intercept, and the line passes through $(\bar x,\bar y)$ — breaking when $x$ has no spread. Many features: absorb $b$ as a constant-$1$ feature, stack samples as rows of $X$, and $\theta^\star=(X^\top X)^{-1}X^\top y$ is the orthogonal projection of $y$ onto $\mathrm{col}(X)$ — breaking when $\mathrm{rank}(X) < d$, where the pseudo-inverse takes over.