Linear least squares (mathematics)

From Wikipedia, the free encyclopedia
Jump to: navigation, search
This article is about the mathematics that underlie curve fitting using linear least squares. For statistical regression analysis using least squares, see linear regression. For linear regression on a single variable, see simple linear regression. For other uses, see ordinary least squares and regression analysis.

In statistics and mathematics, linear least squares is an approach fitting a mathematical or statistical model to data in cases where the idealized value provided by the model for any data point is expressed linearly in terms of the unknown parameters of the model. The resulting fitted model can be used to summarize the data, to predict unobserved values from the same system, and to understand the mechanisms that may underlie the system.

Mathematically, linear least squares is the problem of approximately solving an overdetermined system of linear equations, where the best approximation is defined as that which minimizes the sum of squared differences between the data values and their corresponding modeled values. The approach is called linear least squares since the assumed function is linear in the parameters to be estimated. Linear least squares problems are convex and have a closed-form solution that is unique, provided that the number of data points used for fitting equals or exceeds the number of unknown parameters, except in special degenerate situations. In contrast, non-linear least squares problems generally must be solved by an iterative procedure, and the problems can be non-convex with multiple optima for the objective function. If prior distributions are available, then even an underdetermined system can be solved using the Bayesian MMSE estimator.

In statistics, linear least squares problems correspond to a particularly important type of statistical model called linear regression which arises as a particular form of regression analysis. One basic form of such a model is an ordinary least squares model. The present article concentrates on the mathematical aspects of linear least squares problems, with discussion of the formulation and interpretation of statistical regression models and statistical inferences related to these being dealt with in the articles just mentioned. See outline of regression analysis for an outline of the topic.


A plot of the data points (in red), the least squares line of best fit (in blue), and the residuals (in green).

As a result of an experiment, four (x, y) data points were obtained, (1, 6), (2, 5), (3, 7), and (4, 10) (shown in red in the picture on the right). We hope to find a line y=\beta_1+\beta_2 x that best fits these four points. In other words, we would like to find the numbers \beta_1 and \beta_2 that approximately solve the overdetermined linear system

\beta_1  +  1\beta_2 &&\; = \;&& 6 & \\
\beta_1  +  2\beta_2 &&\; = \;&& 5 & \\
\beta_1  +  3\beta_2 &&\; = \;&& 7 & \\
\beta_1  +  4\beta_2 &&\; = \;&& 10 & \\

of four equations in two unknowns in some "best" sense.

The "error", at each point, between the curve fit and the data is the difference between the right- and left-hand sides of the equations above. The least squares approach to solving this problem is to try to make as small as possible the sum of the squares of these errors; that is, to find the minimum of the function

\begin{align}S(\beta_1, \beta_2) =&
+\left[5-(\beta_1+2\beta_2)   \right]^2 \\
&+\left[7-(\beta_1 +  3\beta_2)\right]^2
+\left[10-(\beta_1  +  4\beta_2)\right]^2 \\
&= 4\beta_1^2 + 30\beta_2^2 + 20\beta_1\beta_2 - 56\beta_1 - 154\beta_2 + 210 .\end{align}

The minimum is determined by calculating the partial derivatives of S(\beta_1, \beta_2) with respect to \beta_1 and \beta_2 and setting them to zero

\frac{\partial S}{\partial \beta_1}=0=8\beta_1 + 20\beta_2 -56
\frac{\partial S}{\partial \beta_2}=0=20\beta_1 + 60\beta_2 -154.

This results in a system of two equations in two unknowns, called the normal equations, which give, when solved


and the equation y=3.5+1.4x of the line of best fit. The residuals, that is, the discrepancies between the y values from the experiment and the y values calculated using the line of best fit are then found to be 1.1, -1.3, -0.7, and 0.9 (see the picture on the right). The minimum value of the sum of squares of the residuals is S(3.5, 1.4)=1.1^2+(-1.3)^2+(-0.7)^2+0.9^2=4.2.

More generally, one can have n regressors x_j, and a linear model

y = \beta_1 + \sum_{i=2}^{n+1} \beta_i x_{i-1} .

Using a quadratic model[edit]

The result of fitting a quadratic function y=\beta_1+\beta_2x+\beta_3x^2\, (in blue) through a set of data points (x_i, y_i) (in red). In linear least squares the function need not be linear in the argument x, but only in the parameters \beta_j that are determined to give the best fit.

Importantly, in "linear least squares", we are not restricted to using a line as the model as in the above example. For instance, we could have chosen the restricted quadratic model y=\beta_1 x^2. This model is still linear in the \beta_1 parameter, so we can still perform the same analysis, constructing a system of equations from the data points:

6 &&\; = \beta_1 (1)^2 \\
5 &&\; = \beta_1 (2)^2 \\
7 &&\; = \beta_1 (3)^2 \\
10 &&\; = \beta_1 (4)^2 \\

The partial derivatives with respect to the parameters (this time there is only one) are again computed and set to 0:

