Linear least squares

Other formulations include:

• Iteratively reweighted least squares (IRLS) is used when heteroscedasticity, or correlations, or both are present among the error terms of the model, but where little is known about the covariance structure of the errors independently of the data.[2] In the first iteration, OLS, or GLS with a provisional covariance structure is carried out, and the residuals are obtained from the fit. Based on the residuals, an improved estimate of the covariance structure of the errors can usually be obtained. A subsequent GLS iteration is then performed using this estimate of the error structure to define the weights. The process can be iterated to convergence, but in many cases, only one iteration is sufficient to achieve an efficient estimate of β.[3][4]
• Instrumental variables regression (IV) can be performed when the regressors are correlated with the errors. In this case, we need the existence of some auxiliary instrumental variables z_i such that E[z_iε_i] = 0. If Z is the matrix of instruments, then the estimator can be given in closed form as
  $${\displaystyle {\hat {\boldsymbol {\beta }}}=(\mathbf {X} ^{\mathsf {T}}\mathbf {Z} (\mathbf {Z} ^{\mathsf {T}}\mathbf {Z} )^{-1}\mathbf {Z} ^{\mathsf {T}}\mathbf {X} )^{-1}\mathbf {X} ^{\mathsf {T}}\mathbf {Z} (\mathbf {Z} ^{\mathsf {T}}\mathbf {Z} )^{-1}\mathbf {Z} ^{\mathsf {T}}\mathbf {y} .}$$

  
  Optimal instruments regression is an extension of classical IV regression to the situation where E[ε_i | z_i] = 0.
• Total least squares (TLS)[5] is an approach to least squares estimation of the linear regression model that treats the covariates and response variable in a more geometrically symmetric manner than OLS. It is one approach to handling the "errors in variables" problem, and is also sometimes used even when the covariates are assumed to be error-free.

• Linear Template Fit (LTF)[6] combines a linear regression with (generalized) least squares in order to determine the best estimator. The Linear Template Fit addresses the frequent issue, when the residuals cannot be expressed analytically or are too time consuming to be evaluate repeatedly, as it is often the case in iterative minimization algorithms. In the Linear Template Fit, the residuals are estimated from the random variables and from a linear approximation of the underlying true model, while the true model needs to be provided for at least ${\displaystyle n+1}$ (were ${\displaystyle n}$ is the number of estimators) distinct reference values β. The true distribution is then approximated by a linear regression, and the best estimators are obtained in closed form as
  $${\displaystyle {\hat {\boldsymbol {\beta }}}=((\mathbf {Y{\tilde {M}}} )^{\mathsf {T}}{\boldsymbol {\Omega }}^{-1}\mathbf {Y{\tilde {M}}} )^{-1}(\mathbf {Y{\tilde {M}}} )^{\mathsf {T}}{\boldsymbol {\Omega }}^{-1}(\mathbf {d} -\mathbf {Y{\bar {m}})} ,}$$

  
  where ${\displaystyle \mathbf {Y} }$ denotes the template matrix with the values of the known or previously determined model for any of the reference values β, ${\displaystyle \mathbf {d} }$ are the random variables (e.g. a measurement), and the matrix ${\displaystyle \mathbf {\tilde {M}} }$ and the vector ${\displaystyle \mathbf {\tilde {m}} }$ are calculated from the values of β. The LTF can also be expressed for Log-normal distribution distributed random variables. A generalization of the LTF is the Quadratic Template Fit, which assumes a second order regression of the model, requires predictions for at least ${\displaystyle n^{2}+2n}$ distinct values β, and it finds the best estimator using Newton's method.

• Percentage least squares focuses on reducing percentage errors, which is useful in the field of forecasting or time series analysis. It is also useful in situations where the dependent variable has a wide range without constant variance, as here the larger residuals at the upper end of the range would dominate if OLS were used. When the percentage or relative error is normally distributed, least squares percentage regression provides maximum likelihood estimates. Percentage regression is linked to a multiplicative error model, whereas OLS is linked to models containing an additive error term.[7]

• Constrained least squares, indicates a linear least squares problem with additional constraints on the solution.

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.[11][12]

