Least Squares and Linear Regression

Least squares is the workhorse of quantitative science: given a cloud of measurements and a model, it produces the best set of parameters and, just as importantly, honest error bars on them. Its reach is enormous — calibrating a detector, fitting a spectral line, reconstructing an image, training the linear layer of a network — yet underneath every application sits one clean piece of geometry. Stack the observations into a vector $\mathbf{y}$ in $n$ -dimensional data space, and collect every fit the model can produce into a $p$ -dimensional plane, the column space of the design matrix $X$ . Least squares drops a perpendicular from $\mathbf{y}$ to that plane: the foot is the fit, the perpendicular itself is the residual. Minimizing squared error is nothing more exotic than "the shortest way to a plane is straight at it". That geometric fact has three payoffs. First, the minimization has a closed form — the normal equations $X^{\mathsf{T}}X\hat{\boldsymbol\beta}=X^{\mathsf{T}}\mathbf{y}$ — because "linear regression" means linear in the parameters, so the objective is a quadratic bowl with a unique bottom. Second, the projection is not merely convenient but optimal: the Gauss–Markov theorem certifies it as the minimum-variance unbiased estimator using only three weak assumptions about the noise, with no probability distribution assumed. Third, when the noise is Gaussian, least squares is promoted to maximum likelihood and the estimator becomes exactly normal, handing us $t$ -intervals and $\chi^2$ goodness-of-fit tests for free. Linearity refers to the parameters, not the predictors — a cubic in $x$ is still linear regression — and that distinction is the source of most of the method's power and most of its misuse.

Key takeaways

Least squares fits a model by minimizing the sum of squared residuals; for models linear in the parameters this collapses to the normal equations $X^{\mathsf{T}}X\,\hat{\boldsymbol\beta}=X^{\mathsf{T}}\mathbf{y}$ with a unique closed-form solution.
The geometric content is orthogonal projection: $\hat{\mathbf{y}}=X\hat{\boldsymbol\beta}=P\mathbf{y}$ projects the data onto the column space of $X$ , so the residual is orthogonal to every regressor, $X^{\mathsf{T}}\hat{\mathbf{r}}=0$ .
Under the Gauss–Markov assumptions — zero-mean, uncorrelated, equal-variance errors — ordinary least squares is BLUE, the Best Linear Unbiased Estimator, with covariance $\sigma^2(X^{\mathsf{T}}X)^{-1}$ ; no Gaussianity is needed.
If the errors are i.i.d. Gaussian, least squares coincides with maximum likelihood and $\hat{\boldsymbol\beta}$ is itself multivariate Gaussian, giving exact $t$ -confidence intervals and $\chi^2$ goodness-of-fit.
Unequal or correlated errors are handled by weighted and generalized least squares, $\hat{\boldsymbol\beta}=(X^{\mathsf{T}}\Sigma^{-1}X)^{-1}X^{\mathsf{T}}\Sigma^{-1}\mathbf{y}$ , which restores optimality — the physicist's $\chi^2$ fit with error bars.
Key idea: Linear least squares is orthogonal projection in disguise: $\hat{\mathbf y}=X(X^{\mathsf T}X)^{-1}X^{\mathsf T}\mathbf y=P\mathbf y$ projects the data onto $\operatorname{col}(X)$ , the residual is orthogonal to all regressors, and the Gauss–Markov theorem certifies this projection as the minimum-variance linear unbiased estimate.

Where this lives today

The normal equations are quietly everywhere. Every laboratory calibration curve — a photomultiplier's counts against a known flux, a thermocouple's voltage against temperature — is a least-squares line whose slope, intercept, and error bars come from equations (3) and (5). In particle physics and astronomy the weighted $\chi^2$ fit of equation (6) is the standard tool for combining measurements with heterogeneous error bars, and the reduced $\chi^2$ is the first number anyone checks; precision experiments such as the CODATA adjustment of the fundamental constants are giant correlated least-squares problems blending thousands of measurements through a full covariance matrix. In geophysics and medical imaging, tomographic reconstruction is a least-squares inversion so ill-conditioned that Tikhonov regularization — the ridge term $\lambda I$ — is mandatory to keep noise from exploding. The method reaches into machine learning too: the final linear read-out layer of a neural network is ordinary least squares, and Gaussian-process regression is GLS with an infinite-dimensional feature map. A modern frontier is high-dimensional regression where $p$ exceeds $n$ : there $X^{\mathsf{T}}X$ is singular, OLS fails, and sparse and regularized variants (LASSO, ridge, elastic net) take over — each still descended from the projection geometry Gauss set down.

The physical picture

