ordinary-least-squares

The fitting criterion for linear regression: choose the weights $\mathbf{w}$ that minimise the sum of squared residuals between predictions and observed targets. Yields a closed-form solution via the normal equation ${\mathbf{w}}^{\ast }=({\mathbf{\Phi }}^{\top }\mathbf{\Phi }{)}^{-1}{\mathbf{\Phi }}^{\top }\mathbf{y}$, computable in one matrix inversion — no iterations.

The Objective

Given training data $\{({\mathbf{x}}_i,y_i){\}}_{i=1}^N$ and a linear model ${\hat{y}}_i={\mathbf{w}}^{\top }\mathbf{\varphi }({\mathbf{x}}_i)$, define the residual for example $i$ as:

$r_i=y_i-{\hat{y}}_i=y_i-{\mathbf{w}}^{\top }\mathbf{\varphi }({\mathbf{x}}_i)$

OLS picks the $\mathbf{w}$ that minimises the sum of squared residuals:

$\boxed{{\mathbf{w}}_{\text{OLS}}=\arg \underset{\mathbf{w}}{\min }R(\mathbf{w}),R(\mathbf{w})=\sum _{i=1}^N(y_i-{\mathbf{w}}^{\top }\mathbf{\varphi }({\mathbf{x}}_i){)}^2}$

The objective is convex and quadratic in $\mathbf{w}$, so its unique minimum is found by setting $\nabla R=0$.

Deriving the Normal Equation

Stack the basis-function evaluations into the design matrix $\mathbf{\Phi }$ (an $N\times (M+1)$ matrix whose row $i$ is $\mathbf{\varphi }({\mathbf{x}}_i{)}^{\top }$). Then:

$R(\mathbf{w})=\Vert \mathbf{y}-\mathbf{\Phi }\mathbf{w}{\Vert }^2$

Differentiate w.r.t. $\mathbf{w}$ and set to zero:

${\nabla }_{\mathbf{w}}R=-2{\mathbf{\Phi }}^{\top }(\mathbf{y}-\mathbf{\Phi }\mathbf{w})=0$

Rearranging gives the normal equation:

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

Solving (assuming ${\mathbf{\Phi }}^{\top }\mathbf{\Phi }$ is invertible):

$\boxed{{\mathbf{w}}_{\text{OLS}}=({\mathbf{\Phi }}^{\top }\mathbf{\Phi }{)}^{-1}{\mathbf{\Phi }}^{\top }\mathbf{y}}$

The matrix ${\mathbf{\Phi }}^{†}=({\mathbf{\Phi }}^{\top }\mathbf{\Phi }{)}^{-1}{\mathbf{\Phi }}^{\top }$ is the Moore–Penrose pseudoinverse of $\mathbf{\Phi }$. So OLS is just ${\mathbf{w}}^{\ast }={\mathbf{\Phi }}^{†}\mathbf{y}$.

Worked Derivation (Degree-1 Polynomial)

Take $\hat{y}=w_0+w_{1}x$. The objective is:

$R(w_0,w_1)={\sum }_i(y_i-w_0-w_1{x}_i{)}^2$

Setting $\partial R/\partial w_0=0$ and $\partial R/\partial w_1=0$ gives a $2\times 2$ linear system:

${\sum }_i{y}_i=w_{0}N+w_1{\sum }_i{x}_i$ ${\sum }_i{x}_i{y}_i=w_0{\sum }_i{x}_i+w_1{\sum }_i{x}_i^2$

In matrix form:

$\left(\begin{array}{c}{\sum }_i{y}_i\\ {\sum }_i{x}_i{y}_i\end{array}\right)=\left(\begin{array}{cc}N& {\sum }_i{x}_i\\ {\sum }_i{x}_i& {\sum }_i{x}_i^2\end{array}\right)\left(\begin{array}{c}{w}_0\\ w_1\end{array}\right)$

Inverting the $2\times 2$ matrix gives the closed-form solution. The general matrix-form derivation uses exactly the same logic, just with $\mathbf{\Phi }$ in place of the $[1,\mathbf{x}]$ matrix.

MLE Equivalence — The Probabilistic Justification

OLS picks “minimise squared residuals” without saying why squared, as opposed to absolute or any other power. The justification comes from a probabilistic model: assume

$y_i={\mathbf{w}}^{\top }\mathbf{\varphi }({\mathbf{x}}_i)+{\epsilon }_i,{\epsilon }_i\sim \mathcal{N}(0,{\sigma }^2)$

Then $y_i\mid {\mathbf{x}}_i\sim \mathcal{N}({\mathbf{w}}^{\top }\mathbf{\varphi }({\mathbf{x}}_i),{\sigma }^2)$. The likelihood of the dataset is:

$\mathcal{L}(\mathbf{w},{\sigma }^2)={\prod }_{i=1}^N\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(y_i-{\mathbf{w}}^{\top }\mathbf{\varphi }({\mathbf{x}}_i){)}^2}{2{\sigma }^2}\right)$

Taking the log:

$\ln \mathcal{L}=-\frac{N}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}{\sum }_{i=1}^N(y_i-{\mathbf{w}}^{\top }\mathbf{\varphi }({\mathbf{x}}_i){)}^2$

The only $\mathbf{w}$-dependent term is the sum of squared residuals — and it appears with a negative coefficient. So:

$\arg {\max }_{\mathbf{w}}\mathcal{L}(\mathbf{w})=\arg {\min }_{\mathbf{w}}{\sum }_i(y_i-{\mathbf{w}}^{\top }\mathbf{\varphi }({\mathbf{x}}_i){)}^2={\mathbf{w}}_{\text{OLS}}$

MLE = OLS under additive Gaussian noise

Under the assumption of i.i.d. Gaussian noise, the maximum-likelihood weights are exactly the OLS weights. This is the missing justification for “why squared error” — squared error isn’t arbitrary, it’s the negative log-likelihood of a Gaussian noise model (up to constants).

The MLE for ${\sigma }^2$ falls out of the same derivation: ${\hat{\sigma }}_{\text{MLE}}^2=\text{RSS}/N$, where $\text{RSS}=\Vert \mathbf{y}-\mathbf{\Phi }{\mathbf{w}}_{\text{OLS}}{\Vert }^2$.

OLS makes no distributional assumption — but its justification does

The OLS objective itself is just “minimise squared residuals” — no statistics required. The probabilistic justification (Gaussian noise → squared loss is optimal) is layered on top. So you can use OLS without believing in Gaussian noise; you just lose the optimality argument.

What Could Go Wrong

• Ill-conditioned ${\mathbf{\Phi }}^{\top }\mathbf{\Phi }$. If columns of $\mathbf{\Phi }$ are near-collinear (e.g., redundant features, polynomial basis evaluated on a narrow range), the matrix is nearly singular. Inverting it is numerically unstable; small data changes swing $\mathbf{w}$ wildly. Use ridge regression or pseudoinverse via SVD.
• $M>N$. When the number of basis functions exceeds the number of examples, ${\mathbf{\Phi }}^{\top }\mathbf{\Phi }$ is rank-deficient and not invertible. Either reduce $M$, regularise, or use a method that doesn’t require inversion.
• Outliers. Squared loss heavily penalises large residuals. A single mislabelled point can shift the entire fit. Robust alternatives: MAE, Huber loss, RANSAC.
• Wrong noise model. If the true noise is heteroscedastic (variance depends on $\mathbf{x}$) or heavy-tailed, the MLE justification breaks down. OLS still minimises squared residuals, but is no longer the optimal estimator.

Computational Cost

• Forming ${\mathbf{\Phi }}^{\top }\mathbf{\Phi }$: $O(N{M}^2)$.
• Inverting (or solving via Cholesky): $O(M^3)$.
• Total: $O(N{M}^2+M^3)$.

For modest $M$ (≤ a few thousand), OLS is extremely fast — one shot, no learning rate. For large $M$ or $N$, gradient descent or stochastic methods scale better, at the cost of iteration.

When OLS Beats Iterative Methods

• $M$ is small enough that ${\mathbf{\Phi }}^{\top }\mathbf{\Phi }$ fits in memory and inverts cheaply.
• The data is well-conditioned (no severe collinearity).
• You want exact, deterministic weights without learning-rate tuning.

When not to use OLS:

• $M$ is huge (e.g., kernel methods with $N=M$) — the matrix doesn’t fit.
• The design matrix is ill-conditioned — use ridge regression or SVD-based pseudoinverse.
• You need online updates — gradient descent on incoming data is more natural.

Connections

• linear-regression — the model OLS fits.
• design-matrix — the matrix $\mathbf{\Phi }$ that appears in the normal equation.
• maximum likelihood estimation — the principle that OLS implements under a Gaussian noise assumption.
• gaussian-distribution — the assumed noise distribution that makes MLE = OLS.
• gradient descent — iterative alternative when the closed form is impractical.
• hessian-matrix — for quadratic objectives, $H=2{\mathbf{\Phi }}^{\top }\mathbf{\Phi }$, so Newton’s method converges in one step (and that step is exactly the normal equation).

Active Recall

Write the OLS solution in a single equation. What does ${\mathbf{\Phi }}^{†}$ stand for?

${\mathbf{w}}_{\text{OLS}}=({\mathbf{\Phi }}^{\top }\mathbf{\Phi }{)}^{-1}{\mathbf{\Phi }}^{\top }\mathbf{y}={\mathbf{\Phi }}^{†}\mathbf{y}$, where ${\mathbf{\Phi }}^{†}$ is the Moore–Penrose pseudoinverse of $\mathbf{\Phi }$.

Why is "minimise squared error" a reasonable choice of objective?

Under the assumption that observations are corrupted by additive i.i.d. Gaussian noise, the MLE for the regression weights is exactly the squared-error minimiser. So squared loss is the negative log-likelihood (up to constants) of a Gaussian noise model — minimising it is principled, not arbitrary.

True or false: OLS requires an assumption about the noise distribution.

False. OLS is purely an optimisation criterion (“minimise sum of squared residuals”) — no distributional assumption. The MLE-equivalence justification requires Gaussian noise, but you can still apply OLS without believing it.

When does the normal equation fail?

When ${\mathbf{\Phi }}^{\top }\mathbf{\Phi }$ is singular or ill-conditioned. This happens with collinear features, redundant basis functions, or when there are more features than examples ($M>N$). Mitigations: drop redundant features, regularise (ridge), or use the SVD-based pseudoinverse.