Linear Model Theory: Gauss-Markov, the Hat Matrix, and BLUE

Intuition Beginner

You have a stack of measurements, and you suspect that one number (the response) is built up out of a few other numbers (the predictors) plus some random noise. The linear model is the simplest honest way to write that down: the response is a weighted sum of the predictors, and whatever the weighted sum misses is chalked up to error. The whole game is to recover the weights from noisy data and to say how much you should trust them.

The recipe for the weights is least squares: pick the weights that make the total squared miss as small as possible. Squaring the misses means that being off by 4 on one point costs the same as being off by 2 on each of two points only when the misses are equal; large misses are punished much harder than small ones. This single choice turns out to have a clean geometry and a clean statistical pay-off.

The geometry is a shadow. Imagine your response values as the coordinates of one point floating in a high-dimensional space, and imagine all the response patterns your model could possibly produce as a flat sheet (a plane) inside that space. Least squares drops a perpendicular from your point straight down onto the sheet. The foot of that perpendicular is the model's best guess; the leftover vertical drop is the part the model cannot reach, the residual.

Why squared error and not, say, absolute error? Because under a mild and very common set of assumptions about the noise, least squares is provably the best you can do within a fair contest. The contest restricts attention to weight-recipes that are unbiased (right on average) and that combine the data in a straight linear way. Among all such recipes, least squares has the smallest spread. That is the promise this unit makes precise, and it is the reason linear models are the bedrock of statistics.

Visual Beginner

Object	Plain meaning	Lives where
Response vector	the column of measured outcomes	a point in $n$-dimensional space
Model plane	every outcome pattern the weights could produce	a flat sheet through the origin
Fitted values	the model's best reconstruction of the outcomes	the foot of the perpendicular
Residual	the part the model cannot reach	the perpendicular drop itself
Hat matrix	the machine that drops the perpendicular	a square table of numbers
Leverage	how much one point pulls its own fit	one number per data point

The single picture to carry forward: the fitted values are a shadow of the data cast straight down onto the sheet of reachable patterns, and the residual is the gap between the data and its shadow. Everything else in this unit is bookkeeping about that shadow.

Worked example Beginner

Take four data points with one predictor and an intercept. The predictor values are $x=0,1,2,3$ and the responses are $y=1,3,4,6$. We fit $\hat{y}=a+bx$.

First the averages: $\bar{x}=(0+1+2+3)/4=1.5$ and $\bar{y}=(1+3+4+6)/4=3.5$.

Now the slope. The numerator is the sum of $(x_i-\bar{x})(y_i-\bar{y})$:

$(0-1.5)(1-3.5)+(1-1.5)(3-3.5)+(2-1.5)(4-3.5)+(3-1.5)(6-3.5)$ $=(-1.5)(-2.5)+(-0.5)(-0.5)+(0.5)(0.5)+(1.5)(2.5)=3.75+0.25+0.25+3.75=8$.

The denominator is the sum of $(x_i-\bar{x}{)}^2=2.25+0.25+0.25+2.25=5$.

So the slope is $b=8/5=1.6$, and the intercept is $a=\bar{y}-b\bar{x}=3.5-1.6\times 1.5=3.5-2.4=1.1$.

The fitted line is $\hat{y}=1.1+1.6x$. The fitted values at the four predictors are $1.1,2.7,4.3,5.9$. The residuals (observed minus fitted) are $-0.1,0.3,-0.3,0.1$.

Two checks confirm the geometry. The residuals add to $-0.1+0.3-0.3+0.1=0$: the perpendicular drop has no leftover component along the all-ones direction. And the sum of $x_i$ times the residuals is $0(-0.1)+1(0.3)+2(-0.3)+3(0.1)=0+0.3-0.6+0.3=0$: no leftover component along the predictor either. What this tells us is that the residual is genuinely perpendicular to the model plane, exactly as the shadow picture demands.

Check your understanding Beginner

Formal definition Intermediate+

Fix a design matrix $X\in {\mathbb{R}}^{n\times p}$ of full column rank $p<n$, and consider the linear model