\frac{\partial S}{\partial \beta_1} = 0 = 708 \beta_1 - 498

and solved

\beta_1 = .703

leading to the resulting best fit model y = .703 x^2

The general problem[edit]

Consider an overdetermined system

\sum_{j=1}^{n} X_{ij}\beta_j = y_i,\ (i=1, 2, \dots, m),

of m linear equations in n unknown coefficients, β1,β2,…,βn, with m > n. (Note: for a linear model as above, not all of X contains information on the data points. The first column is populated with ones, X_{j1} = 1, only the other columns contain actual data, and n = number of regressors + 1). This can be written in matrix form as

\mathbf {X} \boldsymbol {\beta} = \mathbf {y},


\mathbf {X}=\begin{bmatrix}
X_{11} & X_{12} & \cdots & X_{1n} \\
X_{21} & X_{22} & \cdots & X_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
X_{m1} & X_{m2} & \cdots & X_{mn}
\end{bmatrix} , 
\qquad \boldsymbol \beta = \begin{bmatrix} 
\beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{bmatrix} , 
\qquad \mathbf y = \begin{bmatrix} 
y_1 \\ y_2 \\ \vdots \\ y_m

Such a system usually has no solution, so the goal is instead to find the coefficients β which fit the equations "best," in the sense of solving the quadratic minimization problem

\hat{\boldsymbol{\beta}} = \underset{\boldsymbol{\beta}}{\operatorname{arg\,min}}\,S(\boldsymbol{\beta}),

where the objective function S is given by

S(\boldsymbol{\beta}) = \sum_{i=1}^{m}\bigl| y_i - \sum_{j=1}^{n} X_{ij}\beta_j\bigr|^2 = \bigl\|\mathbf y - \mathbf X \boldsymbol \beta \bigr\|^2.

A justification for choosing this criterion is given in properties below. This minimization problem has a unique solution, provided that the n columns of the matrix X are linearly independent, given by solving the normal equations

(\mathbf X^{\rm T} \mathbf X )\hat{\boldsymbol{\beta}}= \mathbf X^{\rm T} \mathbf y.

The matrix \mathbf X^{\rm T} \mathbf X is known as the Gramian matrix of \mathbf X, which possesses several nice properties such as being a positive semi-definite matrix.

Following is one generalized example which shows how to find a best fit line using least squares method (in two dimension). Later, a MATLAB code implementation of the same is also written which can help one visualize graphically.

Let us consider 3 points (x1,y1), (x2,y2) and (x3,y3) in a two dimensional frame. Now, assume that the best fit line is in the form: y = C + Dx. Substituting the three points in this equation will give us three equations.

These three equations can be written in the form of a matrix equation of the type Ax=b. Projection of matrix A on b gives us x. This can be represented by the formula :

(\mathbf A^{\rm T} \mathbf A )\hat{\boldsymbol{\beta}}= \mathbf A^{\rm T} \mathbf y.

We get beta from this expression; i.e. \hat{\boldsymbol{\beta}} is x. Hence, this follows that we have obtained C and D of the equation y = C + Dx. Hence, the best fit line obtained.

Matlab code implementation for better understanding in visuals:

 %MATLAB code for finding the best fit line using least squares method
 x=input('enter a')                      %input in the form of matrix, rows contain points
     a=[1,x(1,1);1,x(2,1);1,x(3,1)]          %forming A of Ax=b
     b=[x(1,2);x(2,2);x(3,2)]                %forming b of Ax=b
     yy=inv(transpose(a)*a)*transpose(a)*b   %computing projection of matrix A on b, giving x
 %plotting the best fit line
     %plotting the points(data) for which we found the best fit line
 hold on
 hold on
 hold on
 hold off

Python code using the same variable naming as the Matlab code:

import numpy as np
import matplotlib.pyplot as plt

m = 3
x = np.random.rand(m, 2) * 10
a = np.matrix([np.ones(m), x[:,0]]).T
b = np.matrix(x[:,1]).T
yy = (a.T * a).I * a.T * b
xx = np.linspace(1, 10, 50)
y = np.array(yy[0] + yy[1] * xx)

plt.plot(xx, y.T, color='r')
plt.scatter(x[:,0], x[:,1])

Derivation of the normal equations[edit]

Define the ith residual to be

r_i= y_i - \sum_{j=1}^{n} X_{ij}\beta_j.

Then S can be rewritten

S = \sum_{i=1}^m r_i^2.

S is minimized when its gradient vector is zero. (This follows by definition: if the gradient vector is not zero, there is a direction in which we can move to minimize it further - see maxima and minima.) The elements of the gradient vector are the partial derivatives of S with respect to the parameters:

\frac{\partial S}{\partial \beta_j}=2\sum_{i = 1}^m r_i\frac{\partial r_i}{\partial \beta_j} \ (j=1,2,\dots, n).

The derivatives are

\frac{\partial r_i}{\partial \beta_j}=-X_{ij}.

Substitution of the expressions for the residuals and the derivatives into the gradient equations gives

