The Statistics of Least Squares

Math 123: Mathematical Aspects of Data Analysis – Fall 2026

Published

October 1, 2026

Recap

Last time we set up the linear model \(\bfy=\bfX\bftheta+\bfvarepsilon\) and fit it by least squares. We discovered that geometrically, least squares is orthogonal projection. That is, the fitted values \(\widehat{\bfy}\) are the projection of \(\bfy\) onto \(\Ccal(\bfX)\), and the residual is orthogonal to \(\Ccal(\bfX)\). From there, we derived the normal equations, which told us that \(\widehat{\bftheta}_{\mathrm{LS}}=(\bfX^\top\bfX)^{-1}\bfX^\top\bfy\) when \(\bfX\) has full column rank.

The pseudoinverse and minimum-norm least squares

When \(\operatorname{rank}(\bfX)<d\) there are infinitely many least squares solutions, forming the affine set \(\widehat{\bftheta}+\Ncal(\bfX)\) (the set of least squares solutions). A natural way to pick one is to take the shortest, and the SVD produces it. Throughout this section, \(\bfX\in\Rbb^{n\times d}\) has rank \(r\geq1\) and compact SVD \[ \bfX=\bfU_r\Sigma_r\bfV_r^\top=\sum_{i=1}^r\sigma_i\bfu_i\bfv_i^\top, \qquad \sigma_1\geq\cdots\geq\sigma_r>0. \]

Definition 1 (Moore-Penrose pseudoinverse) The (Moore-Penrose) pseudoinverse of \(\bfX\) is \[ \bfX^\dagger=\bfV_r\Sigma_r^{-1}\bfU_r^\top =\sum_{i=1}^r\frac{1}{\sigma_i}\bfv_i\bfu_i^\top \in\Rbb^{d\times n}. \] (For \(\bfX=\mathbf{0}\) we set \(\bfX^\dagger=\mathbf{0}\).)

The pseudoinverse undoes \(\bfX\) one singular direction at a time: \[ \bfX\bfv_i=\sigma_i\bfu_i, \qquad \bfX^\dagger\bfu_i=\frac{1}{\sigma_i}\bfv_i, \qquad \bfX^\dagger\bfw=\mathbf{0}\ \text{for}\ \bfw\perp\Ccal(\bfX). \]

Theorem 1 (Minimum-norm least squares) For every \(\bfy\in\Rbb^n\), the vector \[ \bftheta^\dagger =\bfX^\dagger\bfy =\sum_{i=1}^r\frac{\bfu_i^\top\bfy}{\sigma_i}\,\bfv_i \] is a least squares solution, and it is the unique least squares solution of smallest Euclidean norm: \(\|\bftheta^\dagger\|_2<\|\bftheta\|_2\) for every other least squares solution \(\bftheta\). In particular, if \(\operatorname{rank}(\bfX)=d\), then \(\bfX^\dagger=(\bfX^\top\bfX)^{-1}\bfX^\top\) and \(\bftheta^\dagger=\widehat{\bftheta}_{\mathrm{LS}}\).

Proof

Proof. First, \(\bftheta^\dagger\) is a least squares solution, since \[ \bfX\bftheta^\dagger =\bfU_r\Sigma_r\bfV_r^\top\bfV_r\Sigma_r^{-1}\bfU_r^\top\bfy =\bfU_r\bfU_r^\top\bfy =\bfP\bfy, \] using \(\bfV_r^\top\bfV_r=\bfI_r\) and the column space from the SVD, part 1, and then least squares is orthogonal projection, part 2.

By the set of least squares solutions, every other least squares solution has the form \(\bftheta^\dagger+\bfz\) with \(\bfz\in\Ncal(\bfX)\), \(\bfz\neq\mathbf{0}\). Now \(\bftheta^\dagger\) is a combination of \(\bfv_1,\ldots,\bfv_r\), while \(\bfz\) is orthogonal to all of them by the null space from the SVD, part 2. So by the Pythagorean theorem, \[ \|\bftheta^\dagger+\bfz\|_2^2 =\|\bftheta^\dagger\|_2^2+\|\bfz\|_2^2 >\|\bftheta^\dagger\|_2^2. \]