$$y=X\beta +\epsilon ,\mathbb{E}[\epsilon ]=0,\mathrm{Cov}(\epsilon )={\sigma }^2{I}_n,$$

where $y\in {\mathbb{R}}^n$ is the observed response, $\beta \in {\mathbb{R}}^p$ is the unknown parameter, and $\epsilon$ is a random error vector. The two moment conditions — zero mean and the scalar covariance ${\sigma }^2{I}_n$ — are the second-order or Gauss-Markov assumptions: the errors are uncorrelated (off-diagonal covariances vanish) and homoscedastic (equal variances ${\sigma }^2$ down the diagonal). No distributional form is assumed yet.

The ordinary least-squares (OLS) estimator is any minimizer of the residual sum of squares $\mathrm{RSS}(\beta )=\Vert y-X\beta {\Vert }^2$. Setting the gradient to zero gives the normal equations $X^{\top }X\hat{\beta }=X^{\top }y$. Full column rank makes $X^{\top }X$ invertible, so

$$\hat{\beta }=(X^{\top }X{)}^{-1}{X}^{\top }y.$$

The numerical algorithms that solve this system stably (QR, SVD) belong to the computational treatment 43.04.01; here the object of study is the estimator's statistical behaviour, sharpening the elementary survey 26.06.01.

The fitted values are $\hat{y}=X\hat{\beta }=Hy$ where

$$H=X(X^{\top }X{)}^{-1}{X}^{\top }$$

is the hat matrix (it puts the hat on $y$). $H$ is symmetric and idempotent, $H^2=H=H^{\top }$, hence it is the orthogonal projection onto the column space $\mathrm{col}(X)$ — the same projection structure as the spectral/projection theory 01.01.13. Its complement $M=I-H$ projects onto the orthogonal complement, and the residual vector is $\hat{\epsilon }=y-\hat{y}=My$. The diagonal entries $h_{ii}$ of $H$ are the leverages; since $H$ is a projection, $0\le h_{ii}\le 1$ and ${\sum }_i{h}_{ii}=\mathrm{tr}(H)=p$.

The residual sum of squares is $\mathrm{RSS}=\Vert \hat{\epsilon }{\Vert }^2=y^{\top }My$, and the unbiased variance estimator is ${\hat{\sigma }}^2=\mathrm{RSS}/(n-p)$. The coefficient of determination is $R^2=1-\mathrm{RSS}/\mathrm{TSS}$ with $\mathrm{TSS}={\sum }_i(y_i-\bar{y}{)}^2$, the fraction of response variance the model captures.

Counterexamples to common slips

• Idempotence is not optionality. $H$ being a projection does not make $Hy=y$; that holds only when $y$ already lies in $\mathrm{col}(X)$, i.e. the residual is zero. In general $Hy$ is the nearest point of the column space.
• BLUE is not BUE. Dropping the linearity restriction breaks the optimality claim: under normal errors a nonlinear estimator cannot beat OLS, but the Gauss-Markov theorem itself never promises that — it only ranks OLS first among linear unbiased estimators.
• Leverage is not residual size. A point with $h_{ii}$ near $1$ has high leverage (an extreme predictor configuration) regardless of whether its residual is large; high leverage with a small residual still signals an influential point.
• $R^2$ does not validate a model. A high $R^2$ is compatible with a badly mis-specified mean function; it measures variance explained under the fitted model, not correctness of the model.

Key theorem with proof Intermediate+

Gauss-Markov theorem. Under the linear model $y=X\beta +\epsilon$ with $\mathbb{E}[\epsilon ]=0$ and $\mathrm{Cov}(\epsilon )={\sigma }^2{I}_n$ and $X$ of full column rank, the OLS estimator $\hat{\beta }=(X^{\top }X{)}^{-1}{X}^{\top }y$ is the best linear unbiased estimator (BLUE): for every linear unbiased estimator $\tilde{\beta }=Cy$, the difference $\mathrm{Cov}(\tilde{\beta })-\mathrm{Cov}(\hat{\beta })$ is positive semidefinite. The proof follows Hastie, Tibshirani & Friedman ^{[Hastie 2009]} §3.2.2 and Seber & Lee ^{[Seber 2003]} Ch. 3.