\frac{\partial S}{\partial \beta_j} = 2\sum_{i=1}^{m} \left( y_i-\sum_{k=1}^{n} X_{ik}\beta_k \right) (-X_{ij})\ (j=1,2,\dots, n).

Thus if \hat \beta minimizes S, we have

2\sum_{i=1}^{m} \left( y_i-\sum_{k=1}^{n} X_{ik}\hat \beta_k \right) (-X_{ij}) = 0\ (j=1,2,\dots, n).

Upon rearrangement, we obtain the normal equations:

\sum_{i=1}^{m}\sum_{k=1}^{n} X_{ij}X_{ik}\hat \beta_k=\sum_{i=1}^{m} X_{ij}y_i\ (j=1,2,\dots, n).

The normal equations are written in matrix notation as

(\mathbf X^\mathrm{T} \mathbf X) \hat{\boldsymbol{\beta}} = \mathbf X^\mathrm{T} \mathbf y (where XT is the matrix transpose of X).

The solution of the normal equations yields the vector \hat{\boldsymbol{\beta}} of the optimal parameter values.

Derivation directly in terms of matrices[edit]

The normal equations can be derived directly from a matrix representation of the problem as follows. The objective is to minimize

= \bigl\|\mathbf y - \mathbf X \boldsymbol \beta \bigr\|^2 
= (\mathbf y-\mathbf X \boldsymbol \beta)^{\rm T}(\mathbf y-\mathbf X \boldsymbol \beta) 
= \mathbf y ^{\rm T} \mathbf y - \boldsymbol \beta ^{\rm T} \mathbf X ^{\rm T} \mathbf y - \mathbf y ^{\rm T} \mathbf X \boldsymbol \beta + \boldsymbol \beta ^{\rm T} \mathbf X ^{\rm T} \mathbf X \boldsymbol \beta .

Note that  :( \boldsymbol \beta ^{\rm T} \mathbf X ^{\rm T} \mathbf y ) ^{\rm T} = \mathbf y ^{\rm T} \mathbf X \boldsymbol \beta has the dimension 1x1 (the number of columns of \mathbf y), so it is a scalar and equal to its own transpose, hence \boldsymbol \beta ^{\rm T} \mathbf X ^{\rm T} \mathbf y = \mathbf y ^{\rm T} \mathbf X \boldsymbol \beta and the quantity to minimize becomes

S(\boldsymbol{\beta}) = \mathbf y ^{\rm T} \mathbf y - 2\boldsymbol \beta ^{\rm T} \mathbf X ^{\rm T} \mathbf y + \boldsymbol \beta ^{\rm T} \mathbf X ^{\rm T} \mathbf X \boldsymbol \beta .

Differentiating this with respect to \boldsymbol \beta and equating to zero to satisfy the first-order conditions gives

- \mathbf X^{\rm T} \mathbf y+ (\mathbf X^{\rm T} \mathbf X ){\boldsymbol{\beta}} = 0,

which is equivalent to the above-given normal equations. A sufficient condition for satisfaction of the second-order conditions for a minimum is that \mathbf X have full column rank, in which case \mathbf X^{\rm T} \mathbf X is positive definite.

Derivation without calculus[edit]

When \mathbf X^{\rm T} \mathbf X is positive definite, the formula for the minimizing value of  \boldsymbol \beta can be derived without the use of derivatives. The quantity

S(\boldsymbol{\beta}) = \mathbf y ^{\rm T} \mathbf y - 2\boldsymbol \beta ^{\rm T} \mathbf X ^{\rm T} \mathbf y + \boldsymbol \beta ^{\rm T} \mathbf X ^{\rm T} \mathbf X \boldsymbol \beta

can be written as

 \langle \boldsymbol \beta, \boldsymbol \beta \rangle - 2\langle \boldsymbol \beta, (\mathbf X^{\rm T} \mathbf X)^{-1}\mathbf X ^{\rm T} \mathbf y \rangle + \langle(\mathbf X^{\rm T} \mathbf X)^{-1}\mathbf X ^{\rm T} \mathbf y,(\mathbf X^{\rm T} \mathbf X)^{-1}\mathbf X ^{\rm T} \mathbf y \rangle+ C,

where  C depends only on  \mathbf y and  \mathbf X , and  \langle \cdot, \cdot \rangle is the inner product defined by

 \langle x, y \rangle = x ^{\rm T} (\mathbf X^{\rm T} \mathbf X) y.

It follows that  S(\boldsymbol{\beta}) is equal to

\langle \boldsymbol \beta - (\mathbf X^{\rm T} \mathbf X)^{-1}\mathbf X ^{\rm T} \mathbf y,\boldsymbol \beta - (\mathbf X^{\rm T} \mathbf X)^{-1}\mathbf X ^{\rm T} \mathbf y \rangle+ C

and therefore minimized exactly when

\boldsymbol \beta - (\mathbf X^{\rm T} \mathbf X)^{-1}\mathbf X ^{\rm T} \mathbf y = 0.


