Linear Models in Matrix Form

<p class="kicker">MST0052 · Lecture 4 supplement · Optional</p>

## Linear Models<br>in Matrix Form

### Vegard H. Larsen — BI Norwegian Business School

---

## What this deck is

- **Optional** — nothing here is required for the project
- The **matrix algebra** behind Lecture 4's formulas
- Depth for the **oral exam** — *why* ridge works, not just when
- New required material

</div>

---

## The model, in matrices

---

## All n equations at once

L4 wrote one prediction. Stack all **n observations**:

$$y = X\beta + \varepsilon$$

$$\underbrace{y}\_{n \times 1} \qquad \underbrace{X}\_{n \times (p+1)} \qquad \underbrace{\beta}\_{(p+1) \times 1} \qquad \underbrace{\varepsilon}\_{n \times 1}$$

---

## Anatomy of the design matrix

$$X = \begin{bmatrix} 1 & x\_{11} & \cdots & x\_{1p} \\\\ 1 & x\_{21} & \cdots & x\_{2p} \\\\ \vdots & \vdots & & \vdots \\\\ 1 & x\_{n1} & \cdots & x\_{np} \end{bmatrix}$$

- Each **row** — one observation
- Each **column** — one feature
- The column of **ones** carries the intercept

---

## The objective, compactly

The sum of squared residuals from L4 is a **squared length**:

$$\hat{\beta} = \arg\min\_{\beta} \lVert y - X\beta \rVert^2$$

Find the β that makes the **error vector shortest**.

---

## The solution

---

## Set the gradient to zero

$$\nabla\_{\beta} \lVert y - X\beta \rVert^2 = -2X^\top(y - X\beta) = 0$$

which gives the **normal equations**:

$$X^\top X \hat{\beta} = X^\top y$$

---

## The closed-form solution

$$\hat{\beta} = (X^\top X)^{-1} X^\top y$$

- One formula, **no iteration** — this is `LinearRegression().fit()`
- Valid only when X<sup>⊤</sup>X is **invertible**

<p class="muted">scikit-learn actually solves this via SVD — same answer, better numerics (backup).</p>

---

## The geometry

![y projected onto the column space of X](../figures/ols-projection.svg)

ŷ = Xβ̂ is the **projection** of y onto the columns of X — the residual is **orthogonal**.

---

## When the formula breaks

---

## When can't you invert X<sup>⊤</sup>X?

- A feature is an exact **combination** of others — the dummy trap
- **More features than observations** — p > n
- Then many different β fit the training data **equally well**

---

## Near-singular is almost as bad

- Highly **correlated features** → X<sup>⊤</sup>X *barely* invertible
- β̂ still exists, but its **variance explodes**
- Tiny changes in the data → wildly different coefficients

<p class="muted">This is the "variance" in the bias-variance tradeoff — Lecture 6.</p>

---

## Ridge, mechanically

$$\hat{\beta}\_{\text{ridge}} = (X^\top X + \lambda I)^{-1} X^\top y$$

- The **+λI** makes the matrix invertible for any λ > 0
- Shrinkage is not a hack — it is the **fix for the broken inverse**

<p class="muted">Closed form of the penalised objective from the L4 deck.</p>

---

## A linear model is a tiny neural network

---

## One neuron, identity activation

![A linear model drawn as a single neuron](../figures/linear-neuron.svg)

---

## Two roads to the same β̂

<div class="cols">
<div class="callout">
<h3>Closed form</h3>
<p>Solve the normal equations in one shot. What <code>LinearRegression</code> does.</p>
</div>
<div class="callout">
<h3>Gradient descent</h3>
<p>Step downhill along −∇ until it flattens. How neural networks train.</p>
</div>
</div>

The loss surface is a **convex bowl** — both land in the same place.

---

## See it in code

```python
from sklearn.linear_model import LinearRegression, SGDRegressor

ols = LinearRegression().fit(X_s, y)
sgd = SGDRegressor(loss="squared_error", penalty=None,
                   max_iter=10_000, tol=1e-12,
                   random_state=42).fit(X_s, y)

print(ols.coef_.round(2))  # [27.48 79.   31.05 71.29 10.13]
print(sgd.coef_.round(2))  # [27.45 79.03 31.05 71.29 10.14]
```

<p class="muted">Pre-tested: <code>scripts/examples/lecture_04_linear_models_supplement.py</code>. Standardise features first — gradient descent assumes comparable scales.</p>

---

## Why this matters for Lecture 14

- Add a **sigmoid** to the neuron → logistic regression (L5)
- Add **hidden layers** → the network learns its own features
- **Weight decay** in neural nets is exactly ridge's λ
- You have been training neural networks **since Lecture 4**

---

## Questions

Backup slides below ↓

## Backup — the derivation in full

$$\lVert y - X\beta \rVert^2 = y^\top y - 2\beta^\top X^\top y + \beta^\top X^\top X \beta$$

$$\nabla\_{\beta} = -2X^\top y + 2X^\top X\beta = 0 \quad\Rightarrow\quad X^\top X\hat{\beta} = X^\top y$$

<p class="muted">The Hessian 2X<sup>⊤</sup>X is positive semi-definite — the bowl is convex, the stationary point is the global minimum.</p>

## Backup — the variance of β̂

$$\mathrm{Var}(\hat{\beta}) = \sigma^2 (X^\top X)^{-1}$$

- Nearly collinear columns → near-zero **eigenvalues** of X<sup>⊤</sup>X
- Inverting divides by them → **huge** coefficient variance

## Backup — what +λI does to the eigenvalues

$$\text{eigenvalues of } X^\top X: d\_j \quad\longrightarrow\quad \text{of } X^\top X + \lambda I: d\_j + \lambda$$

- Every eigenvalue lifted **away from zero**
- Inverse entries bounded by **1/λ** — variance tamed
- Directions with small d<sub>j</sub> get shrunk **hardest**

## Backup — what scikit-learn actually computes

- Nobody inverts X<sup>⊤</sup>X in practice — **numerically fragile**
- `LinearRegression` → SVD-based least squares (`scipy.lstsq`)
- `Ridge` → Cholesky / sparse solvers depending on data
- Same β̂, far better behaviour with near-collinear columns

---

## What's next

Back to the **main track** — this deck returns twice:

- **Lecture 6:** the variance you just saw, as a tradeoff
- **Lecture 14:** the neuron, with depth