In some cases the (weighted) normal equations matrix X^TX 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. In these cases, the least squares estimate amplifies the measurement noise and may be grossly inaccurate. 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 ${\displaystyle \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, ${\displaystyle \|\mathbf {y} -\mathbf {X} {\hat {\boldsymbol {\beta }}}\|}$ is minimized, whereas in some cases one is truly interested in obtaining small error in the parameter ${\displaystyle \mathbf {\hat {\boldsymbol {\beta }}} }$, e.g., a small value of ${\displaystyle \|{\boldsymbol {\beta }}-{\hat {\boldsymbol {\beta }}}\|}$. However, since the true parameter ${\displaystyle {\boldsymbol {\beta }}}$ is necessarily unknown, this quantity cannot be directly minimized. If a prior probability on ${\displaystyle {\hat {\boldsymbol {\beta }}}}$ is known, then a Bayes estimator can be used to minimize the mean squared error, ${\displaystyle 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.

The primary application of linear least squares is in data fitting. Given a set of m data points ${\displaystyle y_{1},y_{2},\dots ,y_{m},}$ consisting of experimentally measured values taken at m values ${\displaystyle x_{1},x_{2},\dots ,x_{m}}$ of an independent variable (${\displaystyle x_{i}}$ may be scalar or vector quantities), and given a model function ${\displaystyle y=f(x,{\boldsymbol {\beta }}),}$ with ${\displaystyle {\boldsymbol {\beta }}=(\beta _{1},\beta _{2},\dots ,\beta _{n}),}$ it is desired to find the parameters ${\displaystyle \beta _{j}}$ such that the model function "best" fits the data. In linear least squares, linearity is meant to be with respect to parameters ${\displaystyle \beta _{j},}$ so

$${\displaystyle f(x,{\boldsymbol {\beta }})=\sum _{j=1}^{n}\beta _{j}\varphi _{j}(x).}$$

Here, the functions ${\displaystyle \varphi _{j}}$ may be nonlinear with respect to the variable x.

Ideally, the model function fits the data exactly, so

$${\displaystyle y_{i}=f(x_{i},{\boldsymbol {\beta }})}$$

for all ${\displaystyle 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

$${\displaystyle r_{i}({\boldsymbol {\beta }})=y_{i}-f(x_{i},{\boldsymbol {\beta }}),\ (i=1,2,\dots ,m)}$$

so to minimize the function

$${\displaystyle S({\boldsymbol {\beta }})=\sum _{i=1}^{m}r_{i}^{2}({\boldsymbol {\beta }}).}$$

After substituting for ${\displaystyle r_{i}}$ and then for ${\displaystyle f}$, this minimization problem becomes the quadratic minimization problem above with

$${\displaystyle X_{ij}=\varphi _{j}(x_{i}),}$$

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

One of the simplest possible relationships between ${\displaystyle x}$ and ${\displaystyle y}$ is a line ${\displaystyle y=\beta _{1}+\beta _{2}x}$. The intercept ${\displaystyle \beta _{1}}$ and the slope ${\displaystyle \beta _{2}}$ are initially unknown. The researcher would like to find values of ${\displaystyle \beta _{1}}$ and ${\displaystyle \beta _{2}}$ that cause the line to pass through the four data points. In other words, the researcher would like to solve the system of linear equations

$${\displaystyle {\begin{alignedat}{3}\beta _{1}+1\beta _{2}&&\;=\;&&6,&\\\beta _{1}+2\beta _{2}&&\;=\;&&5,&\\\beta _{1}+3\beta _{2}&&\;=\;&&7,&\\\beta _{1}+4\beta _{2}&&\;=\;&&10.&\\\end{alignedat}}}$$

With four equations in two unknowns, this system is overdetermined. There is no exact solution. To consider approximate solutions, one introduces residuals ${\displaystyle r_{1}}$, ${\displaystyle r_{2}}$, ${\displaystyle r_{3}}$, ${\displaystyle r_{4}}$ into the equations:

$${\displaystyle {\begin{alignedat}{3}\beta _{1}+1\beta _{2}+r_{1}&&\;=\;&&6,&\\\beta _{1}+2\beta _{2}+r_{2}&&\;=\;&&5,&\\\beta _{1}+3\beta _{2}+r_{3}&&\;=\;&&7,&\\\beta _{1}+4\beta _{2}+r_{4}&&\;=\;&&10.&\\\end{alignedat}}}$$

The ${\displaystyle i}$th residual ${\displaystyle r_{i}}$ is the misfit between the ${\displaystyle i}$th observation ${\displaystyle y_{i}}$ and the ${\displaystyle i}$th prediction ${\displaystyle \beta _{1}+\beta _{2}x_{i}}$:

$${\displaystyle {\begin{alignedat}{3}r_{1}&&\;=\;&&6-(\beta _{1}+1\beta _{2}),&\\r_{2}&&\;=\;&&5-(\beta _{1}+2\beta _{2}),&\\r_{3}&&\;=\;&&7-(\beta _{1}+3\beta _{2}),&\\r_{4}&&\;=\;&&10-(\beta _{1}+4\beta _{2}).&\\\end{alignedat}}}$$

Among all approximate solutions, the researcher would like to find the one that is "best" in some sense.

In least squares, one focuses on the sum ${\displaystyle S}$ of the squared residuals:

$${\displaystyle {\begin{aligned}S(\beta _{1},\beta _{2})&=r_{1}^{2}+r_{2}^{2}+r_{3}^{2}+r_{4}^{2}\\[6pt]&=[6-(\beta _{1}+1\beta _{2})]^{2}+[5-(\beta _{1}+2\beta _{2})]^{2}+[7-(\beta _{1}+3\beta _{2})]^{2}+[10-(\beta _{1}+4\beta _{2})]^{2}\\[6pt]&=4\beta _{1}^{2}+30\beta _{2}^{2}+20\beta _{1}\beta _{2}-56\beta _{1}-154\beta _{2}+210.\\[6pt]\end{aligned}}}$$

The best solution is defined to be the one that minimizes ${\displaystyle S}$ with respect to ${\displaystyle \beta _{1}}$ and ${\displaystyle \beta _{2}}$. The minimum can be calculated by setting the partial derivatives of ${\displaystyle S}$ to zero:

$${\displaystyle 0={\frac {\partial S}{\partial \beta _{1}}}=8\beta _{1}+20\beta _{2}-56,}$$

$${\displaystyle 0={\frac {\partial S}{\partial \beta _{2}}}=20\beta _{1}+60\beta _{2}-154.}$$

These normal equations constitute a system of two linear equations in two unknowns. The solution is ${\displaystyle \beta _{1}=3.5}$ and ${\displaystyle \beta _{2}=1.4}$, and the best-fit line is therefore ${\displaystyle y=3.5+1.4x}$. The residuals are ${\displaystyle 1.1,}$ ${\displaystyle -1.3,}$ ${\displaystyle -0.7,}$ and ${\displaystyle 0.9}$ (see the diagram on the right). The minimum value of the sum of squared residuals is

$${\displaystyle S(3.5,1.4)=1.1^{2}+(-1.3)^{2}+(-0.7)^{2}+0.9^{2}=4.2.}$$

This calculation can be expressed in matrix notation as follows. The original system of equations is ${\displaystyle \mathbf {y} =\mathbf {X} \mathbf {\beta } }$, where

$${\displaystyle \mathbf {y} =\left[{\begin{array}{c}6\\5\\7\\10\end{array}}\right],\;\;\;\;\mathbf {X} =\left[{\begin{array}{cc}1&1\\1&2\\1&3\\1&4\end{array}}\right],\;\;\;\;\mathbf {\beta } =\left[{\begin{array}{c}\beta _{1}\\\beta _{2}\end{array}}\right].}$$

Intuitively,

$${\displaystyle \mathbf {y} =\mathbf {X} \mathbf {\beta } \;\;\;\;\Rightarrow \;\;\;\;\mathbf {X} ^{\top }\mathbf {y} =\mathbf {X} ^{\top }\mathbf {X} \mathbf {\beta } \;\;\;\;\Rightarrow \;\;\;\;\mathbf {\beta } =\left(\mathbf {X} ^{\top }\mathbf {X} \right)^{-1}\mathbf {X} ^{\top }\mathbf {y} =\left[{\begin{array}{c}3.5\\1.4\end{array}}\right].}$$

More rigorously, if ${\displaystyle \mathbf {X} ^{\top }\mathbf {X} }$ is invertible, then the matrix ${\displaystyle \mathbf {X} \left(\mathbf {X} ^{\top }\mathbf {X} \right)^{-1}\mathbf {X} ^{\top }}$ represents orthogonal projection onto the column space of ${\displaystyle \mathbf {X} }$. Therefore, among all vectors of the form ${\displaystyle \mathbf {X} \mathbf {\beta } }$, the one closest to ${\displaystyle \mathbf {y} }$ is ${\displaystyle \mathbf {X} \left(\mathbf {X} ^{\top }\mathbf {X} \right)^{-1}\mathbf {X} ^{\top }\mathbf {y} }$. Setting

$${\displaystyle \mathbf {X} \left(\mathbf {X} ^{\top }\mathbf {X} \right)^{-1}\mathbf {X} ^{\top }\mathbf {y} =\mathbf {X} \mathbf {\beta } ,}$$

it is evident that ${\displaystyle \mathbf {\beta } =\left(\mathbf {X} ^{\top }\mathbf {X} \right)^{-1}\mathbf {X} ^{\top }\mathbf {y} }$ is a solution.

The result of fitting a quadratic function ${\displaystyle y=\beta _{1}+\beta _{2}x+\beta _{3}x^{2}\,}$ (in blue) through a set of data points ${\displaystyle (x_{i},y_{i})}$ (in red). In linear least squares the function need not be linear in the argument ${\displaystyle x,}$ but only in the parameters ${\displaystyle \beta _{j}}$ that are determined to give the best fit.

Suppose that the hypothetical researcher wishes to fit a parabola of the form ${\displaystyle y=\beta _{1}x^{2}}$. Importantly, this model is still linear in the unknown parameters (now just ${\displaystyle \beta _{1}}$), so linear least squares still applies. The system of equations incorporating residuals is

$${\displaystyle {\begin{alignedat}{2}6&&\;=\beta _{1}(1)^{2}+r_{1}\\5&&\;=\beta _{1}(2)^{2}+r_{2}\\7&&\;=\beta _{1}(3)^{2}+r_{3}\\10&&\;=\beta _{1}(4)^{2}+r_{4}\\\end{alignedat}}}$$

The sum of squared residuals is

$${\displaystyle S(\beta _{1})=(6-\beta _{1})^{2}+(5-4\beta _{1})^{2}+(7-9\beta _{1})^{2}+(10-16\beta _{1})^{2}.}$$

There is just one partial derivative to set to 0:

$${\displaystyle 0={\frac {\partial S}{\partial \beta _{1}}}=708\beta _{1}-498.}$$

The solution is ${\displaystyle \beta _{1}=0.703}$, and the fit model is ${\displaystyle y=0.703x^{2}}$.

In matrix notation, the equations without residuals are again ${\displaystyle \mathbf {y} =\mathbf {X} \mathbf {\beta } }$, where now

$${\displaystyle \mathbf {y} =\left[{\begin{array}{c}6\\5\\7\\10\end{array}}\right],\;\;\;\;\mathbf {X} =\left[{\begin{array}{c}1\\4\\9\\16\end{array}}\right],\;\;\;\;\mathbf {\beta } =\left[{\begin{array}{c}\beta _{1}\end{array}}\right].}$$

By the same logic as above, the solution is

$${\displaystyle \mathbf {\beta } =\left(\mathbf {X} ^{\top }\mathbf {X} \right)^{-1}\mathbf {X} ^{\top }\mathbf {y} =\left[{\begin{array}{c}0.703\end{array}}\right].}$$

The figure shows an extension to fitting the three parameter parabola using a design matrix ${\displaystyle \mathbf {X} }$ with three columns (one for ${\displaystyle x^{0}}$, ${\displaystyle x^{1}}$, and ${\displaystyle x^{2}}$), and one row for each of the red data points.

More generally, one can have ${\displaystyle n}$ regressors ${\displaystyle x_{j}}$, and a linear model

$${\displaystyle y=\beta _{0}+\sum _{j=1}^{n}\beta _{j}x_{j}.}$$