A general approach to the least squares problem \operatorname{\,min} \, \big\|\mathbf y - X \boldsymbol \beta \big\|^2 can be described as follows. Suppose that we can find an n by m matrix S such that XS is an orthogonal projection onto the image of X. Then a solution to our minimization problem is given by

\boldsymbol \beta  =  S   \mathbf y

simply because

  X \boldsymbol \beta  =  X ( S   \mathbf y) = (X  S)  \mathbf y

is exactly a sought for orthogonal projection of   \mathbf y onto an image of X (see the picture below and note that as explained in the next section the image of X is just a subspace generated by column vectors of X). A few popular ways to find such a matrix S are described below.

Inverting the matrix of the normal equations[edit]

The algebraic solution of the normal equations can be written as

 \hat{\boldsymbol{\beta}} = (\mathbf X^ {\rm T} \mathbf X )^{-1} \mathbf X^ {\rm T} \mathbf y 
= \mathbf X^+ \mathbf y

where X + is the Moore–Penrose pseudoinverse of X. Although this equation is correct, and can work in many applications, it is not computationally efficient to invert the normal equations matrix (the Gramian matrix). An exception occurs in numerical smoothing and differentiation where an analytical expression is required.

If the matrix XTX is well-conditioned and positive definite, implying that it has full rank, the normal equations can be solved directly by using the Cholesky decomposition RTR, where R is an upper triangular matrix, giving:

 R^{\rm T} R \hat{\boldsymbol{\beta}} =  X^{\rm T} \mathbf y.

The solution is obtained in two stages, a forward substitution step, solving for z:

 R^{\rm T} \mathbf z = X^{\rm T} \mathbf y,

followed by a backward substitution, solving for \hat{\boldsymbol{\beta}}

R \hat{\boldsymbol{\beta}}= \mathbf z.

Both substitutions are facilitated by the triangular nature of R.

See example of linear regression for a worked-out numerical example with three parameters.

Orthogonal decomposition methods[edit]

Orthogonal decomposition methods of solving the least squares problem are slower than the normal equations method but are more numerically stable because they avoid forming the product XTX.

The residuals are written in matrix notation as

\mathbf r= \mathbf y - X \hat{\boldsymbol{\beta}}.

The matrix X is subjected to an orthogonal decomposition, e.g., the QR decomposition as follows.

X=QR \ ,

where Q is an m×n orthogonal matrix (QTQ=I) and R is an n×n upper triangular matrix with r_{ii}>0.

The residual vector is left-multiplied by QT.

Q^{\rm T} \mathbf r = Q^{\rm T} \mathbf y - \left( Q^{\rm T} Q \right) R \hat{\boldsymbol{\beta}}= \begin{bmatrix}
\left(Q^{\rm T} \mathbf y \right)_n - R \hat{\boldsymbol{\beta}}  \\
\left(Q^{\rm T} \mathbf y \right)_{m-n} 
= \begin{bmatrix}
\mathbf u \\
\mathbf v

Because Q is orthogonal, the sum of squares of the residuals, s, may be written as:

s = \|\mathbf r \|^2 = \mathbf r^{\rm T} \mathbf r = \mathbf r^{\rm T} Q Q^{\rm T} \mathbf r = \mathbf u^{\rm T} \mathbf u + \mathbf v^{\rm T} \mathbf v

Since v doesn't depend on β, the minimum value of s is attained when the upper block, u, is zero. Therefore the parameters are found by solving:

 R \hat{\boldsymbol{\beta}} =\left(Q^{\rm T} \mathbf y \right)_n.

These equations are easily solved as R is upper triangular.

An alternative decomposition of X is the singular value decomposition (SVD)[1]

 X = U \Sigma V^{\rm T} \ ,

where U is m by m orthogonal matrix, V is n by n orthogonal matrix and \Sigma is an m by n matrix with all its elements outside of the main diagonal equal to 0. The pseudoinverse of \Sigma is easily obtained by inverting its non-zero diagonal elements and transposing. Hence,

  \mathbf X \mathbf X^+ = U \Sigma V^{\rm T}  V \Sigma^+  U^{\rm T} =  U P  U^{\rm T},

where P is obtained from \Sigma by replacing its non-zero diagonal elements with ones. Since (\mathbf X \mathbf X^+)^* = \mathbf X \mathbf X^+ (the property of pseudoinverse), the matrix U P U^{\rm T} is an orthogonal projection onto the image (column-space) of X. In accordance with a general approach described in the introduction above (find XS which is an orthogonal projection),

 S = \mathbf X^+ ,

and thus,

 \beta = V\Sigma^+ U^{\rm T} \mathbf y

