Linear Regression OLS: The Normal Equation

Gradient descent finds the optimal weights iteratively. The normal equation finds them in one algebraic step. For small datasets, this is faster and exact — no learning rate to tune, no convergence to wait for. Understanding the derivation also reveals exactly when the normal equation fails and why Ridge regression fixes it.

The Plan — From Calculus to Code

We'll derive the normal equation from scratch (gradient → set to zero → solve), then compute every matrix step by hand on our 6-house dataset, then compare to gradient descent on complexity.

Step 1 — Setting Up the Matrix System

We want to minimize:

$J(\mathbf{w})=\Vert \mathbf{y}-X\mathbf{w}{\Vert }^2=(\mathbf{y}-X\mathbf{w}{)}^{\top }(\mathbf{y}-X\mathbf{w})$

Expand the product:

$J(\mathbf{w})={\mathbf{y}}^{\top }\mathbf{y}-2{\mathbf{w}}^{\top }{X}^{\top }\mathbf{y}+{\mathbf{w}}^{\top }{X}^{\top }X\mathbf{w}$

Take the gradient with respect to $\mathbf{w}$ and set to zero:

$\frac{\partial J}{\partial \mathbf{w}}=-2X^{\top }\mathbf{y}+2X^{\top }X\mathbf{w}=0$

Rearranging gives the normal equations:

$X^{\top }X\mathbf{w}=X^{\top }\mathbf{y}$

If $X^{\top }X$ is invertible, the unique solution is:

${\mathbf{w}}^{\ast }=(X^{\top }X{)}^{-1}{X}^{\top }\mathbf{y}$

The geometric interpretation: ${\mathbf{w}}^{\ast }$ projects $\mathbf{y}$ onto the column space of $X$. The fitted values $X{\mathbf{w}}^{\ast }$ are the closest point in the column space to $\mathbf{y}$.

✓ Step 1 complete. The normal equations $X^{\top }X\mathbf{w}=X^{\top }\mathbf{y}$ come from setting the gradient of $J(\mathbf{w})$ to zero. The solution projects $\mathbf{y}$ onto the column space of $X$.

Anchor dataset (2-feature, design matrix with intercept column):

import numpy as np

X_aug = np.array([
    [1, 650,  2],
    [1, 850,  2],
    [1, 1100, 3],
    [1, 1400, 3],
    [1, 1600, 4],
    [1, 1900, 4]
], dtype=float)

y = np.array([180, 220, 280, 340, 370, 430], dtype=float)

Computing $X^{\top }X$

$X^{\top }$ is $3\times 6$, $X$ is $6\times 3$, so $X^{\top }X$ is $3\times 3$:

XtX = X_aug.T @ X_aug
print(XtX)

[[    6.    7500.   18.  ]
 [ 7500. 10985000. 23400.]
 [   18.   23400.   62.  ]]

Diagonal: $[n=6,\sum x_{i1}^2=10985000,\sum x_{i2}^2=62]$

Off-diagonal: $[\sum x_{i1}=7500,\sum x_{i2}=18,\sum x_{i1}{x}_{i2}=23400]$

✓ Step 2 complete. $X^{\top }X$ is $3\times 3$: diagonal = squared sums, off-diagonal = cross terms.

Step 3 — Computing $X^{\top }\mathbf{y}$

$X^{\top }$ is $3\times 6$, $\mathbf{y}$ is $6\times 1$, so $X^{\top }\mathbf{y}$ is $3\times 1$:

Xty = X_aug.T @ y
print(Xty)

[  1820.    2497000.    5510.   ]

These are $[\sum y_i=1820,\sum x_{i1}{y}_i=2497000,\sum x_{i2}{y}_i=5510]$.

✓ Step 3 complete. $X^{\top }\mathbf{y}$ is the $3\times 1$ vector of target sums, cross-terms with each feature.

Step 4 — Solving the Normal Equations

Two approaches: direct inversion and LU decomposition. Always prefer the latter.

# Method 1: direct inversion (numerically unstable for ill-conditioned matrices)
w_inv = np.linalg.inv(XtX) @ Xty

# Method 2: np.linalg.solve (uses LU decomposition — preferred)
w_ols = np.linalg.solve(XtX, Xty)

print(f"w₀ = {w_ols[0]:.4f}")
print(f"w₁ = {w_ols[1]:.6f}")
print(f"w₂ = {w_ols[2]:.4f}")

