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.

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}$.

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)

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]$

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)

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

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}")

np.linalg.inv computes the full matrix inverse, then multiplies — two operations, twice the floating-point rounding error. np.linalg.solve uses LU decomposition to directly solve the linear system — more numerically stable, especially when $X^{\top }X$ is nearly singular.

Verify: Prediction and Residuals

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

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.

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

OLS Formula Summary

• Normal equations: $X^{\top }X\mathbf{w}=X^{\top }\mathbf{y}$
• Closed form: ${\mathbf{w}}^{\ast }=(X^{\top }X{)}^{-1}{X}^{\top }\mathbf{y}$
• Numerically stable: use np.linalg.solve(XtX, Xty) over np.linalg.inv(XtX) @ Xty
• Geometric interpretation: ${\mathbf{w}}^{\ast }$ projects $\mathbf{y}$ onto the column space of $X$

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. Beyond those thresholds, the $p^3$ inversion becomes the bottleneck.

The deeper limitation: the normal equation tells you the optimal weights for the MSE objective, but it can't tell you if those weights are statistically reliable. When features are nearly collinear, the coefficients are mathematically correct but practically unstable — small data changes cause large coefficient changes. That's the problem Ridge regularization addresses, not by changing the optimal solution, but by changing the objective.

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.