Proof. A linear estimator is $\tilde{\beta }=Cy$ for some fixed $C\in {\mathbb{R}}^{p\times n}$. Unbiasedness requires $\mathbb{E}[Cy]=CX\beta =\beta$ for every $\beta \in {\mathbb{R}}^p$, which forces the constraint $CX=I_p$. Write $A=(X^{\top }X{)}^{-1}{X}^{\top }$, so that $\hat{\beta }=Ay$ and $AX=I_p$, and set $D=C-A$. The unbiasedness constraint gives $DX=CX-AX=I_p-I_p=0$.

Because $\mathrm{Cov}(y)={\sigma }^2{I}_n$, the covariance of any linear estimator is $\mathrm{Cov}(Cy)={\sigma }^{2}C{C}^{\top }$. Expand with $C=A+D$:

$$C{C}^{\top }=(A+D)(A+D{)}^{\top }=A{A}^{\top }+A{D}^{\top }+D{A}^{\top }+D{D}^{\top }.$$

The cross terms vanish. Indeed $A{D}^{\top }=(X^{\top }X{)}^{-1}{X}^{\top }{D}^{\top }=(X^{\top }X{)}^{-1}(DX{)}^{\top }=0$ since $DX=0$, and $D{A}^{\top }$ is its transpose, also zero. Therefore

$$\mathrm{Cov}(\tilde{\beta })={\sigma }^{2}C{C}^{\top }={\sigma }^{2}A{A}^{\top }+{\sigma }^{2}D{D}^{\top }=\mathrm{Cov}(\hat{\beta })+{\sigma }^{2}D{D}^{\top },$$

using $A{A}^{\top }=(X^{\top }X{)}^{-1}{X}^{\top }X(X^{\top }X{)}^{-1}=(X^{\top }X{)}^{-1}$, which is exactly $\mathrm{Cov}(\hat{\beta })$. The matrix $D{D}^{\top }$ is positive semidefinite, so $\mathrm{Cov}(\tilde{\beta })-\mathrm{Cov}(\hat{\beta })={\sigma }^{2}D{D}^{\top }⪰0$. Equality in the Loewner order holds for a given direction iff the corresponding rows of $D$ vanish, and $D=0$ recovers $\tilde{\beta }=\hat{\beta }$, so OLS is the unique BLUE up to almost-sure equality. $\square$

Bridge. The proof's engine is the orthogonal decomposition $C=A+D$ with $DX=0$: the deviation $D$ from the least-squares recipe is forced into the orthogonal complement of the column space, where it can only add variance. This is exactly the projection geometry of 01.01.13 wearing statistical clothing, and it is the foundational reason least squares is optimal rather than merely convenient. The result builds toward the distribution theory below — once normal errors are added, the same $H=X(X^{\top }X{)}^{-1}{X}^{\top }$ splits the variance into independent chi-squared pieces — and it appears again in the regularized estimators of 45.06.03, where deliberately accepting bias (so the unbiasedness constraint $CX=I$ is dropped) buys a variance reduction that beats the BLUE in mean squared error. Putting these together, Gauss-Markov is the optimality benchmark that shrinkage methods are measured against; the central insight is that the theorem's hypotheses, not its conclusion, are where modern high-dimensional statistics intervenes.

Exercises Intermediate+

Lean formalization Intermediate+

Mathlib does not yet host the statistical linear model, so no module is wired in (lean_status: none). The intended statement of the Gauss-Markov theorem, once a multivariate-Gaussian-with-given-covariance API exists, would read roughly as follows.

-- Intended shape; not part of the current Babel Bible Lean build.
-- Requires: a `Cov` operator on random vectors and the Loewner order `⪯`.
variable {n p : ℕ} (X : Matrix (Fin n) (Fin p) ℝ)
variable (hX : (Xᵀ * X).det ≠ 0)  -- full column rank