Think of every measurement as one coordinate. With $n$ data points, the observation vector $\mathbf{y}=(y_1,\dots,y_n)$ is a single arrow in an $n$ -dimensional space — one axis per data point, not per variable. The model cannot reach every arrow: as the $p$ parameters range over all values, the achievable fits $X\boldsymbol\beta$ sweep out only a $p$ -dimensional plane through the origin, the column space of $X$ . Least squares asks for the point on that plane closest to $\mathbf{y}$ — the orthogonal projection. The perpendicular dropped from the data to the plane is the residual, and its foot is the fit $\hat{\mathbf{y}}$ . Because the shortest segment to a plane meets it at a right angle, the residual is orthogonal to every direction the model could have moved in — the statement $X^{\mathsf{T}}\hat{\mathbf{r}}=0$ : the model has drained every drop of linear signal, and what remains is, by hypothesis, pure noise.

Figure

Least squares as orthogonal projection: the data vector y in n dimensions drops a perpendicular onto the column space of X, the fitted vector ŷ = Py is the foot of that perpendicular, and the residual r = (I−P)y is orthogonal to every regressor, so X·r = 0

There is a second view a physicist meets far more often: the scatter plot. Fit a straight line to a set of $(x_i,y_i)$ points and the same projection is at work, only now the residuals are the little vertical stubs from each point to the line. Least squares chooses the line for which those stubs, squared and summed, are smallest — equivalently, the one whose residuals sum to zero and are uncorrelated with $x$ . The second figure is a real fit to six calibration points, carrying the actual slope, intercept, and uncertainties that fall out of the normal equations and the covariance $\sigma^2(X^{\mathsf{T}}X)^{-1}$ .

Figure

Ordinary least squares fit of a straight line to six calibration points: the normal equations give slope β₁ = 1.983 ± 0.052 and intercept β₀ = 0.093 ± 0.202, the residual scatter is s = 0.217, and the red residual stubs sum to zero