is a solution of a least squares problem. This method is the most computationally intensive, but is particularly useful if the normal equations matrix, XTX, is very ill-conditioned (i.e. if its condition number multiplied by the machine's relative round-off error is appreciably large). In that case, including the smallest singular values in the inversion merely adds numerical noise to the solution. This can be cured with the truncated SVD approach, giving a more stable and exact answer, by explicitly setting to zero all singular values below a certain threshold and so ignoring them, a process closely related to factor analysis.

Properties of the least-squares estimators[edit]

The residual vector, y-X \hat{\boldsymbol\beta}, which corresponds to the solution of a least squares system, y=X\boldsymbol \beta +\epsilon, is orthogonal to the column space of the matrix X.

The gradient equations at the minimum can be written as

(\mathbf y - X \hat{\boldsymbol{\beta}})^{\rm T} X=0.

A geometrical interpretation of these equations is that the vector of residuals, \mathbf y - X \hat{\boldsymbol{\beta}} is orthogonal to the column space of X, since the dot product (\mathbf y-X\hat{\boldsymbol{\beta}})\cdot X \mathbf v is equal to zero for any conformal vector, v. This means that \mathbf y - X \boldsymbol{\hat \beta} is the shortest of all possible vectors \mathbf{y}- X \boldsymbol \beta, that is, the variance of the residuals is the minimum possible. This is illustrated at the right.

Introducing \hat{\boldsymbol{\gamma}} and a matrix K with the assumption that a matrix [X \ K] is non-singular and KT X = 0 (cf. Orthogonal projections), the residual vector should satisfy the following equation:

\hat{\mathbf{r}} \triangleq \mathbf{y} - X \hat{\boldsymbol{\beta}} = K \hat{{\boldsymbol{\gamma}}}.

The equation and solution of linear least squares are thus described as follows:

 \mathbf{y} = \begin{bmatrix}X & K\end{bmatrix} \begin{pmatrix} \hat{\boldsymbol{\beta}} \\ \hat{\boldsymbol{\gamma}} \end{pmatrix} ,
 \begin{pmatrix} \hat{\boldsymbol{\beta}} \\ \hat{\boldsymbol{\gamma}} \end{pmatrix} = \begin{bmatrix}X & K\end{bmatrix}^{-1} \mathbf{y} = \begin{bmatrix} (X^{\rm T} X)^{-1} X^{\rm T} \\ (K^{\rm T} K)^{-1} K^{\rm T} \end{bmatrix} \mathbf{y} .

If the experimental errors, \epsilon \,, are uncorrelated, have a mean of zero and a constant variance, \sigma, the Gauss-Markov theorem states that the least-squares estimator, \hat{\boldsymbol{\beta}}, has the minimum variance of all estimators that are linear combinations of the observations. In this sense it is the best, or optimal, estimator of the parameters. Note particularly that this property is independent of the statistical distribution function of the errors. In other words, the distribution function of the errors need not be a normal distribution. However, for some probability distributions, there is no guarantee that the least-squares solution is even possible given the observations; still, in such cases it is the best estimator that is both linear and unbiased.

For example, it is easy to show that the arithmetic mean of a set of measurements of a quantity is the least-squares estimator of the value of that quantity. If the conditions of the Gauss-Markov theorem apply, the arithmetic mean is optimal, whatever the distribution of errors of the measurements might be.

However, in the case that the experimental errors do belong to a normal distribution, the least-squares estimator is also a maximum likelihood estimator.[2]

These properties underpin the use of the method of least squares for all types of data fitting, even when the assumptions are not strictly valid.


An assumption underlying the treatment given above is that the independent variable, x, is free of error. In practice, the errors on the measurements of the independent variable are usually much smaller than the errors on the dependent variable and can therefore be ignored. When this is not the case, total least squares or more generally errors-in-variables models, or rigorous least squares, should be used. This can be done by adjusting the weighting scheme to take into account errors on both the dependent and independent variables and then following the standard procedure.[3][4]

In some cases the (weighted) normal equations matrix XTX is ill-conditioned. When fitting polynomials the normal equations matrix is a Vandermonde matrix. Vandermonde matrices become increasingly ill-conditioned as the order of the matrix increases.[citation needed] In these cases, the least squares estimate amplifies the measurement noise and may be grossly inaccurate.[citation needed] Various regularization techniques can be applied in such cases, the most common of which is called ridge regression. If further information about the parameters is known, for example, a range of possible values of \mathbf{\hat{\boldsymbol{\beta}}}, then various techniques can be used to increase the stability of the solution. For example, see constrained least squares.

Another drawback of the least squares estimator is the fact that the norm of the residuals, \| \mathbf y - X\hat{\boldsymbol{\beta}} \| is minimized, whereas in some cases one is truly interested in obtaining small error in the parameter \mathbf{\hat{\boldsymbol{\beta}}}, e.g., a small value of \|{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}\|.[citation needed] However, since the true parameter {\boldsymbol{\beta}} is necessarily unknown, this quantity cannot be directly minimized. If a prior probability on \hat{\boldsymbol{\beta}} is known, then a Bayes estimator can be used to minimize the mean squared error, E \left\{ \| {\boldsymbol{\beta}} - \hat{\boldsymbol{\beta}} \|^2 \right\} . The least squares method is often applied when no prior is known. Surprisingly, when several parameters are being estimated jointly, better estimators can be constructed, an effect known as Stein's phenomenon. For example, if the measurement error is Gaussian, several estimators are known which dominate, or outperform, the least squares technique; the best known of these is the James–Stein estimator. This is an example of more general shrinkage estimators that have been applied to regression problems.

Weighted linear least squares[edit]

See also: Weighted mean

In some cases the observations may be weighted—for example, they may not be equally reliable. In this case, one can minimize the weighted sum of squares:

\underset{\boldsymbol \beta}{ \operatorname{arg\,min} }\, \sum_{i=1}^{m} w_i \left|y_i - \sum_{j=1}^{n} X_{ij}\beta_j\right|^2 = \underset{\boldsymbol \beta}{ \operatorname{arg\,min} } \, \big\|W^{1/2} (\mathbf y - X \boldsymbol \beta) \big\|^2.

where wi > 0 is the weight of the ith observation, and W is the diagonal matrix of such weights.

The weights should, ideally, be equal to the reciprocal of the variance of the measurement.[5] [6] The normal equations are then:

\left(X^{\rm T} W X \right)\hat{\boldsymbol{\beta}} = X^{\rm T} W \mathbf y.

This method is used in iteratively reweighted least squares.

Parameter errors and correlation[edit]

The estimated parameter values are linear combinations of the observed values

\hat{\boldsymbol{\beta}} = (X^{\rm T} W X)^{-1} X^{\rm T} W \mathbf y. \,

Therefore an expression for the residuals (i.e., the estimated errors in the parameters) can be obtained by error propagation from the errors in the observations. Let the variance-covariance matrix for the observations be denoted by M and that of the parameters by Mβ. Then,

M^\beta= (X^{\rm T} W X)^{-1} X^{\rm T} W M W^{\rm T} X (X^{\rm T} W^{\rm T} X)^{-1}.

When W = M −1 this simplifies to

M^\beta=(X^{\rm T} W X)^{-1}.

When unit weights are used (W = I) it is implied that the experimental errors are uncorrelated and all equal: M = σ2I, where σ2 is the variance of an observation, and I is the identity matrix. In this case σ2 is approximated by \frac{S}{m-n}, where S is the minimum value of the objective function

M^\beta=\frac{S}{m-n}(X^{\rm T} X)^{-1}.

The denominator, mn, is the number of degrees of freedom; see effective degrees of freedom for generalizations for the case of correlated observations. In all cases, the variance of the parameter \beta_i is given by M^\beta_{ii} and the covariance between parameters \beta_i and \beta_j is given by M^\beta_{ij}. Standard deviation is the square root of variance, and the correlation coefficient is given by \rho_{ij} = M^\beta_{ij}/(\sigma_i \sigma_j). These error estimates reflect only random errors in the measurements. The true uncertainty in the parameters is larger due to the presence of systematic errors which, by definition, cannot be quantified. Note that even though the observations may be un-correlated, the parameters are typically correlated.

Parameter confidence limits[edit]

Main article: Confidence interval

It is often assumed, for want of any concrete evidence but often appealing to the central limit theorem -- see Normal distribution#Occurrence -- that the error on each observation belongs to a normal distribution with a mean of zero and standard deviation \sigma. Under that assumption the following probabilities can be derived for a single scalar parameter estimate in terms of its estimated standard error se_{\beta} (given here):

68% that the interval \hat \beta \pm se_{\beta} encompasses the true coefficient value
95% that the interval \hat \beta \pm 2se_{\beta} encompasses the true coefficient value
99% that the interval \hat \beta \pm 2.5se_{\beta} encompasses the true coefficient value

The assumption is not unreasonable when m >> n. If the experimental errors are normally distributed the parameters will belong to a Student's t-distribution with m − n degrees of freedom. When m >> n Student's t-distribution approximates a normal distribution. Note, however, that these confidence limits cannot take systematic error into account. Also, parameter errors should be quoted to one significant figure only, as they are subject to sampling error.[7]

When the number of observations is relatively small, Chebychev's inequality can be used for an upper bound on probabilities, regardless of any assumptions about the distribution of experimental errors: the maximum probabilities that a parameter will be more than 1, 2 or 3 standard deviations away from its expectation value are 100%, 25% and 11% respectively.

Residual values and correlation[edit]

The residuals are related to the observations by

\mathbf{\hat r} = \mathbf y- X \hat{\boldsymbol{\beta}}= \mathbf y- H \mathbf y =  (I - H) \mathbf y

where H is the idempotent matrix known as the hat matrix:

H = X \left(X^{\rm T} W X \right)^{-1}X^{\rm T} W

and I is the identity matrix. The variance-covariance matrix of the residuals, Mr is given by

M^{\mathbf r}=\left(I-H \right) M \left(I-H \right)^{\rm T}.

Thus the residuals are correlated, even if the observations are not.

When W = M^{-1},

M^{\mathbf r}=\left(I-H \right) M.

The sum of residual values is equal to zero whenever the model function contains a constant term. Left-multiply the expression for the residuals by XT:

X^{\rm T} \hat{\mathbf r}=X^{\rm T} \mathbf y- X^{\rm T} X \hat{\boldsymbol{\beta}} = X^{\rm T} \mathbf y- (X^{\rm T} X)(X^{\rm T} X)^{-1}X^{\rm T} \mathbf y= \mathbf 0

Say, for example, that the first term of the model is a constant, so that X_{i1}=1 for all i. In that case it follows that

\sum_i^m X_{i1} \hat r_i=\sum_i^m \hat r_i=0.

Thus, in the motivational example, above, the fact that the sum of residual values is equal to zero it is not accidental but is a consequence of the presence of the constant term, α, in the model.

If experimental error follows a normal distribution, then, because of the linear relationship between residuals and observations, so should residuals,[8] but since the observations are only a sample of the population of all possible observations, the residuals should belong to a Student's t-distribution. Studentized residuals are useful in making a statistical test for an outlier when a particular residual appears to be excessively large.

Objective function[edit]

The optimal value of the objective function, found by substituting in the optimal expression for the coefficient vector, can be written as

S=\mathbf y^{\rm T} (I-H)^{\rm T} (I-H) \mathbf y= \mathbf y^{\rm T} (I-H) \mathbf y,

the latter equality holding since (IH) is symmetric and idempotent. It can be shown from this[9] that under an appropriate assignment of weights the expected value of S is m-n. If instead unit weights are assumed, the expected value of S is (m-n)\sigma^2, where \sigma^2 is the variance of each observation.

If it is assumed that the residuals belong to a normal distribution, the objective function, being a sum of weighted squared residuals, will belong to a chi-squared (\chi ^2) distribution with m-n degrees of freedom. Some illustrative percentile values of \chi ^2 are given in the following table.[10]

m-n \chi ^2 _{0.50} \chi ^2 _{0.95} \chi ^2 _{0.99}
10 9.34 18.3 23.2
25 24.3 37.7 44.3
100 99.3 124 136

These values can be used for a statistical criterion as to the goodness-of-fit. When unit weights are used, the numbers should be divided by the variance of an observation.

Constrained linear least squares[edit]

Often it is of interest to solve a linear least squares problem with an additional constraint on the solution. With constrained linear least squares, the original equation

\mathbf {X} \boldsymbol {\beta} = \mathbf {y}

must be fit as closely as possible (in the least squares sense) while ensuring that some other property of \boldsymbol {\beta} is maintained. There are often special purpose algorithms for solving such problems efficiently. Some examples of constraints are given below:

When the constraint only applies to some of the variables, the mixed problem may be solved using separable least squares by letting \mathbf {X} = [\mathbf {X_1} \mathbf {X_2} ] and \mathbf {\beta}^{\rm T} = [\mathbf {\beta_1}^{\rm T} \mathbf {\beta_2}^{\rm T}] represent the unconstrained (1) and constrained (2) components. Then substituting the least squares solution for \mathbf {\beta_1}, i.e.

\hat{\boldsymbol {\beta_1}} = \mathbf {X_1}^+ (\mathbf {y} - \mathbf {X_2} \boldsymbol {\beta_2})

back into the original expression gives (following some rearrangement) an equation that can be solved as a purely constrained problem in \mathbf {\beta_2}.

 \mathbf{P} \mathbf {X_2} \boldsymbol {\beta_2} = \mathbf{P}\mathbf {y}

where \mathbf{P}:=\mathbf{I}-\mathbf {X_1} \mathbf {X_1}^+ is a projection matrix. Following the constrained estimation of \hat{\boldsymbol {\beta_2}} the vector \hat{\boldsymbol {\beta_1}} is obtained from the expression above.

Typical uses and applications[edit]

Uses in data fitting[edit]

The primary application of linear least squares is in data fitting. Given a set of m data points y_1, y_2,\dots, y_m, consisting of experimentally measured values taken at m values x_1, x_2,\dots, x_m of an independent variable (x_i may be scalar or vector quantities), and given a model function y=f(x, \boldsymbol \beta), with \boldsymbol \beta = (\beta_1, \beta_2, \dots, \beta_n), it is desired to find the parameters \beta_j such that the model function "best" fits the data. In linear least squares, linearity is meant to be with respect to parameters \beta_j, so

f(x, \boldsymbol \beta) = \sum_{j=1}^{n} \beta_j \phi_j(x).

Here, the functions \phi_j may be nonlinear with respect to the variable x.

Ideally, the model function fits the data exactly, so

y_i = f(x_i, \boldsymbol \beta)

for all i=1, 2, \dots, m. This is usually not possible in practice, as there are more data points than there are parameters to be determined. The approach chosen then is to find the minimal possible value of the sum of squares of the residuals

r_i(\boldsymbol \beta)= y_i - f(x_i, \boldsymbol \beta),\  (i=1, 2, \dots, m)

so to minimize the function

S(\boldsymbol \beta)=\sum_{i=1}^{m}r_i^2(\boldsymbol \beta).

After substituting for r_i and then for f, this minimization problem becomes the quadratic minimization problem above with


and the best fit can be found by solving the normal equations.

Further discussion[edit]

The numerical methods for linear least squares are important because linear regression models are among the most important types of model, both as formal statistical models and for exploration of data-sets. The majority of statistical computer packages contain facilities for regression analysis that make use of linear least squares computations. Hence it is appropriate that considerable effort has been devoted to the task of ensuring that these computations are undertaken efficiently and with due regard to round-off error.

Individual statistical analyses are seldom undertaken in isolation, but rather are part of a sequence of investigatory steps. Some of the topics involved in considering numerical methods for linear least squares relate to this point. Thus important topics can be

  • Computations where a number of similar, and often nested, models are considered for the same data-set. That is, where models with the same dependent variable but different sets of independent variables are to be considered, for essentially the same set of data-points.
  • Computations for analyses that occur in a sequence, as the number of data-points increases.
  • Special considerations for very extensive data-sets.

Fitting of linear models by least squares often, but not always, arise in the context of statistical analysis. It can therefore be important that considerations of computation efficiency for such problems extend to all of the auxiliary quantities required for such analyses, and are not restricted to the formal solution of the linear least squares problem.

Rounding errors[edit]

Matrix calculations, like any other, are affected by rounding errors. An early summary of these effects, regarding the choice of computation methods for matrix inversion, was provided by Wilkinson.[13]

See also[edit]


  1. ^ Lawson, C. L.; Hanson, R. J. (1974). Solving Least Squares Problems. Englewood Cliffs, NJ: Prentice-Hall. ISBN 0-13-822585-0. 
  2. ^ Margenau, Henry; Murphy, George Moseley (1956). The Mathematics of Physics and Chemistry. Princeton: Van Nostrand. 
  3. ^ a b Gans, Peter (1992). Data fitting in the Chemical Sciences. New York: Wiley. ISBN 0-471-93412-7. 
  4. ^ Deming, W. E. (1943). Statistical adjustment of Data. New York: Wiley. 
  5. ^ This implies that the observations are uncorrelated. If the observations are correlated, the expression \textstyle S=\sum_k \sum_j r_k W_{kj} r_j\, applies. In this case the weight matrix should ideally be equal to the inverse of the variance-covariance matrix of the observations.
  6. ^ Strutz, T. (2016). Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond). Springer Vieweg. ISBN 978-3-658-11455-8. , chapter 3
  7. ^ Mandel, John (1964). The Statistical Analysis of Experimental Data. New York: Interscience. 
  8. ^ Mardia, K. V.; Kent, J. T.; Bibby, J. M. (1979). Multivariate analysis. New York: Academic Press. ISBN 0-12-471250-9. 
  9. ^ Hamilton, W. C. (1964). Statistics in Physical Science. New York: Ronald Press. 
  10. ^ Spiegel, Murray R. (1975). Schaum's outline of theory and problems of probability and statistics. New York: McGraw-Hill. ISBN 0-585-26739-1. 
  11. ^ Acton, F. S. (1959). Analysis of Straight-Line Data. New York: Wiley. 
  12. ^ Guest, P. G. (1961). Numerical Methods of Curve Fitting. Cambridge: Cambridge University Press. [page needed]
  13. ^ Wilkinson, J.H. (1963) "Chapter 3: Matrix Computations", Rounding Errors in Algebraic Processes, London: Her Majesty's Stationery Office (National Physical Laboratory, Notes in Applied Science, No.32)

Further reading[edit]

  • Bevington, Philip R; Robinson, Keith D (2003). Data Reduction and Error Analysis for the Physical Sciences. McGraw Hill. ISBN 0-07-247227-8. 
  • Barlow, Jesse L. (1993), "Chapter 9: Numerical aspects of Solving Linear Least Squares Problems", in Rao, C.R., Computational Statistics, Handbook of Statistics 9, North-Holland, ISBN 0-444-88096-8 
  • Björck, Åke (1996). Numerical methods for least squares problems. Philadelphia: SIAM. ISBN 0-89871-360-9. 
  • Goodall, Colin R. (1993), "Chapter 13: Computation using the QR decomposition", in Rao, C.R., Computational Statistics, Handbook of Statistics 9, North-Holland, ISBN 0-444-88096-8 
  • National Physical Laboratory (1961), "Chapter 1: Linear Equations and Matrices: Direct Methods", Modern Computing Methods, Notes on Applied Science 16 (2nd ed.), Her Majesty's Stationery Office 
  • National Physical Laboratory (1961), "Chapter 2: Linear Equations and Matrices: Direct Methods on Automatic Computers", Modern Computing Methods, Notes on Applied Science 16 (2nd ed.), Her Majesty's Stationery Office 

External links[edit]