w₀ = 58.3271
w₁ = 0.175000
w₂ = 9.6053

np.linalg.solve(XtX, Xty) uses LU decomposition — avoids computing the inverse directly. Always prefer solve over inv unless you explicitly need the inverse matrix (almost never).

✓ Step 4 complete. OLS weights: $w_0=58.33$, $w_1=0.175$, $w_2=9.61$. SSE = 743.2 — higher than simple regression due to near-multicollinearity.

Step 5 — Verify: Prediction and Residuals

y_pred = X_aug @ w_ols
sse = ((y - y_pred) ** 2).sum()
print(f"SSE (OLS): {sse:.4f}")

SSE (OLS): 743.2191

This is larger than simple regression's SSE of 133.3 — not because multiple regression is worse, but because the near-multicollinearity ($r=0.997$ between sqft and bedrooms) destabilizes the coefficients. The OLS solution is still mathematically optimal for these features, but the features themselves add little independent signal.

✓ Step 5 complete. Predictions with $\hat{y}=X_{\text{aug}}\mathbf{w}$ give SSE = 743.2. The residuals show larger errors in the tails — a symptom of multicollinearity.

Step 6 — Why $(X^{\top }X{)}^{-1}$ Exists — and When It Doesn't

$(X^{\top }X{)}^{-1}$ exists if and only if $X$ has full column rank — no perfect linear dependence among the feature columns.

Three situations where $X^{\top }X$ is singular:

1. More features than samples ($p>n$): $X^{\top }X$ is $p\times p$ but rank $\le n<p$ — underdetermined system.
2. Perfect multicollinearity: feature A = $2\times$ feature B → the two columns are linearly dependent → singular.
3. Zero-variance column: a feature with identical values for all samples has zero variance → column of constants (besides the bias) → dependent with the bias column.

Fix when singular: Ridge regression adds $\lambda I$ to $X^{\top }X$:

${\mathbf{w}}^{\ast }=(X^{\top }X+\lambda I{)}^{-1}{X}^{\top }\mathbf{y}$

Adding $\lambda I$ to the diagonal guarantees positive-definiteness and invertibility for any $\lambda >0$, regardless of multicollinearity.

Computational Complexity

For $n=1,000,000$ samples and $p=1,000$ features: $p^3=10^9$ operations — hours of computation. Gradient descent costs $O(np)$ per iteration — far better for large $p$.

Normal Equation vs Gradient Descent

✓ Step 6 complete. Normal equation is closed-form and exact but costs $O(n{p}^2+p^3)$. Gradient descent is iterative but scales to large $n$ and $p$.

Related Concepts and Honest Limitations

The normal equation is the analytical solution to the OLS problem — it's not an approximation. For pure linear regression on small datasets ($n<100k$, $p<1000$), it's the right choice.

A caveat that practitioners discover the hard way: the normal equation gives you the optimal weights for MSE, but "optimal" doesn't mean "useful." When features are nearly collinear, the coefficients minimize train error but are statistically unstable — small data changes cause large coefficient swings. Ridge regularization addresses this by changing the objective (adding $\lambda \Vert \mathbf{w}{\Vert }^2$), not by finding a different solution to the original objective.

Another blind spot: the normal equation assumes $X^{\top }X$ is invertible. If $X$ has duplicate columns or $p>n$, the solution literally doesn't exist without Ridge. Gradient descent, by contrast, will converge to some solution (the minimum-norm least-squares solution, depending on initialization), but the coefficients still aren't uniquely determined.

Test Your Understanding

1. Compute $X^{\top }X$ by hand for the 1-feature anchor ($n=6$, design matrix with intercept column). Verify the diagonal entries are $[n,\sum x_i^2]$.

2. The np.linalg.solve function uses LU decomposition. Why does this avoid the numerical instability of inv()? What does "ill-conditioned" mean for $X^{\top }X$?

3. You have $n=500$ samples and $p=600$ features. What does $X^{\top }X$ look like, and why is $(X^{\top }X{)}^{-1}$ undefined without Ridge?

4. Adding $\lambda I$ to $X^{\top }X$ makes it invertible. What happens to the OLS solution as $\lambda \to \infty$? As $\lambda \to 0$?

5. Two datasets: Dataset A has $n=1000$, $p=5$. Dataset B has $n=100000$, $p=50$. For which would you use the normal equation, and for which gradient descent? Compute the approximate operation count for each.