-- OLS estimator as a fixed linear map applied to y
noncomputable def olsHat : Matrix (Fin p) (Fin n) ℝ :=
  (Xᵀ * X)⁻¹ * Xᵀ

-- Gauss-Markov: for any linear unbiased C (i.e. C * X = 1),
-- Cov (C y) - Cov (olsHat y) is positive semidefinite.
theorem gauss_markov
    (C : Matrix (Fin p) (Fin n) ℝ) (hC : C * X = 1)
    (σ2 : ℝ) (hσ : 0 ≤ σ2) :
    (olsHat X) * (olsHat X)ᵀ ⪯ C * Cᵀ := by
  sorry

the Mathlib gap analysis records what is missing: the random-vector covariance operator, the Loewner-order API on covariance matrices (Mathlib has Matrix.PosSemidef but not the named statistical model), and the chi-squared/F distribution layer needed for the inference half.

Advanced results Master

Assume henceforth normal errors, $\epsilon \sim N(0,{\sigma }^2{I}_n)$. The estimator $\hat{\beta }=(X^{\top }X{)}^{-1}{X}^{\top }y$ is then a linear image of a Gaussian, hence Gaussian:

$$\hat{\beta }\sim N(\beta ,{\sigma }^2(X^{\top }X{)}^{-1}).$$

The residual quadratic form is governed by the projection $M=I-H$, which is idempotent of rank $n-p$. Writing $\mathrm{RSS}={\epsilon }^{\top }M\epsilon$ and diagonalizing $M$ (its eigenvalues are $n-p$ ones and $p$ zeros), the standardized residual sum of squares is a sum of $n-p$ independent squared standard normals:

$$\frac{\mathrm{RSS}}{{\sigma }^2}=\frac{(n-p){\hat{\sigma }}^2}{{\sigma }^2}\sim {\chi }_{n-p}^2.$$

This is the distributional core of the theory: a quadratic form in a Gaussian vector with idempotent kernel of rank $r$ is exactly ${\sigma }^2{\chi }_r^2$. Crucially $\hat{\beta }$ and $\mathrm{RSS}$ are independent. The estimator depends on $y$ only through $Hy$ and the residuals only through $My$; the two projections are orthogonal ($HM=0$), and orthogonal linear images of a single Gaussian are independent. This is Cochran's theorem specialized to the regression decomposition $y=Hy+My$.

Independence is what makes the $t$-test valid. For a single coordinate, ${\hat{\beta }}_j\sim N({\beta }_j,{\sigma }^2[(X^{\top }X{)}^{-1}{]}_{jj})$. Dividing the standardized coefficient (which uses the unknown $\sigma$) by the independent $\sqrt{{\hat{\sigma }}^2/{\sigma }^2}$ converts a normal over a scaled independent chi root into Student's distribution:

$$t_j=\frac{{\hat{\beta }}_j-{\beta }_j}{\hat{\sigma }\sqrt{[(X^{\top }X{)}^{-1}{]}_{jj}}}\sim t_{n-p}.$$

For a linear hypothesis $H_0:A\beta =c$ with $A\in {\mathbb{R}}^{q\times p}$ of rank $q$, the same machinery yields the $F$-test. The numerator measures how far the constrained fit lies from the unconstrained fit, scaled by $q$; the denominator is ${\hat{\sigma }}^2$. Concretely,

$$F=\frac{(A\hat{\beta }-c{)}^{\top }[A(X^{\top }X{)}^{-1}{A}^{\top }{]}^{-1}(A\hat{\beta }-c)/q}{{\hat{\sigma }}^2}\sim F_{q,n-p}\text{under }{H}_0,$$

because the numerator quadratic form is ${\sigma }^2{\chi }_q^2$, independent of the denominator's ${\sigma }^2{\chi }_{n-p}^2$. Setting $A=I$ and $c=0$ recovers the overall regression $F$-statistic, and a single-row $A$ recovers $t_j^2=F_{1,n-p}$.