Finally, if \(\operatorname{rank}(\bfX)=d\), the least squares solution is unique, so \(\bfX^\dagger\bfy=\widehat{\bftheta}_{\mathrm{LS}}=(\bfX^\top\bfX)^{-1}\bfX^\top\bfy\) for every \(\bfy\). Taking \(\bfy=\bfe_j\) shows the two matrices have the same \(j\)-th column, for every \(j\).

The catch is the factor \(1/\sigma_i\): it is harmless when \(\sigma_i\) is large, but when some \(\sigma_i\) is small, a little noise in \(\bfy\) along \(\bfu_i\) becomes a large change in \(\bftheta^\dagger\) along \(\bfv_i\) (small singular values amplify noise). We now make this precise statistically.

The Statistics of Least Squares

So far there has been no probability: least squares is a geometric recipe that makes sense for any \(\bfX\) and \(\bfy\). Now we ask a statistical question. If the data really were generated by \(\bfy=\bfX\bftheta+\bfvarepsilon\) with random noise, how good is \(\widehat{\bftheta}_{\mathrm{LS}}\) as an estimate of the true \(\bftheta\)?

We will treat the design matrix \(\bfX\) as fixed (not random), so all of the randomness comes from the noise \(\bfvarepsilon\). This matches the regression assumption that the features are measured exactly. The textbook allows \(\bfX\) to be random and writes everything “given \(\bfX\)”, as in \(\Ebb[\bfvarepsilon\mid\bfX]\); if you have not seen conditional expectation, read “given \(\bfX\)” as “with \(\bfX\) held fixed”, which is what we do here.

Random vectors

We need one small extension of Lectures 1 and 2: means and covariances of random vectors.

Definition 2 (Mean and covariance of a random vector) A random vector \(\bfz=(z_1,\ldots,z_m)^\top\) is a vector whose entries are random variables. Its mean is the vector of means, and its covariance matrix collects all the pairwise covariances: \[ \Ebb[\bfz]= \begin{pmatrix} \Ebb[z_1]\\ \vdots\\ \Ebb[z_m] \end{pmatrix}\in\Rbb^m, \qquad \Cov(\bfz)=\Ebb\big[(\bfz-\Ebb[\bfz])(\bfz-\Ebb[\bfz])^\top\big]\in\Rbb^{m\times m}, \] where the expectation of a random matrix is taken entry by entry. The \((i,j)\) entry of \(\Cov(\bfz)\) is \(\Cov(z_i,z_j)\), and its diagonal entries are the variances \(\Var(z_i)\).

Proposition 1 (Means and covariances under linear maps) Let \(\bfz\in\Rbb^m\) be a random vector, and let \(\bfA\in\Rbb^{k\times m}\) and \(\bfb\in\Rbb^k\) be fixed. Then \[ \Ebb[\bfA\bfz+\bfb]=\bfA\,\Ebb[\bfz]+\bfb, \qquad \Cov(\bfA\bfz+\bfb)=\bfA\Cov(\bfz)\bfA^\top. \] In particular, for a fixed \(\bfa\in\Rbb^m\), \(\Var(\bfa^\top\bfz)=\bfa^\top\Cov(\bfz)\bfa\geq0\), so \(\Cov(\bfz)\) is positive semidefinite.

Proof

Proof. See Homework 4, Problem 3.

The standard assumptions

Definition 3 (Standard assumptions for the linear model) The design matrix \(\bfX\in\Rbb^{n\times d}\) is fixed, \(\bfy=\bfX\bftheta+\bfvarepsilon\) for a fixed unknown \(\bftheta\in\Rbb^d\), and:

(A1) Zero-mean noise: \(\Ebb[\bfvarepsilon]=\mathbf{0}\).
(A2) Homoscedastic, uncorrelated noise: \(\Cov(\bfvarepsilon)=\sigma^2\bfI_n\) for some \(\sigma^2>0\). That is, \(\Var(\varepsilon_i)=\sigma^2\) for every \(i\), and \(\Cov(\varepsilon_i,\varepsilon_j)=0\) for \(i\neq j\).
(A3) Full column rank: \(\operatorname{rank}(\bfX)=d\) (so \(n\geq d\)).

Each assumption has a job.

(A1) says the noise has no systematic part: on average the model is right, \(\Ebb[y_i]=\bfx_i^\top\bftheta\). This is the most important assumption, and it fails if we leave out a relevant feature. If we drop the \(t^2\) term from the Keeling model, the acceleration is swept into \(\varepsilon_i\) and the noise no longer averages to zero: the data sit above a straight-line fit at both ends and below it in the middle. (In the random-\(\bfX\) language of the textbook, this is called exogeneity: the noise is not correlated with the features.)
(A2) says every observation is equally noisy, and the noise in one observation says nothing about another. This is an idealization. For example, if we regress income on years of education, the spread of incomes is much larger among people with 16 years of education than among those with 4; this is called heteroscedasticity.
(A3) we met already: it makes the least squares solution unique.

If the \(\varepsilon_i\) are independent with mean \(0\) and variance \(\sigma^2\), then (A1) and (A2) hold, since independent random variables have covariance \(0\).

By Proposition 1, the response is itself a random vector with \[ \Ebb[\bfy]=\bfX\bftheta, \qquad \Cov(\bfy)=\Cov(\bfvarepsilon)=\sigma^2\bfI_n. \]

Least squares is unbiased

Theorem 2 (Mean and covariance of the least squares estimator) Under (A1)–(A3), the least squares estimator \(\widehat{\bftheta}_{\mathrm{LS}}=(\bfX^\top\bfX)^{-1}\bfX^\top\bfy\) satisfies

(Unbiasedness) \(\Ebb[\widehat{\bftheta}_{\mathrm{LS}}]=\bftheta\);
(Covariance) \(\Cov(\widehat{\bftheta}_{\mathrm{LS}})=\sigma^2(\bfX^\top\bfX)^{-1}\).

Proof

Proof. Write \(\widehat{\bftheta}_{\mathrm{LS}}=\bfA\bfy\) with the fixed matrix \(\bfA=(\bfX^\top\bfX)^{-1}\bfX^\top\), and apply Proposition 1. For the mean, \[ \Ebb[\widehat{\bftheta}_{\mathrm{LS}}] =\bfA\,\Ebb[\bfy] =(\bfX^\top\bfX)^{-1}\bfX^\top\bfX\bftheta =\bftheta. \] For the covariance, using that \((\bfX^\top\bfX)^{-1}\) is symmetric, \[ \Cov(\widehat{\bftheta}_{\mathrm{LS}}) =\bfA(\sigma^2\bfI_n)\bfA^\top =\sigma^2(\bfX^\top\bfX)^{-1}\bfX^\top\bfX(\bfX^\top\bfX)^{-1} =\sigma^2(\bfX^\top\bfX)^{-1}. \]

Unbiased means that if we could repeat the experiment many times with fresh noise and average the resulting estimates, we would land exactly on \(\bftheta\). It does not mean any single estimate is accurate; the covariance tells us how far a typical estimate wanders.

Example 1 (The variance of the sample mean, again) In fitting a constant, \(\bfX=\mathbf{1}\), so \(\bfX^\top\bfX=n\) and \(\widehat{\theta}_{\mathrm{LS}}=\bar{y}\). Theorem 2 says \(\Ebb[\bar{y}]=\theta\) and \(\Var(\bar{y})=\sigma^2/n\). This is averaging away independent noise from Lecture 2, which the theorem generalizes to any design matrix.

