Gauss-Markov and the Bias-Variance Tradeoff
Math 123: Mathematical Aspects of Data Analysis – Fall 2026
Recap
Last time we put probability into the linear model \(\bfy=\bfX\bftheta+\bfvarepsilon\). Under the standard assumptions (A1)–(A3), we showed that least squares is unbiased with covariance \[ \Ebb[\widehat{\bftheta}_{\mathrm{LS}}]=\bftheta, \qquad \Cov(\widehat{\bftheta}_{\mathrm{LS}})=\sigma^2(\bfX^\top\bfX)^{-1}. \] We ended by asking where in \(\Rbb^d\) this uncertainty lives, and the answer is again the SVD.
Throughout, \(\bfX\in\Rbb^{n\times d}\) has rank \(r\geq1\) and compact SVD \(\bfX=\sum_{i=1}^r\sigma_i\bfu_i\bfv_i^\top\) with \(\sigma_1\geq\cdots\geq\sigma_r>0\). Under (A3), \(r=d\) and \(\bfv_1,\ldots,\bfv_d\) is an orthonormal basis of \(\Rbb^d\).
Where the least squares error lives
We will use the following fact twice today, so we isolate it.
Proof. By linearity of expectation, \(\Ebb\big[\|\bfz-\bfmu\|_2^2\big]=\sum_{j=1}^m\Ebb[(z_j-\mu_j)^2]=\sum_{j=1}^m\Var(z_j)\). These variances are the diagonal entries of \(\Cov(\bfz)\), so their sum is its trace.
Proof. Since \(r=d\), the matrix \(\bfV_d\) is square with orthonormal columns, and by the SVD and eigendecomposition, \[ \bfX^\top\bfX=\bfV_d\Sigma_d^2\bfV_d^\top, \qquad (\bfX^\top\bfX)^{-1}=\bfV_d\Sigma_d^{-2}\bfV_d^\top=\sum_{j=1}^d\frac{1}{\sigma_j^2}\bfv_j\bfv_j^\top, \] exactly as in the proof of the hat matrix. The scalar \(\bfv_i^\top\widehat{\bftheta}_{\mathrm{LS}}\) is a fixed linear map applied to \(\widehat{\bftheta}_{\mathrm{LS}}\), so by means and covariances under linear maps, \[ \Var(\bfv_i^\top\widehat{\bftheta}_{\mathrm{LS}}) =\bfv_i^\top\Cov(\widehat{\bftheta}_{\mathrm{LS}})\bfv_i =\sigma^2\sum_{j=1}^d\frac{(\bfv_i^\top\bfv_j)^2}{\sigma_j^2} =\frac{\sigma^2}{\sigma_i^2}, \] since \(\bfv_i^\top\bfv_j\) is \(1\) if \(j=i\) and \(0\) otherwise.
For the second claim, \(\widehat{\bftheta}_{\mathrm{LS}}\) has mean \(\bftheta\), so by Lemma 1 and linearity of the trace, \[ \Ebb\big[\|\widehat{\bftheta}_{\mathrm{LS}}-\bftheta\|_2^2\big] =\operatorname{tr}\Cov(\widehat{\bftheta}_{\mathrm{LS}}) =\sigma^2\sum_{j=1}^d\frac{1}{\sigma_j^2}\operatorname{tr}(\bfv_j\bfv_j^\top) =\sigma^2\sum_{j=1}^d\frac{1}{\sigma_j^2}, \] since \(\operatorname{tr}(\bfv_j\bfv_j^\top)=\|\bfv_j\|_2^2=1\).
This is the statistical face of the warning about small singular values. 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.
The Gauss-Markov theorem
Least squares is unbiased, but so are other estimators. Is it the best one? We compare it against the following family.
Least squares is linear and unbiased, with \(\bfB=(\bfX^\top\bfX)^{-1}\bfX^\top\).
Proof. Write \(\widetilde{\bftheta}=\bfB\bfy\). By means and covariances under linear maps, \(\Ebb[\widetilde{\bftheta}]=\bfB\bfX\bftheta\), and this must equal \(\bftheta\) for every \(\bftheta\). Taking \(\bftheta=\bfe_j\) shows that column \(j\) of \(\bfB\bfX\) is \(\bfe_j\), so \(\bfB\bfX=\bfI_d\).
Now write \(\bfB=(\bfX^\top\bfX)^{-1}\bfX^\top+\bfC\), where \(\bfC\) is just the difference. Then \[ \bfC\bfX=\bfB\bfX-(\bfX^\top\bfX)^{-1}\bfX^\top\bfX=\bfI_d-\bfI_d=\mathbf{0}. \] By means and covariances under linear maps again, \[ \begin{aligned} \Cov(\widetilde{\bftheta}) &=\sigma^2\bfB\bfB^\top\\ &=\sigma^2\big[(\bfX^\top\bfX)^{-1}\bfX^\top+\bfC\big]\big[\bfX(\bfX^\top\bfX)^{-1}+\bfC^\top\big]\\ &=\sigma^2\big[(\bfX^\top\bfX)^{-1}+(\bfX^\top\bfX)^{-1}(\bfC\bfX)^\top+\bfC\bfX(\bfX^\top\bfX)^{-1}+\bfC\bfC^\top\big]\\ &=\Cov(\widehat{\bftheta}_{\mathrm{LS}})+\sigma^2\bfC\bfC^\top, \end{aligned} \] since \(\bfC\bfX=\mathbf{0}\) kills the two middle terms. Finally, for every \(\bfa\in\Rbb^d\), \[ \Var(\bfa^\top\widetilde{\bftheta})-\Var(\bfa^\top\widehat{\bftheta}_{\mathrm{LS}}) =\sigma^2\bfa^\top\bfC\bfC^\top\bfa =\sigma^2\|\bfC^\top\bfa\|_2^2 \geq0. \]
Taking \(\bfa=\bfe_j\) shows each coefficient of least squares has the smallest possible variance, and summing over \(j\) (Lemma 1) shows least squares also has the smallest expected squared error \(\Ebb\|\widetilde{\bftheta}-\bftheta\|_2^2\) among linear unbiased estimators.
For example, suppose we throw away some rows of \(\bfX\) and \(\bfy\) and run least squares on the rest. This is again linear and unbiased (as long as what remains has full column rank), so Gauss-Markov says it can only be worse. Throwing away data never helps, which is reassuring.
But the fine print matters: the competition is only among estimators that are both linear and unbiased. And by Corollary 1, the winner can still be terrible when some \(\sigma_i\) is small. The rest of today is about what we gain by giving up unbiasedness.
Bias, variance, and mean squared error
The following definitions make sense for any estimator, linear or not.
So Corollary 1 computed \(\operatorname{MSE}(\widehat{\bftheta}_{\mathrm{LS}})=\sigma^2\sum_i\sigma_i^{-2}\).
Proof. Let \(\bar{\bftheta}=\Ebb[\widehat{\bftheta}]\) and split the error into a random part and a fixed part, \[ \widehat{\bftheta}-\bftheta=(\widehat{\bftheta}-\bar{\bftheta})+(\bar{\bftheta}-\bftheta). \] Expanding the square, \[ \|\widehat{\bftheta}-\bftheta\|_2^2 =\|\widehat{\bftheta}-\bar{\bftheta}\|_2^2 +2(\bar{\bftheta}-\bftheta)^\top(\widehat{\bftheta}-\bar{\bftheta}) +\|\bar{\bftheta}-\bftheta\|_2^2. \] Take expectations term by term. The first term gives \(\operatorname{tr}\Cov(\widehat{\bftheta})\) by Lemma 1. In the middle term \(\bar{\bftheta}-\bftheta\) is a fixed vector, so by linearity its expectation is \(2(\bar{\bftheta}-\bftheta)^\top(\Ebb[\widehat{\bftheta}]-\bar{\bftheta})=0\). The last term is not random at all, and equals \(\|\operatorname{Bias}(\widehat{\bftheta})\|_2^2\).
This is the same computation as \(\Ebb[(Y-c)^2]=\Var(Y)+(\Ebb[Y]-c)^2\) for a scalar random variable, done one coordinate at a time. Nothing about the linear model was used: the decomposition holds for any estimator and any noise.
An unbiased estimator spends its whole MSE budget on variance. The question is whether we can buy a large drop in variance for a small amount of bias. Here is the simplest possible case.
The catch is that \(c^\ast\) depends on the unknown \(\theta\). This is the general shape of the story: regularization buys lower error, at the price of a tuning parameter that we have to choose from data. Note the ratio \(c^\ast=\mathrm{SNR}/(1+\mathrm{SNR})\), where \(\mathrm{SNR}=\theta^2/(\sigma^2/n)\) is a signal-to-noise ratio; we will see it again.
Prediction error
Often we care less about \(\bftheta\) itself than about predicting a new response. Suppose a new observation arrives at a fixed feature vector \(\bfx_0\), with \(y_0=\bfx_0^\top\bftheta+\varepsilon_0\), where \(\varepsilon_0\) has mean \(0\) and variance \(\sigma^2\) and is independent of the training noise \(\bfvarepsilon\). We predict it with \(\widehat{y}_0=\bfx_0^\top\widehat{\bftheta}\).
Proof. Write \(y_0-\bfx_0^\top\widehat{\bftheta}=\varepsilon_0-\bfx_0^\top(\widehat{\bftheta}-\bftheta)\) and expand the square. The cross term has expectation \(-2\,\Ebb[\varepsilon_0]\,\Ebb[\bfx_0^\top(\widehat{\bftheta}-\bftheta)]=0\), by independence and expectations of products. What remains is \(\Ebb[\varepsilon_0^2]=\sigma^2\) plus \(\Ebb[(\bfx_0^\top\widehat{\bftheta}-\bfx_0^\top\bftheta)^2]\), which is the MSE of the scalar estimator \(\bfx_0^\top\widehat{\bftheta}\) of \(\bfx_0^\top\bftheta\). Apply Theorem 2 with \(d=1\).
The three terms are three different sources of error.
- Bias is systematic error from the method, or from a model that is wrong, for instance by leaving out the \(t^2\) term in the Keeling curve.
- Variance is error from the particular noisy sample we happened to get. For least squares it is \(\sigma^2\bfx_0^\top(\bfX^\top\bfX)^{-1}\bfx_0=\sigma^2\sum_i(\bfv_i^\top\bfx_0)^2/\sigma_i^2\), which is large when \(\bfx_0\) points along a poorly measured direction (extrapolation).
- Irreducible error \(\sigma^2\) is noise in the new response itself, and no estimator can remove it.
Underfitting and overfitting
A flexible model (many features, for instance a high-degree polynomial in \(t\)) has small bias but large variance: it can bend to fit the noise in the particular sample, which is overfitting. A rigid model (few features) has small variance but may have large bias: it cannot capture the real structure, which is underfitting.
The training error \(\frac1n\|\bfy-\bfX\widehat{\bftheta}\|_2^2\) is computed on the same data used to fit \(\widehat{\bftheta}\), so it is optimistic. With \(n\) distinct points, a polynomial of degree \(n-1\) interpolates the data, so its training error is exactly \(0\), while its predictions between and beyond the data can be wild. What we actually care about is the test error, the expected error on a new observation as in Corollary 2. The difference between the two is called the generalization gap.
Summary
- Least squares has variance \(\sigma^2/\sigma_i^2\) along each \(\bfv_i\) (Corollary 1), so small singular values make it noisy.
- Even so, it is the best linear unbiased estimator (Theorem 1).
- For any estimator, the mean squared error splits into squared bias plus variance (Theorem 2), and prediction error adds an irreducible \(\sigma^2\) (Corollary 2).
- Giving up unbiasedness can lower the mean squared error, as in Example 2.
Next time, ridge regression does this systematically for the linear model, by shrinking the poorly measured directions.