When the homoscedasticity assumption fails — $\mathrm{Cov}(\epsilon )={\sigma }^{2}V$ for known positive-definite $V$ — OLS remains unbiased but loses optimality. Aitken's theorem restores it: the GLS estimator ${\hat{\beta }}_{\mathrm{GLS}}=(X^{\top }{V}^{-1}X{)}^{-1}{X}^{\top }{V}^{-1}y$ is BLUE, with covariance ${\sigma }^2(X^{\top }{V}^{-1}X{)}^{-1}⪯{\sigma }^2(X^{\top }X{)}^{-1}{X}^{\top }VX(X^{\top }X{)}^{-1}$, the latter being the (inflated) true covariance of the OLS estimator under correlated errors. Feasible GLS plugs an estimate $\hat{V}$ into the formula and is the workhorse for serially correlated and heteroscedastic data.

Synthesis. The hat matrix $H$ is the foundational reason the whole edifice coheres: it is simultaneously the least-squares solver, the variance-decomposition operator, and the independence-generator. Putting these together, the decomposition $y=Hy+My$ into orthogonal projections is dual to the additive split $\hat{\beta }\perp \mathrm{RSS}$, and this is exactly what licenses every $t$- and $F$-statistic — the central insight is that normality plus orthogonality equals independence, so the inference layer is forced by the geometry rather than assumed. The same projection generalises: replacing $I$ by $V^{-1}$ in the inner product turns OLS into GLS and Gauss-Markov into Aitken, and the bridge is the whitening transform $L^{-1}$ that maps the correlated problem back to the canonical one. This theory grounds the regularized estimators of 45.06.03, where shrinking $\hat{\beta }$ trades the BLUE's unbiasedness for lower variance, and the discriminant and additive-model methods of 45.08.01 that reuse the projection and quadratic-form distribution theory verbatim.

Full proof set Master

Proposition 1 (Sampling distribution of $\hat{\beta }$). Under $\epsilon \sim N(0,{\sigma }^2{I}_n)$ and full column rank $X$, $\hat{\beta }\sim N(\beta ,{\sigma }^2(X^{\top }X{)}^{-1})$.

Proof. $\hat{\beta }=(X^{\top }X{)}^{-1}{X}^{\top }y=(X^{\top }X{)}^{-1}{X}^{\top }(X\beta +\epsilon )=\beta +(X^{\top }X{)}^{-1}{X}^{\top }\epsilon$. The second term is a fixed linear map applied to the Gaussian $\epsilon$, so $\hat{\beta }$ is Gaussian. Its mean is $\beta +(X^{\top }X{)}^{-1}{X}^{\top }\mathbb{E}[\epsilon ]=\beta$, and its covariance is $(X^{\top }X{)}^{-1}{X}^{\top }({\sigma }^{2}I)X(X^{\top }X{)}^{-1}={\sigma }^2(X^{\top }X{)}^{-1}$. $\square$

Proposition 2 (Distribution of RSS). $\mathrm{RSS}/{\sigma }^2\sim {\chi }_{n-p}^2$.

Proof. From Exercise 5, $\mathrm{RSS}={\epsilon }^{\top }M\epsilon$ with $M=I-H$ symmetric idempotent of rank $n-p$ (rank equals trace for a projection, and $\mathrm{tr}M=n-p$). Write the spectral decomposition $M=Q\Lambda Q^{\top }$ with $Q$ orthogonal and $\Lambda =\mathrm{diag}(1,\dots ,1,0,\dots ,0)$ carrying $n-p$ ones. Put $z=Q^{\top }\epsilon /\sigma$; since $Q$ is orthogonal and $\epsilon /\sigma \sim N(0,I)$, also $z\sim N(0,I)$. Then $\mathrm{RSS}/{\sigma }^2=z^{\top }\Lambda z={\sum }_{k=1}^{n-p}{z}_k^2$, a sum of $n-p$ independent squared standard normals, which is ${\chi }_{n-p}^2$. $\square$

Proposition 3 (Independence of $\hat{\beta }$ and RSS). Under normality, $\hat{\beta }$ and $\mathrm{RSS}$ are independent.