The SVD tells us where the uncertainty lives.

Corollary 1 (Where the least squares error lives) Under (A1)–(A3), with compact SVD \(\bfX=\sum_{i=1}^d\sigma_i\bfu_i\bfv_i^\top\), \[ \Var(\bfv_i^\top\widehat{\bftheta}_{\mathrm{LS}})=\frac{\sigma^2}{\sigma_i^2}, \qquad \Ebb\big[\|\widehat{\bftheta}_{\mathrm{LS}}-\bftheta\|_2^2\big]=\sigma^2\sum_{i=1}^d\frac{1}{\sigma_i^2}. \]

Proof

Proof. From the proof of the hat matrix, \(\Cov(\widehat{\bftheta}_{\mathrm{LS}})=\sigma^2\bfV_d\Sigma_d^{-2}\bfV_d^\top\). By Proposition 1 and \(\bfV_d^\top\bfv_i=\bfe_i\), \[ \Var(\bfv_i^\top\widehat{\bftheta}_{\mathrm{LS}}) =\sigma^2\bfv_i^\top\bfV_d\Sigma_d^{-2}\bfV_d^\top\bfv_i =\sigma^2\bfe_i^\top\Sigma_d^{-2}\bfe_i =\frac{\sigma^2}{\sigma_i^2}. \] For the second claim, \(\widehat{\bftheta}_{\mathrm{LS}}-\bftheta\) has mean zero, so its covariance is \(\Ebb[(\widehat{\bftheta}_{\mathrm{LS}}-\bftheta)(\widehat{\bftheta}_{\mathrm{LS}}-\bftheta)^\top]\). Using \(\|\bfw\|_2^2=\operatorname{tr}(\bfw\bfw^\top)\) (properties of the trace, part 4), and that expectation and trace are both linear and so can be swapped, \[ \Ebb\big[\|\widehat{\bftheta}_{\mathrm{LS}}-\bftheta\|_2^2\big] =\operatorname{tr}\Cov(\widehat{\bftheta}_{\mathrm{LS}}) =\sigma^2\operatorname{tr}(\bfV_d\Sigma_d^{-2}\bfV_d^\top) =\sigma^2\operatorname{tr}(\Sigma_d^{-2}\bfV_d^\top\bfV_d) =\sigma^2\sum_{i=1}^d\frac{1}{\sigma_i^2}, \] by the cyclic property of the trace.

This is the statistical face of the warning about small singular values from last time. Along a direction \(\bfv_i\) where the data are spread out (\(\sigma_i\) large), the estimate is precise. Along a direction where the data barely vary (\(\sigma_i\) small), the estimate is very noisy, and a single tiny \(\sigma_d\) is enough to make the total error huge. Collecting more data helps because it makes the \(\sigma_i\) grow: duplicating every row of \(\bfX\) doubles \(\bfX^\top\bfX\) and halves the error. See the example The Variance of Least Squares for a simulation of deblurring a photograph, where most of the singular values are tiny.

Summary: why squares?

Viewpoint	What least squares is
Geometry	The orthogonal projection of \(\bfy\) onto \(\Ccal(\bfX)\) (least squares is orthogonal projection).
Algebra	A solution of the normal equations \(\bfX^\top\bfX\bftheta=\bfX^\top\bfy\) (normal equations), with minimum-norm solution \(\bfX^\dagger\bfy\) (Theorem 1).
Statistics	Unbiased (Theorem 2), with variance \(\sigma^2/\sigma_i^2\) along each \(\bfv_i\) (Corollary 1).

Small singular values of \(\bfX\) are trouble, though. They blow up \(1/\sigma_i\) in the pseudoinverse, they square into a terrible condition number in the normal equations, and they produce estimates with variance \(\sigma^2/\sigma_i^2\). Next time we will prove that least squares is nonetheless the best linear unbiased estimator, and after that we will treat small singular values with ridge regression.