In statistics, the Gauss–Markov theorem (or simply Gauss theorem for some authors) states that the ordinary least squares (OLS) estimator has the lowest sampling variance (variance of the estimator across samples) within the class of linear unbiased estimators, if the errors in the linear regression model are uncorrelated, have equal variances and expectation value of zero. The errors do not need to be normal, nor do they need to be independent and identically distributed (only uncorrelated with mean zero and homoscedastic with finite variance). The requirement that the estimator be unbiased cannot be dropped, since biased estimators exist with lower variance. See, for example, the James–Stein estimator (which also drops linearity), ridge regression, or simply any degenerate estimator. The theorem was named after Carl Friedrich Gauss and Andrey Markov, although Gauss' work significantly predates Markov's. But while Gauss derived the result under the assumption of independence and normality, Markov reduced the assumptions to the form stated above. A further generalization to non-spherical errors was given by Alexander Aitken.
Scalar case statement Suppose we are given two random variables X , Y {\displaystyle X,Y} and that we want to find the best linear estimator of Y {\displaystyle Y} given X {\displaystyle X} , using the best linear estimator
Y ^ = α X + μ , {\displaystyle {\hat {Y}}=\alpha X+\mu ,} where the parameters α {\displaystyle \alpha } and μ {\displaystyle \mu } are both real numbers. Such an estimator Y ^ {\displaystyle {\hat {Y}}} would have the same mean and standard deviation as Y {\displaystyle Y} , that is, μ Y ^ = μ Y , σ Y ^ = σ Y {\displaystyle \mu _{\hat {Y}}=\mu _{Y},\sigma _{\hat {Y}}=\sigma _{Y}} . Therefore, if the vector X {\displaystyle X} has respective mean and standard deviation μ x , σ x {\displaystyle \mu _{x},\sigma _{x}} , the best linear estimator would be
Y ^ = σ y ( X − μ x ) σ x + μ y , {\displaystyle {\hat {Y}}=\sigma _{y}{\frac {(X-\mu _{x})}{\sigma _{x}}}+\mu _{y},} since Y ^ {\displaystyle {\hat {Y}}} has the same mean and standard deviation as Y {\displaystyle Y} .
Statement Suppose we have, in matrix notation, the linear relationship
y = X β + ε , ( y , ε ∈ R n , β ∈ R K and X ∈ R n × K ) , {\displaystyle y=X\beta +\varepsilon ,\quad (y,\varepsilon \in \mathbb {R} ^{n},\beta \in \mathbb {R} ^{K}{\text{ and }}X\in \mathbb {R} ^{n\times K}),}
where β j {\displaystyle \beta _{j}} are non-random but unobservable parameters, X i j {\displaystyle X_{ij}} are non-random and observable (called the "explanatory variables"), ε i {\displaystyle \varepsilon _{i}} are random, and so y i {\displaystyle y_{i}} are random. The random variables ε i {\displaystyle \varepsilon _{i}} are called the "disturbance", "noise" or simply "error" (will be contrasted with "residual" later in the article; see errors and residuals in statistics). Note that to include a constant in the model above, one can choose to introduce the constant as a variable β K + 1 {\displaystyle \beta _{K+1}} with a newly introduced last column of X {\displaystyle X} being unity i.e., X i ( K + 1 ) = 1 {\displaystyle X_{i(K+1)}=1} for all i {\displaystyle i} . Note that though y i {\displaystyle y_{i}} , as sample responses, are observable, the following statements and arguments including assumptions, proofs and the others assume under the only condition of knowing X i j , {\displaystyle X_{ij},} but not y i . {\displaystyle y_{i}.}
The Gauss–Markov assumptions concern the set of error random variables ε i {\displaystyle \varepsilon _{i}} :
They have mean zero: E [ ε i ] = 0. {\displaystyle \operatorname {E} [\varepsilon _{i}]=0.}
They are homoscedastic, that is, all have the same finite variance: Var ( ε i ) = σ 2 < ∞ {\displaystyle \operatorname {Var} (\varepsilon _{i})=\sigma ^{2}<\infty } for all i {\displaystyle i} . Distinct error terms are uncorrelated: Cov ( ε i , ε j ) = 0 , ∀ i ≠ j . {\displaystyle \operatorname {Cov} (\varepsilon _{i},\varepsilon _{j})=0,\forall i\neq j.}
A linear estimator of β j {\displaystyle \beta _{j}} is a linear combination
β ^ j = c 1 j y 1 + ⋯ + c n j y n , {\displaystyle {\widehat {\beta }}_{j}=c_{1j}y_{1}+\cdots +c_{nj}y_{n},}
in which the coefficients c i j {\displaystyle c_{ij}} are not allowed to depend on the underlying coefficients β j {\displaystyle \beta _{j}} , since those are not observable, but are allowed to depend on the values X i j {\displaystyle X_{ij}} , since these data are observable. (The dependence of the coefficients on each X i j {\displaystyle X_{ij}} is typically nonlinear; the estimator is linear in each y i {\displaystyle y_{i}} and hence in each random ε , {\displaystyle \varepsilon ,} which is why this is "linear" regression.) The estimator is said to be unbiased if and only if
E [ β ^ j ] = β j {\displaystyle \operatorname {E} [{\hat {\beta }}_{j}]=\beta _{j}}
regardless of the values of X i j {\displaystyle X_{ij}} . Now, let ∑ j = 1 K λ j β j {\textstyle \sum _{j=1}^{K}\lambda _{j}\beta _{j}} be some linear combination of the coefficients. Then the mean squared error of the corresponding estimation is
E [ ( ∑ j = 1 K λ j ( β ^ j − β j ) ) 2 ] , {\displaystyle \operatorname {E} \left[\left(\sum _{j=1}^{K}\lambda _{j}({\hat {\beta }}_{j}-\beta _{j})\right)^{2}\right],}
in other words, it is the expectation of the square of the weighted sum (across parameters) of the differences between the estimators and the corresponding parameters to be estimated. (Since we are considering the case in which all the parameter estimates are unbiased, this mean squared error is the same as the variance of the linear combination.) The best linear unbiased estimator (BLUE) of the vector β {\displaystyle \beta } of parameters β j {\displaystyle \beta _{j}} is one with the smallest mean squared error for every vector λ {\displaystyle \lambda } of linear combination parameters. This is equivalent to the condition that
Var ( β ~ ) − Var ( β ^ ) {\displaystyle \operatorname {Var} ({\tilde {\beta }})-\operatorname {Var} ({\hat {\beta }})}
is a positive semi-definite matrix for every other linear unbiased estimator β ~ {\displaystyle {\widetilde {\beta }}} . The ordinary least squares estimator (OLS) is the function
β ^ = ( X T X ) − 1 X T y {\displaystyle {\widehat {\beta }}=(X^{\mathsf {T}}X)^{-1}X^{\mathsf {T}}y}
of y {\displaystyle y} and X {\displaystyle X} (where X T {\displaystyle X^{\mathsf {T}}} denotes the transpose of X {\displaystyle X} ) that minimizes the sum of squares of residuals (misprediction amounts):
∑ i = 1 n ( y i − y ^ i ) 2 = ∑ i = 1 n ( y i − ∑ j = 1 K β ^ j X i j ) 2 . {\displaystyle \sum _{i=1}^{n}(y_{i}-{\hat {y}}_{i})^{2}=\sum _{i=1}^{n}\left(y_{i}-\sum _{j=1}^{K}{\hat {\beta }}_{j}X_{ij}\right)^{2}.}
The theorem now states that the OLS estimator is a best linear unbiased estimator (BLUE). The main idea of the proof is that the least-squares estimator is uncorrelated with every linear unbiased estimator of zero, i.e., with every linear combination a 1 y 1 + ⋯ + a n y n {\displaystyle a_{1}y_{1}+\cdots +a_{n}y_{n}} whose coefficients do not depend upon the unobservable β {\displaystyle \beta } but whose expected value is always zero.
Remark Proof that the OLS indeed minimizes the sum of squares of residuals may proceed as follows with a calculation of the Hessian matrix and showing that it is positive definite. The MSE function we want to minimize is
f ( β 0 , β 1 , … , β p ) = ∑ i = 1 n ( y i − β 0 − β 1 x i 1 − ⋯ − β p x i p ) 2 {\displaystyle f(\beta _{0},\beta _{1},\dots ,\beta _{p})=\sum _{i=1}^{n}(y_{i}-\beta _{0}-\beta _{1}x_{i1}-\dots -\beta _{p}x_{ip})^{2}}
for a multiple regression model with p variables. The first derivative is
d d β f = − 2 X T ( y − X β ) = − 2 [ ∑ i = 1 n ( y i − ⋯ − β p x i p ) ∑ i = 1 n x i 1 ( y i − ⋯ − β p x i p ) ⋮ ∑ i = 1 n x i p ( y i − ⋯ − β p x i p ) ] = 0 p + 1 , {\displaystyle {\begin{aligned}{\frac {d}{d{\boldsymbol {\beta }}}}f&=-2X^{\operatorname {T} }\left(\mathbf {y} -X{\boldsymbol {\beta }}\right)\\&=-2{\begin{bmatrix}\sum _{i=1}^{n}(y_{i}-\dots -\beta _{p}x_{ip})\\\sum _{i=1}^{n}x_{i1}(y_{i}-\dots -\beta _{p}x_{ip})\\\vdots \\\sum _{i=1}^{n}x_{ip}(y_{i}-\dots -\beta _{p}x_{ip})\end{bmatrix}}\\&=\mathbf {0} _{p+1},\end{aligned}}}
where X T {\displaystyle X^{\operatorname {T} }} is the design matrix
X = [ 1 x 11 ⋯ x 1 p 1 x 21 ⋯ x 2 p ⋮ 1 x n 1 ⋯ x n p ] ∈ R n × ( p + 1 ) ; n ≥ p + 1 {\displaystyle X={\begin{bmatrix}1&x_{11}&\cdots &x_{1p}\\1&x_{21}&\cdots &x_{2p}\\&&\vdots \\1&x_{n1}&\cdots &x_{np}\end{bmatrix}}\in \mathbb {R} ^{n\times (p+1)};\qquad n\geq p+1}
The Hessian matrix of second derivatives is
H = 2 [ n ∑ i = 1 n x i 1 ⋯ ∑ i = 1 n x i p ∑ i = 1 n x i 1 ∑ i = 1 n x i 1 2 ⋯ ∑ i = 1 n x i 1 x i p ⋮ ⋮ ⋱ ⋮ ∑ i = 1 n x i p ∑ i = 1 n x i p x i 1 ⋯ ∑ i = 1 n x i p 2 ] = 2 X T X {\displaystyle {\mathcal {H}}=2{\begin{bmatrix}n&\sum _{i=1}^{n}x_{i1}&\cdots &\sum _{i=1}^{n}x_{ip}\\\sum _{i=1}^{n}x_{i1}&\sum _{i=1}^{n}x_{i1}^{2}&\cdots &\sum _{i=1}^{n}x_{i1}x_{ip}\\\vdots &\vdots &\ddots &\vdots \\\sum _{i=1}^{n}x_{ip}&\sum _{i=1}^{n}x_{ip}x_{i1}&\cdots &\sum _{i=1}^{n}x_{ip}^{2}\end{bmatrix}}=2X^{\operatorname {T} }X}
Assuming the columns of X {\displaystyle X} are linearly independent so that X T X {\displaystyle X^{\operatorname {T} }X} is invertible, let X = [ v 1 v 2 ⋯ v p + 1 ] {\displaystyle X={\begin{bmatrix}\mathbf {v_{1}} &\mathbf {v_{2}} &\cdots &\mathbf {v} _{p+1}\end{bmatrix}}} , then
k 1 v 1 + ⋯ + k p + 1 v p + 1 = 0 ⟺ k 1 = ⋯ = k p + 1 = 0 {\displaystyle k_{1}\mathbf {v_{1}} +\dots +k_{p+1}\mathbf {v} _{p+1}=\mathbf {0} \iff k_{1}=\dots =k_{p+1}=0}
Now let k = ( k 1 , … , k p + 1 ) T ∈ R ( p + 1 ) × 1 {\displaystyle \mathbf {k} =(k_{1},\dots ,k_{p+1})^{T}\in \mathbb {R} ^{(p+1)\times 1}} be an eigenvector of H {\displaystyle {\mathcal {H}}} .
k ≠ 0 ⟹ ( k 1 v 1 + ⋯ + k p + 1 v p + 1 ) 2 > 0 {\displaystyle \mathbf {k} \neq \mathbf {0} \implies \left(k_{1}\mathbf {v_{1}} +\dots +k_{p+1}\mathbf {v} _{p+1}\right)^{2}>0}
In terms of vector multiplication, this means
[ k 1 ⋯ k p + 1 ] [ v 1 ⋮ v p + 1 ] [ v 1 ⋯ v p + 1 ] [ k 1 ⋮ k p + 1 ] = k T H k = λ k T k > 0 {\displaystyle {\begin{bmatrix}k_{1}&\cdots &k_{p+1}\end{bmatrix}}{\begin{bmatrix}\mathbf {v_{1}} \\\vdots \\\mathbf {v} _{p+1}\end{bmatrix}}{\begin{bmatrix}\mathbf {v_{1}} &\cdots &\mathbf {v} _{p+1}\end{bmatrix}}{\begin{bmatrix}k_{1}\\\vdots \\k_{p+1}\end{bmatrix}}=\mathbf {k} ^{\operatorname {T} }{\mathcal {H}}\mathbf {k} =\lambda \mathbf {k} ^{\operatorname {T} }\mathbf {k} >0}
where λ {\displaystyle \lambda } is the eigenvalue corresponding to k {\displaystyle \mathbf {k} } . Moreover,
k T k = ∑ i = 1 p + 1