Two things carry over from the geometry to the statistics. The residual, being orthogonal to the plane, "sees" only noise — which is why the residual sum of squares cleanly estimates the noise variance $\sigma^2$ . And the error bars on $\hat{\boldsymbol\beta}$ come from pushing the noise through the linear map $\hat{\boldsymbol\beta}=(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\mathbf{y}$ : if the regressor columns are nearly parallel — collinear — then $X^{\mathsf{T}}X$ is nearly singular, its inverse blows up, and the parameters become wildly uncertain. Ill-conditioning of the design is the disease of which large error bars are the symptom.

The core equations

Least squares begins with a single scalar to minimize: the sum of squared gaps between observation and model. In matrix form it is the squared length of the residual vector, which is why the whole theory is geometric.

S(\boldsymbol\beta)=\sum_{i=1}^{n}\Bigl(y_i-\textstyle\sum_{k=1}^{p}X_{ik}\beta_k\Bigr)^2=\lVert\mathbf{y}-X\boldsymbol\beta\rVert^2 \qquad (1)

where: $y_i$ — the $i$ -th observed response; $X_{ik}$ — design-matrix entry, the $k$ -th basis function evaluated at $x_i$ ; $\beta_k$ — the $k$ -th model parameter; $n,\,p$ — number of data points and parameters, with $n\ge p$ . Reading it: $S$ is a quadratic bowl in $\boldsymbol\beta$ , with a unique minimum where its gradient vanishes. Setting that gradient to zero gives the normal equations — a $p\times p$ linear system, the computational heart of the method. Their name is geometric: they state that the residual is normal (orthogonal) to the column space.

X^{\mathsf{T}}X\,\hat{\boldsymbol\beta}=X^{\mathsf{T}}\mathbf{y} \qquad (2)

where: $X^{\mathsf{T}}X$ — the $p\times p$ Gram matrix of regressors, symmetric and positive semidefinite; $\hat{\boldsymbol\beta}$ — the least-squares estimate. Reading it: if $X$ has full column rank the Gram matrix is invertible and the fit is unique; otherwise the columns are redundant and infinitely many $\hat{\boldsymbol\beta}$ tie. Solving the system gives the estimator in closed form, and multiplying back by $X$ exposes the projection matrix $P$ that carries any data vector to its shadow on the model plane.

\hat{\boldsymbol\beta}=(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\mathbf{y},\qquad \hat{\mathbf{y}}=X\hat{\boldsymbol\beta}=P\mathbf{y} \qquad (3)

where: $P=X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$ — the hat matrix projecting onto $\operatorname{col}(X)$ , obeying $P^2=P=P^{\mathsf{T}}$ ; $\hat{\mathbf{y}}$ — the fitted values. Reading it: $P$ is idempotent (projecting twice changes nothing) and symmetric (the projection is orthogonal). The residual is what is left after projection; it lies entirely in the orthogonal complement of the model plane, and is the algebraic engine of every result that follows.

\hat{\mathbf{r}}=\mathbf{y}-\hat{\mathbf{y}}=(I-P)\mathbf{y},\qquad X^{\mathsf{T}}\hat{\mathbf{r}}=0 \qquad (4)

where: $\hat{\mathbf{r}}$ — the residual vector, living in $\operatorname{col}(X)^{\perp}$ ; $I-P$ — the complementary projection, with $\operatorname{tr}(I-P)=n-p$ . Reading it: the trace counts the residual degrees of freedom — the dimension the noise still occupies after $p$ directions were spent fitting. Propagating the noise covariance through the linear estimator gives the parameter covariance — the error bars on the fit.

\operatorname{Cov}(\hat{\boldsymbol\beta})=\sigma^2(X^{\mathsf{T}}X)^{-1},\qquad s^2=\frac{\hat{\mathbf{r}}^{\mathsf{T}}\hat{\mathbf{r}}}{n-p} \qquad (5)

where: $\sigma^2$ — the common error variance, $\operatorname{Cov}(\boldsymbol\varepsilon)=\sigma^2 I$ ; $s^2$ — the unbiased estimate of $\sigma^2$ from the residuals, with $\mathbb{E}[s^2]=\sigma^2$ . Reading it: the standard error of $\hat\beta_k$ is $s\sqrt{[(X^{\mathsf{T}}X)^{-1}]_{kk}}$ , and the divisor is $n-p$ , not $n$ , precisely because of the trace in (4). When the errors are not equal and uncorrelated, OLS is still unbiased but no longer best. Generalized least squares restores optimality by weighting with the inverse covariance — the physicist's $\chi^2$ fit.

\hat{\boldsymbol\beta}=(X^{\mathsf{T}}\Sigma^{-1}X)^{-1}X^{\mathsf{T}}\Sigma^{-1}\mathbf{y},\qquad \chi^2=(\mathbf{y}-X\boldsymbol\beta)^{\mathsf{T}}\Sigma^{-1}(\mathbf{y}-X\boldsymbol\beta) \qquad (6)

where: $\Sigma$ — the full $n\times n$ error covariance; $\chi^2$ — the objective GLS minimizes. Reading it: for independent points $\Sigma=\operatorname{diag}(\sigma_i^2)$ , so each term is weighted by $w_i=1/\sigma_i^2$ — inverse-variance weighting, the rule that a precise point counts for more. Finally, if the errors are Gaussian the estimator is exactly normal, and each parameter's studentized deviation follows a $t$ -distribution — exact confidence intervals at any sample size.

\frac{\hat\beta_k-\beta_k}{s\sqrt{[(X^{\mathsf{T}}X)^{-1}]_{kk}}}\sim t_{\,n-p} \qquad (7)

where: $t_{n-p}$ — Student's $t$ with $n-p$ degrees of freedom; $[(X^{\mathsf{T}}X)^{-1}]_{kk}$ — the $k$ -th diagonal of the inverse Gram matrix. Reading it: as $n-p\to\infty$ the $t$ -distribution approaches a standard normal, recovering the familiar $\pm1.96\sigma$ interval.

Derivation AFrom the squared-error objective to the normal equations and the projection

The goal is to turn the minimization of $S(\boldsymbol\beta)=\lVert\mathbf{y}-X\boldsymbol\beta\rVert^2$ into an explicit linear system, prove that its solution is a genuine minimum rather than a saddle, and read off the geometric meaning — that the fit is the orthogonal projection of the data onto the column space and the residual is orthogonal to every regressor. We expand the quadratic, differentiate, confirm the second-order condition, and interpret.

Step 1. Expand the norm using $\lVert\mathbf{v}\rVert^2=\mathbf{v}^{\mathsf{T}}\mathbf{v}$ , giving $S=\mathbf{y}^{\mathsf{T}}\mathbf{y}-2\boldsymbol\beta^{\mathsf{T}}X^{\mathsf{T}}\mathbf{y}+\boldsymbol\beta^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol\beta$ . Why: the cross term collapses because $\mathbf{y}^{\mathsf{T}}X\boldsymbol\beta$ is a scalar and therefore equals its own transpose $\boldsymbol\beta^{\mathsf{T}}X^{\mathsf{T}}\mathbf{y}$ , so the two cross terms combine into one.

Step 2. Differentiate with respect to the vector $\boldsymbol\beta$ and set the gradient to zero: $\nabla_{\boldsymbol\beta}S=-2X^{\mathsf{T}}\mathbf{y}+2X^{\mathsf{T}}X\boldsymbol\beta=\mathbf{0}$ . Why: at a smooth minimum the gradient vanishes; the linear term gives $-2X^{\mathsf{T}}\mathbf{y}$ and the quadratic term $2X^{\mathsf{T}}X\boldsymbol\beta$ , using the symmetry of the Gram matrix.

Step 3. Rearrange to the normal equations $X^{\mathsf{T}}X\hat{\boldsymbol\beta}=X^{\mathsf{T}}\mathbf{y}$ , and if $X$ has full column rank invert to get $\hat{\boldsymbol\beta}=(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\mathbf{y}$ . Why: full column rank makes the $p\times p$ matrix $X^{\mathsf{T}}X$ positive definite and hence invertible, so the stationary point is unique.

Step 4. Confirm it is a minimum: the Hessian $\nabla^2 S=2X^{\mathsf{T}}X$ is positive semidefinite because $\mathbf{v}^{\mathsf{T}}X^{\mathsf{T}}X\mathbf{v}=\lVert X\mathbf{v}\rVert^2\ge0$ for all $\mathbf{v}$ , and strictly positive definite under full rank. Why: a positive-definite Hessian everywhere makes $S$ strictly convex, so the single stationary point is the global minimum — no saddles, no local traps.

Step 5. Read off the geometry. The fitted vector is $\hat{\mathbf{y}}=X\hat{\boldsymbol\beta}=P\mathbf{y}$ with $P=X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$ , and the normal equations rewrite as $X^{\mathsf{T}}(\mathbf{y}-X\hat{\boldsymbol\beta})=X^{\mathsf{T}}\hat{\mathbf{r}}=0$ . Why: this last identity says every column of $X$ is orthogonal to the residual, so $\hat{\mathbf{y}}$ is the orthogonal projection of $\mathbf{y}$ onto $\operatorname{col}(X)$ — the perpendicular-dropping picture, now proved. *Verified: $P^2=P$ and $P^{\mathsf{T}}=P$ follow by direct substitution; for the six-point fit the residuals satisfy $\sum_i\hat r_i=5\times10^{-14}$ and $\sum_i x_i\hat r_i=2\times10^{-13}$ numerically — orthogonality confirmed to machine precision.*

Go deeper · Derivation B · 1 minThe Gauss–Markov theorem: ordinary least squares is BLUE

Derivation BThe Gauss–Markov theorem: ordinary least squares is BLUE

The goal is to prove that among all estimators linear in the data and unbiased, ordinary least squares has the smallest covariance — any competitor's covariance exceeds the OLS covariance by a positive-semidefinite matrix — using only that the errors have zero mean, common variance $\sigma^2$ , and zero correlation. No Gaussian assumption enters; the theorem is a statement about second moments alone, which is what makes it so powerful and so often misquoted.

Step 1. Write any competing linear estimator as $\tilde{\boldsymbol\beta}=C\mathbf{y}$ for some fixed $p\times n$ matrix $C$ , and impose unbiasedness: $\mathbb{E}[\tilde{\boldsymbol\beta}]=CX\boldsymbol\beta\stackrel{!}{=}\boldsymbol\beta$ for every $\boldsymbol\beta$ , which forces $CX=I_p$ . Why: unbiasedness must hold for all true parameter values, and a linear identity holding for all $\boldsymbol\beta$ can only be satisfied by the matrix constraint $CX=I$ .

Step 2. Split $C$ around the OLS map $A=(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$ by writing $C=A+D$ . Since $AX=I$ and $CX=I$ , the deviation obeys $DX=0$ . Why: subtracting the two constraints isolates the competitor's departure $D$ from OLS and shows it must annihilate the column space of $X$ .

Step 3. Compute the competitor's covariance using $\operatorname{Cov}(\mathbf{y})=\sigma^2 I$ : $\operatorname{Cov}(\tilde{\boldsymbol\beta})=\sigma^2 CC^{\mathsf{T}}=\sigma^2(A+D)(A+D)^{\mathsf{T}}$ . Why: for a linear map $C$ acting on data of covariance $\sigma^2 I$ , the output covariance is $\sigma^2 CC^{\mathsf{T}}$ — the standard propagation rule.

Step 4. Show the cross terms vanish: $AD^{\mathsf{T}}=(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}D^{\mathsf{T}}=(X^{\mathsf{T}}X)^{-1}(DX)^{\mathsf{T}}=0$ , and likewise $DA^{\mathsf{T}}=0$ . Why: the constraint $DX=0$ from Step 2 is exactly what kills the interference between the OLS part and the deviation, leaving a clean sum of squares.

Step 5. Conclude $\operatorname{Cov}(\tilde{\boldsymbol\beta})=\sigma^2(X^{\mathsf{T}}X)^{-1}+\sigma^2 DD^{\mathsf{T}}\succeq\sigma^2(X^{\mathsf{T}}X)^{-1}$ . Why: $DD^{\mathsf{T}}$ is positive semidefinite, so any deviation $D\ne0$ can only add variance — OLS attains the floor and is therefore the Best Linear Unbiased Estimator. *Verified: cross-checked against Riley–Hobson–Bence ch. 31 and Stone–Goldbart; the single load-bearing fact is $DX=0\Rightarrow AD^{\mathsf{T}}=0$ , and the direction of the matrix inequality follows from positive-semidefiniteness of any $DD^{\mathsf{T}}$ .*

Go deeper · Derivation C · 1 minLeast squares as maximum likelihood, and why the variance divides by n minus p

Derivation CLeast squares as maximum likelihood, and why the variance divides by n minus p

The goal is to show that when the errors are independent Gaussians of common variance, maximizing the likelihood is identical to minimizing the sum of squared residuals, so least squares inherits the full optimality of maximum likelihood; and then to expose why the maximum-likelihood variance uses divisor $n$ while the unbiased estimator uses $n-p$ . This bridges the geometric OLS with the probabilistic machinery of inference, and settles a divisor that trips up nearly everyone.

Step 1. Write the likelihood for $\boldsymbol\varepsilon\sim\mathcal{N}(\mathbf{0},\sigma^2 I)$ as a product of Gaussians, $\mathcal{L}(\boldsymbol\beta,\sigma^2)=\prod_{i=1}^{n}(2\pi\sigma^2)^{-1/2}\exp[-(y_i-\mathbf{x}_i^{\mathsf{T}}\boldsymbol\beta)^2/2\sigma^2]$ , where $\mathbf{x}_i^{\mathsf{T}}$ is the $i$ -th row of $X$ . Why: independence turns the joint density into a product, and each factor is the normal density of one residual.

Step 2. Take the logarithm: $\ln\mathcal{L}=-\tfrac{n}{2}\ln(2\pi\sigma^2)-\tfrac{1}{2\sigma^2}\sum_i(y_i-\mathbf{x}_i^{\mathsf{T}}\boldsymbol\beta)^2$ . Why: the logarithm is monotonic, preserving the maximum's location, and converts the product into a sum that isolates the $\boldsymbol\beta$ -dependence in a single term.

Step 3. Maximize over $\boldsymbol\beta$ : only the last term depends on it, and it enters with a negative sign, so $\max_{\boldsymbol\beta}\ln\mathcal{L}\iff\min_{\boldsymbol\beta}\sum_i(y_i-\mathbf{x}_i^{\mathsf{T}}\boldsymbol\beta)^2=\min_{\boldsymbol\beta}\lVert\mathbf{y}-X\boldsymbol\beta\rVert^2$ . Why: maximizing a negative multiple of the residual sum of squares is the same as minimizing that sum, so the Gaussian MLE is exactly the OLS estimate — least squares is maximum likelihood.

Step 4. Now maximize over $\sigma^2$ by setting $\partial\ln\mathcal{L}/\partial\sigma^2=0$ , which gives $\hat\sigma^2_{\mathrm{MLE}}=\tfrac{1}{n}\sum_i\hat r_i^{\,2}$ . Why: at the joint maximum the residuals are already the OLS residuals, and the stationarity condition in $\sigma^2$ yields their mean square — but with divisor $n$ .

Step 5. Correct the bias. Because the residuals satisfy $p$ linear constraints ( $X^{\mathsf{T}}\hat{\mathbf{r}}=0$ ), they live in an $(n-p)$ -dimensional space and $\mathbb{E}[\hat{\mathbf{r}}^{\mathsf{T}}\hat{\mathbf{r}}]=(n-p)\sigma^2$ , so the unbiased estimate is $s^2=\hat{\mathbf{r}}^{\mathsf{T}}\hat{\mathbf{r}}/(n-p)$ . Why: the MLE divides by $n$ and therefore underestimates $\sigma^2$ ; dividing by the residual degrees of freedom $n-p=\operatorname{tr}(I-P)$ removes the bias exactly. *Verified: the divisor identity uses $\mathbb{E}[\hat{\mathbf{r}}^{\mathsf{T}}\hat{\mathbf{r}}]=\sigma^2\operatorname{tr}(I-P)=\sigma^2(n-p)$ ; the $n$ versus $n-p$ distinction is the standard textbook result (Press et al., Numerical Recipes ch. 15) and reduces to Bessel's $n-1$ correction in the constant-fit case $p=1$ .*

Limiting cases

The single machinery of least squares specializes into all the familiar named methods according to what one assumes about the error covariance $\Sigma$ and any prior on the parameters.

Case	Result	Meaning
Equal i.i.d. errors	$\Sigma=\sigma^2 I \Rightarrow \hat{\boldsymbol\beta}=(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\mathbf{y}$	Ordinary least squares; the Gauss–Markov BLUE with covariance $\sigma^2(X^{\mathsf{T}}X)^{-1}$
Independent unequal errors	$\Sigma=\operatorname{diag}(\sigma_i^2)\Rightarrow w_i=1/\sigma_i^2$	Weighted least squares — the physicist's $\chi^2$ fit; precise points pull harder
Correlated errors	$\hat{\boldsymbol\beta}=(X^{\mathsf{T}}\Sigma^{-1}X)^{-1}X^{\mathsf{T}}\Sigma^{-1}\mathbf{y}$	Generalized least squares; equivalent to whitening by $\Sigma^{-1/2}$ then doing OLS
Single constant, $X=\mathbf{1}$	$\hat\beta=\bar y,\ \operatorname{Var}=\sigma^2/n$	Fitting a constant recovers the sample mean; the standard error of the mean falls as $1/\sqrt{n}$
Gaussian prior on parameters	$\hat{\boldsymbol\beta}=(X^{\mathsf{T}}X+\lambda I)^{-1}X^{\mathsf{T}}\mathbf{y}$	Ridge / Tikhonov regularization — the MAP estimate; tames ill-conditioned $X^{\mathsf{T}}X$

Common trapLinear regression means linear in the parameters, not in the data — and log-linearizing changes the problem

The commonest confusion is that "linear regression" requires a straight-line relationship between $y$ and $x$ . It does not. A cubic $y=\beta_0+\beta_1 x+\beta_2 x^2+\beta_3 x^3$ , a Fourier fit $y=\sum_k(a_k\cos kx+b_k\sin kx)$ , and a fit to any fixed basis functions are all linear regression, because the objective is quadratic in $\boldsymbol\beta$ — the predictors can be as curved as you like. The genuinely dangerous move is the flip side: fitting an exponential $y=A e^{bx}$ by taking logs and running OLS on $\ln y$ against $x$ . That gives a quick straight-line fit but minimizes the wrong objective — the logarithm distorts the noise, turning a symmetric error on $y$ into an asymmetric one on $\ln y$ and over-weighting small- $y$ points, so the log-fit is biased and no longer the maximum-likelihood estimate for Gaussian errors on $y$ . If the noise is on $y$ , fit the nonlinear model to $y$ directly; reserve the log trick for genuinely multiplicative noise.

Worked examples

1 · A straight-line calibration fit with error bars (fully worked)

A detector's response $y$ is measured against six known input standards $x=1,\dots,6$ , giving $y=2.1,\,3.9,\,6.3,\,7.8,\,10.2,\,11.9$ . Fit $y=\beta_0+\beta_1 x$ by ordinary least squares, and report the slope and intercept with their standard errors and a 95% confidence interval on the slope.

Step 1. Assemble the sums the normal equations need. With $n=6$ : $\sum x_i=21$ , $\sum y_i=42.2$ , $\sum x_i^2=91$ , $\sum x_i y_i=182.4$ . The Gram matrix and right-hand side are therefore $X^{\mathsf{T}}X=\begin{pmatrix}6&21\\21&91\end{pmatrix}$ and $X^{\mathsf{T}}\mathbf{y}=\begin{pmatrix}42.2\\182.4\end{pmatrix}$ .

Step 2. Invert the $2\times2$ Gram matrix. Its determinant is $\det=6\cdot91-21^2=546-441=\mathbf{105}$ , so $(X^{\mathsf{T}}X)^{-1}=\frac{1}{105}\begin{pmatrix}91&-21\\-21&6\end{pmatrix}=\begin{pmatrix}0.8667&-0.2000\\-0.2000&0.05714\end{pmatrix}$ .

Step 3. Solve $\hat{\boldsymbol\beta}=(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\mathbf{y}$ : intercept $\hat\beta_0=0.8667(42.2)-0.2000(182.4)=\mathbf{0.0933}$ and slope $\hat\beta_1=-0.2000(42.2)+0.05714(182.4)=\mathbf{1.9829}$ . The scalar formula $\hat\beta_1=[n\sum x y-\sum x\sum y]/[n\sum x^2-(\sum x)^2]=208.2/105$ agrees.

Step 4. Estimate the noise. The fits $\hat y_i=0.0933+1.9829 x_i$ give residuals $\hat r_i=0.024,\,-0.159,\,0.258,\,-0.225,\,0.192,\,-0.091$ , summing to zero. Their squared sum is $\mathrm{RSS}=0.1882$ , and with $n-p=4$ the variance estimate is $s^2=0.1882/4=0.04705$ , so $s=\mathbf{0.217}$ .

Step 5. Turn $\operatorname{Cov}(\hat{\boldsymbol\beta})=s^2(X^{\mathsf{T}}X)^{-1}$ into standard errors: $\mathrm{SE}(\hat\beta_1)=s\sqrt{0.05714}=\mathbf{0.0519}$ and $\mathrm{SE}(\hat\beta_0)=s\sqrt{0.8667}=\mathbf{0.202}$ . With $t_{0.975,\,4}=2.776$ , the 95% slope interval is $1.983\pm2.776\times0.0519=1.983\pm0.144$ , i.e. $\mathbf{[1.839,\,2.127]}$ . The fit is $y=(0.093\pm0.202)+(1.983\pm0.052)\,x$ with $s=0.217$ ; the slope is resolved to about 2.6%.

2 · Weighted least squares: combining two discordant measurements

Two independent experiments measure the same constant: experiment A gives $y_1=10.0\pm0.5$ , experiment B gives $y_2=10.6\pm0.2$ . Combine them by weighted least squares and test whether they are consistent.

Step 1. The model is a single constant, $X=\begin{pmatrix}1\\1\end{pmatrix}$ , with $\Sigma=\operatorname{diag}(0.5^2,0.2^2)$ . Generalized least squares reduces to the inverse-variance weighted mean, $\hat\beta=\sum_i w_i y_i/\sum_i w_i$ with $w_i=1/\sigma_i^2$ .

Step 2. The weights are $w_1=1/0.25=4$ and $w_2=1/0.04=25$ , so the precise measurement B carries over six times the weight of A; the total is $\sum w_i=29$ .

Step 3. The weighted mean is $\hat\beta=(4\cdot10.0+25\cdot10.6)/29=305/29=\mathbf{10.517}$ , pulled toward B — far from the naive unweighted average $10.3$ , which would bias the result.

Step 4. The combined variance is $\operatorname{Var}(\hat\beta)=1/\sum w_i=1/29=0.0345$ , so the standard error is $\sqrt{0.0345}=\mathbf{0.186}$ — tighter than either input, as combining information should be.

Step 5. Test consistency with $\chi^2=\sum_i(y_i-\hat\beta)^2/\sigma_i^2=(10.0-10.517)^2/0.25+(10.6-10.517)^2/0.04=1.07+0.17=\mathbf{1.24}$ on $n-p=1$ degree of freedom. A reduced $\chi^2$ near 1 signals a good fit. The combined value is $10.52\pm0.19$ with $\chi^2=1.24$ for 1 dof — the experiments agree.

Physics of the portrait

Least squares arrived twice, was fought over, and then quietly became the substrate of quantitative science — a small drama of priority between Legendre, who published first, and Gauss, who published the deep theory and claimed to have used it years earlier to recover the lost asteroid Ceres.

Origins · 1 minCarl Friedrich Gauss and how this was discovered

Carl Friedrich Gauss

1777–1855 · the method of least squares and the Gaussian theory of errors

Gauss stands at the origin of least squares as both a computational recipe and a statistical theory. As a young man in Göttingen he used the method to predict where the newly discovered and then-lost dwarf planet Ceres would reappear in 1801 — a public triumph that made his name — though he did not publish the technique until his 1809 Theoria Motus. There he did far more than fit orbits: he asked what error distribution would make the arithmetic mean the most probable value, and derived the bell-shaped normal law now called the Gaussian, tying least squares to a theory of measurement error. Adrien-Marie Legendre had published the bare algorithm first, in an 1805 appendix on comet orbits, coining the name méthode des moindres carrés; a bitter priority dispute followed, and both men have a fair claim — Legendre to the method, Gauss to its probabilistic foundation. In the Theoria Combinationis of 1821–1823 Gauss went further, proving the minimum-variance optimality of least squares without assuming the normal law — the result later rediscovered by Andrei Markov and now bearing both names. Laplace meanwhile tied the edifice to his central limit theorem, explaining why the Gaussian error law is so ubiquitous: errors that are sums of many small independent causes tend to it. The same normal equations Gauss wrote by hand are now solved billions of times a day inside modern computers.

History

1805

— Legendre publishes the method of least squares in an appendix on comet orbits, naming it moindres carrés

1809

— Gauss, in Theoria Motus, links least squares to the normal error law and claims priority from 1795

1810

— Laplace ties least squares to the central limit theorem, explaining the ubiquity of the Gaussian

1823

— Gauss proves the minimum-variance property in Theoria Combinationis without assuming Gaussianity

1935

— Aitken casts generalized least squares for correlated errors in modern matrix form

Going further

Texts. Stone & Goldbart, Mathematics for Physics, chapter on probability and data analysis · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., least-squares fitting section · Hassani, Mathematical Physics, on inner-product spaces and projection · Byron & Fuller, Mathematics of Classical and Quantum Physics, on linear vector spaces · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed., ch. 31 · Press et al., Numerical Recipes, ch. 15 (Modeling of Data).

Provenance. All equations are standard linear-model theory and were re-verified here: (1) the normal equations from $\nabla_{\boldsymbol\beta}\lVert\mathbf{y}-X\boldsymbol\beta\rVert^2=0$ ; (2) $\operatorname{Cov}(\hat{\boldsymbol\beta})=\sigma^2(X^{\mathsf{T}}X)^{-1}$ by error propagation; (3) the Gauss–Markov BLUE argument via $DX=0$ ; (4) the MLE $=$ OLS equivalence and the $n$ versus $n-p$ divisor. The worked numbers were checked numerically: the six-point fit gives $\det(X^{\mathsf{T}}X)=105$ , $\hat\beta_1=1.9829$ , $\hat\beta_0=0.0933$ , $s=0.217$ , $\mathrm{SE}(\hat\beta_1)=0.0519$ , with residuals orthogonal to $\mathbf{1}$ and $\mathbf{x}$ to machine precision; the weighted mean gives $10.517\pm0.186$ with $\chi^2=1.24$ on 1 dof.

Prerequisites. Expectation, Variance and Moments · Multivariate Gaussian Distribution · Random Variables and Probability Distributions · Joint, Marginal and Conditional Distributions. Leads to. Central Limit Theorem · Conditional Probability and Bayes' Theorem.

Ask the tutor: Ask me to run a least-squares fit on your own data, derive any line above step by step, or show exactly why the residual must be orthogonal to the column space.

Is "linear regression" the same as "fitting a straight line"?

No — and mixing them up causes real errors. "Linear" refers to linearity in the parameters $\boldsymbol\beta$ , not in the predictor $x$ . A quadratic, a cubic, a polynomial of any degree, and a sum of sines and cosines are all linear regression, because the objective $\lVert\mathbf{y}-X\boldsymbol\beta\rVert^2$ stays quadratic in $\boldsymbol\beta$ and the normal equations still apply — you just put those curved functions of $x$ into the columns of $X$ . What makes a model nonlinear regression is the parameters entering non-linearly — an unknown rate inside an exponential, $y=Ae^{bx}$ — which has no closed-form solution and needs iterative methods.

Why square the residuals instead of taking absolute values?

Squaring has three justifications. Geometrically, minimizing a sum of squares is minimizing a Euclidean length, which yields the clean projection picture and a linear solution; absolute values give a "least absolute deviations" fit that is more robust but has no closed form. Statistically, the sum of squares is exactly the negative log-likelihood for Gaussian errors, so squaring is the maximum-likelihood choice when the noise is normal — which the central limit theorem often makes a good approximation. Algebraically, the squared objective is differentiable everywhere, giving the tidy normal equations, whereas the absolute value has a kink at zero. The price is sensitivity to outliers, since a squared residual grows fast; when outliers dominate, robust alternatives are preferred.

What is the hat matrix and why does it matter?

The hat matrix $P=X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$ "puts the hat on $\mathbf{y}$ ," turning data into fitted values, $\hat{\mathbf{y}}=P\mathbf{y}$ . It is the orthogonal projector onto $\operatorname{col}(X)$ , hence $P^2=P$ (projecting twice does nothing new) and $P^{\mathsf{T}}=P$ (the projection is perpendicular). Its diagonal entries $P_{ii}$ , the leverages, measure how much each point pulls its own fitted value — a leverage near 1 dominates the fit locally and deserves scrutiny. Its trace equals $p$ , so $\operatorname{tr}(I-P)=n-p$ is the residual degrees of freedom that sets the divisor in $s^2$ . In one matrix, $P$ encodes the fit, the residuals, the degrees of freedom, and the influence of each observation.

When must I use weighted or generalized least squares instead of ordinary OLS?

Whenever the error covariance is not proportional to the identity. If the points have different known uncertainties — heteroscedasticity, common in physics where each measurement carries its own error bar — weight each by $w_i=1/\sigma_i^2$ so precise points count more; this is the ordinary $\chi^2$ fit. If the errors are correlated — a shared systematic, or time-series autocorrelation — you need full generalized least squares with the covariance matrix $\Sigma$ . Ordinary least squares is not wrong here: by Gauss–Markov it stays unbiased. But it is no longer minimum-variance, so it throws away precision, and worse, its naive error bars $\sigma^2(X^{\mathsf{T}}X)^{-1}$ are then incorrect. The right $\Sigma$ restores both optimality and honest uncertainties.

How do I know if my fit is any good?

The residuals are the primary diagnostic — look at them, do not just trust a single number. For a weighted fit with known error bars, compute the reduced chi-squared $\chi^2_\nu=\chi^2_{\min}/(n-p)$ : a value near 1 means the model fits within the stated errors, $\chi^2_\nu\gg1$ signals a wrong model or underestimated errors, and $\chi^2_\nu\ll1$ hints at overestimated errors or overfitting. But a good $\chi^2_\nu$ is necessary, not sufficient: always plot the residuals against $x$ and against the fitted values. Structure in that plot — a trend, a curve, a fan-shape — reveals a missing term, non-linearity, or heteroscedasticity that a summary statistic hides. A good fit shows residuals scattered randomly about zero.

What happens when there are more parameters than data points?

Then $p>n$ , the columns of $X$ cannot all be independent, $X^{\mathsf{T}}X$ is singular, and ordinary least squares breaks down: infinitely many parameter vectors fit the data exactly with zero residual, so the problem is under-determined and the naive inverse does not exist. This is the norm in modern high-dimensional statistics and machine learning, and the cures all inject extra information. Ridge (Tikhonov) regression adds $\lambda I$ to make $X^{\mathsf{T}}X+\lambda I$ invertible, biasing the estimate toward zero — equivalently a Gaussian prior on the parameters. LASSO adds an $\ell_1$ penalty that drives most parameters to exactly zero, selecting a sparse subset. Each replaces the missing data constraints with a modeling assumption, and which one is right depends on what you believe about the parameters.