The normal equation

The Closed-Form Solution

The normal equation is something remarkable: a single matrix formula that computes the exact optimal weights for linear regression in one step, with no iterations, no learning rate to tune, and no convergence to worry about. Understanding it also reveals when gradient descent is the better choice — and why that is almost always in practice.

From the previous lesson, the gradient of MSE loss with respect to weights is:

Set this to zero to find the minimum:

If $\mathbf{X}^\top\mathbf{X}$ is invertible, multiply both sides on the left by $(\mathbf{X}^\top\mathbf{X})^{-1}$:

This is the Normal Equation. It gives the exact optimal weights for multi-feature linear regression in a single computation - no iterations, no learning rate, no approximation.

Anatomy of the Formula

Before the technical breakdown, here is the intuition in plain terms. $\mathbf{X}^\top$ is just $\mathbf{X}$ turned on its side — rows become columns. Multiplying $\mathbf{X}^\top$ by $\mathbf{X}$ creates a compact $p \times p$ matrix that records how every pair of features relates to each other across all training examples (i.e., how correlated they are). Multiplying $\mathbf{X}^\top$ by $\mathbf{y}$ captures how each feature relates to the target. The inverse $(\mathbf{X}^\top\mathbf{X})^{-1}$ then "divides out" those inter-feature correlations — exactly like isolating a variable in algebra — leaving the weights that best explain $\mathbf{y}$.

Read it right to left to understand what each piece computes:

$\mathbf{X}^\top\mathbf{y}$ - a $p \times 1$ vector. Entry $j$ is $\sum_i x_{ij} y_i$: how much feature $j$ co-varies with the target across all training examples. Features that correlate strongly with the target tend to get large weights.

$\mathbf{X}^\top\mathbf{X}$ - a $p \times p$ matrix. Entry $(j, k)$ is $\sum_i x_{ij} x_{ik}$: how correlated features $j$ and $k$ are with each other. Diagonal entry $(j, j) = \sum_i x_{ij}^2$ is the total squared magnitude of feature $j$.

The Why features must be linearly independent ⓘ is the core requirement for the Normal Equation to work.

$(\mathbf{X}^\top\mathbf{X})^{-1}$ - inverts the feature correlation matrix. This is the critical step: it adjusts for redundancy between features. If two features are highly correlated, you should not give both large weights - the inverse handles this by "dividing out" shared information between features.

When the Normal Equation Excels

For small to medium datasets with moderate feature counts ($p$ up to a few thousand):

• One computation - no iterations, no learning rate to tune, no convergence checking
• Exact - not an approximation; the true optimal weights
• Stable - production implementations use QR decomposition ($\texttt{numpy.linalg.lstsq}$ internally) rather than explicit matrix inversion, which is numerically more robust

The QR decomposition vs. explicit inversion ⓘ is why numpy.linalg.solve is preferred over numpy.linalg.inv.

When the Normal Equation Fails

Too many features: inverting a $p \times p$ matrix costs $O(p^3)$. For $p = 10{,}000$: $10^{12}$ operations. For $p = 1{,}000{,}000$: $10^{18}$. Modern NLP models have billions of parameters - the Normal Equation is not an option.

Perfectly correlated features: if two features are identical (or one is a linear combination of others), $\mathbf{X}^\top\mathbf{X}$ is singular. Fix: add L2 regularization, solving $(\mathbf{X}^\top\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^\top\mathbf{y}$ instead. For any \lambda > 0, this matrix is always invertible.

Non-linear models: neural networks have non-convex loss surfaces. There is no algebraic structure to "set derivative to zero and solve." Gradient descent is mandatory.