Proof. $\hat{\beta }=\beta +(X^{\top }X{)}^{-1}{X}^{\top }\epsilon$ depends on $\epsilon$ through $X^{\top }\epsilon =X^{\top }H\epsilon$ (since $X^{\top }H=X^{\top }$), i.e. through $H\epsilon$, while $\mathrm{RSS}={\epsilon }^{\top }M\epsilon =(M\epsilon {)}^{\top }(M\epsilon )$ depends on $\epsilon$ through $M\epsilon$. The pair $(H\epsilon ,M\epsilon )$ is jointly Gaussian with cross-covariance ${\sigma }^{2}H{M}^{\top }={\sigma }^{2}H(I-H)=0$. Jointly Gaussian and uncorrelated implies independent, so $H\epsilon \perp M\epsilon$; any measurable function of one is independent of any measurable function of the other, giving $\hat{\beta }\perp \mathrm{RSS}$. $\square$

Proposition 4 ($t$-distribution of the coefficient statistic). $t_j=({\hat{\beta }}_j-{\beta }_j)/(\hat{\sigma }\sqrt{[(X^{\top }X{)}^{-1}{]}_{jj}})\sim t_{n-p}$.

Proof. Let $v_{jj}=[(X^{\top }X{)}^{-1}{]}_{jj}$. By Proposition 1, $Z:=({\hat{\beta }}_j-{\beta }_j)/(\sigma \sqrt{v_{jj}})\sim N(0,1)$. By Proposition 2, $W:=\mathrm{RSS}/{\sigma }^2\sim {\chi }_{n-p}^2$, and by Proposition 3, $Z\perp W$. Now $t_j=({\hat{\beta }}_j-{\beta }_j)/(\hat{\sigma }\sqrt{v_{jj}})=Z/\sqrt{W/(n-p)}$, a standard normal divided by the square root of an independent chi-squared over its degrees of freedom — the definition of Student's $t_{n-p}$. $\square$

Proposition 5 (Aitken / GLS optimality). If $\mathrm{Cov}(\epsilon )={\sigma }^{2}V$ with $V\succ 0$ known, then ${\hat{\beta }}_{\mathrm{GLS}}=(X^{\top }{V}^{-1}X{)}^{-1}{X}^{\top }{V}^{-1}y$ is BLUE.

Proof. Factor $V=L{L}^{\top }$ (Cholesky, $L$ invertible) and set $y^{\ast }=L^{-1}y$, $X^{\ast }=L^{-1}X$, ${\epsilon }^{\ast }=L^{-1}\epsilon$. Then $y^{\ast }=X^{\ast }\beta +{\epsilon }^{\ast }$ with $\mathrm{Cov}({\epsilon }^{\ast })=L^{-1}({\sigma }^{2}V)L^{-\top }={\sigma }^2{L}^{-1}L{L}^{\top }{L}^{-\top }={\sigma }^{2}I$. The starred model meets the Gauss-Markov assumptions, so by the theorem its OLS estimator $(X^{\ast \top }{X}^{\ast }{)}^{-1}{X}^{\ast \top }{y}^{\ast }$ is BLUE in the starred parametrization. Substituting $X^{\ast }=L^{-1}X$ and $y^{\ast }=L^{-1}y$ and using $L^{-\top }{L}^{-1}=V^{-1}$ gives $(X^{\top }{V}^{-1}X{)}^{-1}{X}^{\top }{V}^{-1}y={\hat{\beta }}_{\mathrm{GLS}}$. Because $y\mapsto L^{-1}y$ is an invertible linear map, the class of linear unbiased estimators is preserved, so BLUE-ness transfers back to the original model. $\square$

Connections Master

• The hat matrix is the orthogonal projection of the spectral/projection theory 01.01.13; the idempotent self-adjoint characterization of projections used in Propositions 2 and 3 is precisely the finite-dimensional spectral statement, here carrying the variance decomposition.

• The elementary regression survey 26.06.01 introduces the slope, intercept, and $R^2$ at the scatter-plot level; this unit supplies the matrix-form estimator, the Gauss-Markov optimality proof, and the exact sampling distributions that the survey states without derivation.

