The Normal Equations vs Iterative Solvers

Ordinary least squares regression has an exact, closed-form solution: given a feature matrix $X$ and target vector $y$, the coefficient vector that minimizes squared error is $\hat\beta = (X^TX)^{-1}X^Ty$ — the normal equations. In plain English, this says: find the coefficients such that the model's residuals are exactly uncorrelated with every input feature, which is what "minimizing squared error" reduces to algebraically. For small problems this is fast and exact, computed in a single shot with no iteration and no convergence tolerance to tune.

The catch is that computing $(X^TX)^{-1}$ directly costs on the order of $p^3$ operations for $p$ features, and it can be numerically unstable when features are highly correlated (the matrix $X^TX$ becomes nearly singular, amplifying floating-point error). For a regression with thousands of features, or one solved repeatedly inside a larger optimization loop, this becomes both slow and unreliable. Iterative solvers — gradient descent, conjugate gradient, or stochastic gradient descent — instead start from a guess and repeatedly nudge the coefficients toward the minimum, never explicitly forming or inverting $X^TX$. They trade an exact answer in one step for an approximate answer reached progressively, at a cost per step that's far cheaper and that scales better to large, sparse, or ill-conditioned problems.

The practical rule of thumb: for a modest number of well-conditioned features, solve the normal equations directly (via a stable factorization like QR, not literal matrix inversion); for large-scale or ill-conditioned regressions, an iterative solver reaches a good-enough answer far faster and more robustly.