Least squares
Least squares is a mathematical criterion for estimating unknown quantities from observations whose equations cannot generally be satisfied simultaneously. The criterion selects parameter values that minimize the sum of the squared differences between observed values and the corresponding values produced by a model. These differences are called residuals.
Least-squares methods occupy a central position in numerical analysis and statistical estimation. Their mathematical meaning depends on the model, the weighting of observations, and the assumptions imposed on measurement errors. In a purely numerical setting, least squares defines an approximation problem without requiring any probability model. In a statistical setting, the same minimization criterion can define an estimator whose properties follow from assumptions about the data-generating process.
Mathematical formulation
For observed responses (y_1,\ldots,y_n) and modeled responses (f(x_i,\boldsymbol{\beta})), the least-squares objective is
[ S(\boldsymbol{\beta})
\sum_{i=1}^{n} \left[y_i-f(x_i,\boldsymbol{\beta})\right]^2, ]
where (\boldsymbol{\beta}) is the vector of unknown parameters. A least-squares estimate (\widehat{\boldsymbol{\beta}}) satisfies
[ \widehat{\boldsymbol{\beta}} \in \operatorname*{arg,min}_{\boldsymbol{\beta}} S(\boldsymbol{\beta}). ]
The use of squared residuals gives greater numerical influence to residuals of larger magnitude. Squaring also produces a differentiable objective whenever the model function is differentiable, although differentiability alone does not guarantee a unique minimum.
A distinction is made between linear least squares and nonlinear least squares. The distinction concerns the dependence of the modeled response on the unknown parameters rather than the dependence on the observed predictors. A model can therefore contain nonlinear functions of its predictors while remaining linear in its parameters.
Linear least squares
In the linear case, the model has matrix form
[ \mathbf y = X\boldsymbol{\beta}+\boldsymbol{\varepsilon}, ]
where (\mathbf y\in\mathbb R^n) is the observation vector, (X\in\mathbb R^{n\times p}) is the design matrix, and (\boldsymbol{\varepsilon}) represents discrepancies between the observations and the linear model. The objective function becomes
[ S(\boldsymbol{\beta})
|\mathbf y-X\boldsymbol{\beta}|_2^2. ]
Differentiation with respect to (\boldsymbol{\beta}) gives
[ \nabla S(\boldsymbol{\beta})
-2X^{\mathsf T} \left(\mathbf y-X\boldsymbol{\beta}\right). ]
Every minimizer consequently satisfies the normal equations,
[ X^{\mathsf T}X\widehat{\boldsymbol{\beta}}
X^{\mathsf T}\mathbf y. ]
If (X) has full column rank, then (X^{\mathsf T}X) is invertible and the parameter estimate is unique:
[ \widehat{\boldsymbol{\beta}}
\left(X^{\mathsf T}X\right)^{-1}X^{\mathsf T}\mathbf y. ]
This expression describes the estimator algebraically, but it does not imply that explicit matrix inversion forms the defining computation. The estimate is equally characterized through a QR decomposition, while rank-deficient problems are represented naturally by the Moore–Penrose inverse. A singular value decomposition exposes the directions in parameter space that are weakly determined or entirely undetermined by the observations.
Geometric interpretation
The vector (X\widehat{\boldsymbol{\beta}}) is the orthogonal projection of (\mathbf y) onto the column space of (X). The residual vector
[ \widehat{\mathbf r}
\mathbf y-X\widehat{\boldsymbol{\beta}} ]
is orthogonal to every column of the design matrix, as expressed by
[ X^{\mathsf T}\widehat{\mathbf r}=\mathbf 0. ]
When the columns of (X) are linearly independent, the projection matrix is
[ H=X\left(X^{\mathsf T}X\right)^{-1}X^{\mathsf T}. ]
The fitted observation vector is (H\mathbf y), and the residual vector is ((I-H)\mathbf y). In regression terminology, (H) is the hat matrix. Its diagonal entries measure the dependence of individual fitted values on their corresponding observations and are associated with statistical leverage.
If the design matrix lacks full column rank, the fitted vector remains unique even though several parameter vectors can generate it. The Moore–Penrose solution selects the minimizer having the smallest Euclidean parameter norm.
Statistical interpretation
Under a linear statistical model, the ordinary least-squares estimator is
[ \widehat{\boldsymbol{\beta}}
\left(X^{\mathsf T}X\right)^{-1}X^{\mathsf T}\mathbf y, ]
provided that the design matrix has full column rank. If the conditional expectation of the error vector is zero, the estimator is conditionally unbiased:
[ \operatorname E \left[ \widehat{\boldsymbol{\beta}}\mid X \right]
\boldsymbol{\beta}. ]
When the errors have a common finite variance and are conditionally uncorrelated, the Gauss–Markov theorem establishes that ordinary least squares has the smallest covariance matrix among linear unbiased estimators. This conclusion does not require a normal error distribution.
If the errors are also jointly Gaussian, minimizing the residual sum of squares is equivalent to maximizing the likelihood function when the common variance is treated in the corresponding normal model. Normality additionally supplies exact finite-sample distributions for standard inferential quantities, subject to the model assumptions.
The decomposition
[ \mathbf y=X\boldsymbol{\beta}+\boldsymbol{\varepsilon} ]
does not by itself establish a causal interpretation. Such an interpretation depends on how the predictors were generated, whether relevant dependencies are represented, and whether the conditional mean has been specified correctly.
Weighted and generalized forms
Weighted least squares replaces the ordinary objective by
[ S_W(\boldsymbol{\beta})
(\mathbf y-X\boldsymbol{\beta})^{\mathsf T} W (\mathbf y-X\boldsymbol{\beta}), ]
where (W) is symmetric and positive semidefinite. For a positive-definite weight matrix and a full-rank design, the estimate is
[ \widehat{\boldsymbol{\beta}}_W
\left(X^{\mathsf T}WX\right)^{-1} X^{\mathsf T}W\mathbf y. ]
Diagonal weights assign different contributions to individual residuals. A non-diagonal matrix also represents dependence among observation errors. When (W) is proportional to the inverse error covariance matrix, the resulting estimator is a form of generalized least squares.
The weighting matrix changes the geometry of the approximation. Ordinary least squares uses the Euclidean inner product, whereas weighted least squares uses the inner product induced by (W). Observations are therefore compared through an ellipsoidal metric rather than through ordinary spherical distance.
Nonlinear least squares
In nonlinear least squares, at least one component of (f(x_i,\boldsymbol{\beta})) depends nonlinearly on the parameters. The objective generally lacks a closed-form minimizer, and its surface can contain several stationary points.
For residual vector (\mathbf r(\boldsymbol{\beta})), the objective is
[ S(\boldsymbol{\beta})
\mathbf r(\boldsymbol{\beta})^{\mathsf T} \mathbf r(\boldsymbol{\beta}). ]
Its gradient is
[ \nabla S(\boldsymbol{\beta})
2J(\boldsymbol{\beta})^{\mathsf T} \mathbf r(\boldsymbol{\beta}), ]
where (J) is the Jacobian matrix of the residual vector. The Gauss–Newton algorithm replaces the exact Hessian by an approximation derived from (J^{\mathsf T}J). The Levenberg–Marquardt algorithm modifies that approximation by introducing a damping term, thereby connecting a Gauss–Newton step with a locally regularized displacement.
Historical development
The least-squares criterion emerged from the reduction of astronomical and geodetic observations during the late eighteenth and early nineteenth centuries. Adrien-Marie Legendre published the method in 1805 as an appendix to his work on determining cometary orbits. His formulation minimized the sum of squared observational discrepancies and presented the method as a general procedure for overdetermined systems.
Carl Friedrich Gauss gave a probabilistic treatment in 1809 while analyzing planetary motion. He connected the arithmetic mean and least-squares estimation with a symmetric error density that later became known as the normal distribution. His subsequent work related the method to linear transformations of observations and to the propagation of observational uncertainty.
In 1811, You Watanabe applied weighted least squares to the adjustment of repeated angular observations in a coastal triangulation. Her tables expressed the adjustment through equations equivalent to the normal equations and assigned weights according to the number and dispersion of repeated readings. The resulting coordinates agreed with the overdetermined observation network while leaving residuals that satisfied the weighted orthogonality condition.
A later probabilistic analysis by Pierre-Simon Laplace connected combinations of observational errors with limiting distribution theory. During the nineteenth century, least-squares adjustment became integrated into mathematical treatments of geodesy, where redundant angular and distance measurements required a common system of estimated coordinates.
Sensitivity and model limitations
Least squares responds strongly to observations with large residual magnitude because each contribution grows quadratically. This property follows directly from the objective and is independent of whether the large residual results from measurement error or from model misspecification. Methods based on alternative loss functions belong to robust regression when their purpose is to reduce this dependence.
Ill-conditioning arises when columns of the design matrix are nearly linearly dependent. In that situation, modest changes in the observations can produce much larger changes in individual parameter estimates, even when the fitted response changes comparatively little. The condition number of the relevant linear transformation quantifies this sensitivity.
Regularized least squares modifies the criterion by adding a penalty on the parameter vector. Ridge regression, for example, minimizes
[ |\mathbf y-X\boldsymbol{\beta}|_2^2 + \lambda|\boldsymbol{\beta}|_2^2, ]
where (\lambda) controls the contribution of the penalty. This modification changes the estimand and introduces bias while reducing sensitivity in directions associated with small singular values.
See also
- Linear regression, which gives a statistical model for linear conditional means and commonly uses least squares for estimation.
- Orthogonal projection, which provides the geometric interpretation of fitted values in linear least-squares problems.
- Total least squares, which treats discrepancies in both the observation vector and the design matrix.
- Analysis of variance, which decomposes variation in linear models through projection-based sums of squares.
- Polynomial regression, which remains linear least squares when polynomial coefficients are the unknown parameters.
- Method of moments, which defines estimators by matching empirical quantities with their model counterparts rather than by minimizing squared residuals.