• The random-variable and multivariate-normal machinery of 26.03.01 — expectation, covariance, the Gaussian transformation rule $\mathrm{Cov}(Ay)=A\mathrm{Cov}(y)A^{\top }$ — is the engine behind every distributional result here.

• The numerical least-squares unit 43.04.01 solves the same normal equations $X^{\top }X\hat{\beta }=X^{\top }y$ but worries about conditioning and backward stability (QR and SVD over normal equations); this unit takes the solution as given and studies its statistical sampling behaviour, so the two are the computational and inferential faces of one object.

• Ridge and lasso regularization 45.06.03 deliberately relax the unbiasedness constraint $CX=I$ that powers the Gauss-Markov proof, trading the BLUE property for a mean-squared-error gain in the high-dimensional regime where $X^{\top }X$ is ill-conditioned or singular.

• Discriminant analysis and generalized additive models 45.08.01 reuse the projection and quadratic-form distribution theory: linear discriminant boundaries and additive-model backfitting both rest on the column-space projection and the chi-squared distribution of residual quadratic forms developed here.

Historical & philosophical context Master

The least-squares principle was published by Legendre in 1805 and developed probabilistically by Gauss, whose Theoria Motus Corporum Coelestium ^{[Gauss 1809]} tied it to the normal error law and to the determination of celestial orbits. Gauss returned to the problem in the Theoria Combinationis Observationum (1823), where he proved the optimality of least squares without assuming normality — the result now called the Gauss-Markov theorem. The attribution to Markov reflects A. A. Markov's early-twentieth-century textbook exposition, which brought Gauss's argument to a wider audience; the optimality claim is Gauss's.

The distributional theory came a century later. R. A. Fisher's work in the 1920s established the sampling distributions of regression statistics and named the analysis of variance, decomposing the total sum of squares into orthogonal components — the decomposition this unit casts as the projections $H$ and $I-H$. W. G. Cochran ^{[Cochran 1934]} proved the theorem on the distribution of sums of quadratic forms in normal variables that underwrites the independence of $\hat{\beta }$ and the residual sum of squares. A. C. Aitken ^{[Aitken 1935]} extended the Gauss-Markov optimality to correlated errors, giving the generalized-least-squares estimator. The modern matrix-projection formulation, the leverage diagnostics from the hat-matrix diagonal, and the unified treatment via the Moore-Penrose inverse for rank-deficient designs were consolidated by C. R. Rao and by Seber in the second half of the century.

Bibliography Master

@book{hastie2009elements,
  author    = {Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome},
  title     = {The Elements of Statistical Learning: Data Mining, Inference, and Prediction},
  edition   = {2},
  publisher = {Springer},
  year      = {2009},
  note      = {Chapter 3, linear methods for regression}
}

@book{seber2003linear,
  author    = {Seber, George A. F. and Lee, Alan J.},
  title     = {Linear Regression Analysis},
  edition   = {2},
  publisher = {Wiley},
  year      = {2003}
}

@book{rao1973linear,
  author    = {Rao, C. Radhakrishna},
  title     = {Linear Statistical Inference and Its Applications},
  edition   = {2},
  publisher = {Wiley},
  year      = {1973}
}

@book{gauss1809theoria,
  author    = {Gauss, Carl Friedrich},
  title     = {Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium},
  publisher = {Perthes and Besser},
  year      = {1809}
}

@article{cochran1934distribution,
  author  = {Cochran, William G.},
  title   = {The distribution of quadratic forms in a normal system, with applications to the analysis of covariance},
  journal = {Mathematical Proceedings of the Cambridge Philosophical Society},
  volume  = {30},
  number  = {2},
  pages   = {178--191},
  year    = {1934}
}

@article{aitken1935least,
  author  = {Aitken, A. C.},
  title   = {On least squares and linear combination of observations},
  journal = {Proceedings of the Royal Society of Edinburgh},
  volume  = {55},
  pages   = {42--48},
  year    = {1935}
}