Ridge Regression
Math 123: Mathematical Aspects of Data Analysis – Fall 2026
Recap
Last time we proved the Gauss–Markov theorem: under the standard assumptions (A1)–(A3), least squares is the best linear unbiased estimator. Even so, its mean squared error is \(\sigma^2\sum_i\sigma_i^{-2}\) (where the least squares error lives), which is huge when some \(\sigma_i\) is small. The bias–variance decomposition \[ \operatorname{MSE}(\widehat{\bftheta}) =\|\operatorname{Bias}(\widehat{\bftheta})\|_2^2 +\operatorname{tr}\Cov(\widehat{\bftheta}) \] suggests giving up a little unbiasedness in exchange for a large drop in variance, as when shrinking the sample mean. Let’s try it and see what happens!
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\).
Ridge regression
Least squares fails, or becomes unstable, in three related situations: redundant features (\(\operatorname{rank}(\bfX)<d\)), more features than observations (\(d>n\)), and nearly redundant features (a tiny \(\sigma_d\)). In each case the problem is directions \(\bfv_i\) in coefficient space that the data barely see. Ridge regression adds a penalty that keeps the coefficients from exploding along those directions.
The objective balances two goals: fit the data (small \(\|\bfy-\bfX\bftheta\|_2^2\)), and keep the coefficients small (small \(\|\bftheta\|_2^2\)). The parameter \(\lambda\) sets the exchange rate between them. In numerical analysis the same idea is called Tikhonov regularization.
The matrix \(\bfX^\top\bfX+\lambda\bfI_d\) will appear everywhere below, so we first work out its SVD.
Proof. From the compact SVD and \(\bfu_i^\top\bfu_j=0\) for \(i\neq j\), \(\bfu_i^\top\bfu_i=1\), \[ \bfX^\top\bfX =\sum_{i,j=1}^r\sigma_i\sigma_j\,\bfv_i(\bfu_i^\top\bfu_j)\bfv_j^\top =\sum_{i=1}^r\sigma_i^2\,\bfv_i\bfv_i^\top =\sum_{i=1}^d\sigma_i^2\,\bfv_i\bfv_i^\top. \] Since \(\bfv_1,\ldots,\bfv_d\) is an orthonormal basis, \(\bfI_d=\sum_{i=1}^d\bfv_i\bfv_i^\top\); adding \(\lambda\) times this gives the first formula. The \(\sigma_i\) are decreasing and the extra ones are \(0\), so the numbers \(\sigma_i^2+\lambda\) are decreasing and at least \(\lambda>0\), and the first formula is an SVD. For the second, multiply the two sums and use \(\bfv_i^\top\bfv_j=0\) for \(i\neq j\) and \(\bfv_i^\top\bfv_i=1\): \[ \Big(\sum_{i=1}^d(\sigma_i^2+\lambda)\,\bfv_i\bfv_i^\top\Big)\Big(\sum_{j=1}^d\frac{1}{\sigma_j^2+\lambda}\,\bfv_j\bfv_j^\top\Big) =\sum_{i=1}^d\frac{\sigma_i^2+\lambda}{\sigma_i^2+\lambda}\,\bfv_i\bfv_i^\top =\bfI_d. \] Listing the terms in reverse order puts the reciprocals in decreasing order, so this is an SVD too.
Proof. Ridge regression is ordinary least squares in disguise. Stack \[ \widetilde{\bfX}=\begin{pmatrix}\bfX\\ \sqrt{\lambda}\,\bfI_d\end{pmatrix}\in\Rbb^{(n+d)\times d}, \qquad \widetilde{\bfy}=\begin{pmatrix}\bfy\\ \mathbf{0}\end{pmatrix}\in\Rbb^{n+d}. \] The squared norm of a stacked vector is the sum of the squared norms of its pieces, so \[ \|\widetilde{\bfy}-\widetilde{\bfX}\bftheta\|_2^2 =\|\bfy-\bfX\bftheta\|_2^2+\|\mathbf{0}-\sqrt{\lambda}\,\bftheta\|_2^2 =\|\bfy-\bfX\bftheta\|_2^2+\lambda\|\bftheta\|_2^2. \] So \(\widehat{\bftheta}_\lambda\) is exactly a least squares solution for \((\widetilde{\bfX},\widetilde{\bfy})\). The stacked matrix has full column rank: if \(\widetilde{\bfX}\bfz=\mathbf{0}\), then its bottom block says \(\sqrt{\lambda}\,\bfz=\mathbf{0}\), so \(\bfz=\mathbf{0}\). By the set of least squares solutions and \(\bfX^\top\bfX\) and \(\bfX\) have the same null space, \(\widetilde{\bfX}^\top\widetilde{\bfX}\) is invertible and the solution is unique, namely \[ (\widetilde{\bfX}^\top\widetilde{\bfX})^{-1}\widetilde{\bfX}^\top\widetilde{\bfy} =(\bfX^\top\bfX+\lambda\bfI_d)^{-1}\bfX^\top\bfy. \]
There’s a good interpretation of what’s happening by looking at the stacked version discussed in the above proof. Ridge regression is least squares after appending \(d\) fake observations that say “each coefficient is \(0\)”, each with weight \(\sqrt{\lambda}\).
Ridge regression through the SVD
Proof. Use the notation of Lemma 1. From the SVD of \(\bfX\), \(\bfX^\top\bfy=\sum_{j=1}^r\sigma_j(\bfu_j^\top\bfy)\,\bfv_j\). So by Theorem 1 and the formula for the inverse in Lemma 1, using \(\bfv_i^\top\bfv_j=0\) for \(i\neq j\) and \(\bfv_i^\top\bfv_i=1\), \[ \widehat{\bftheta}_\lambda =(\bfX^\top\bfX+\lambda\bfI_d)^{-1}\bfX^\top\bfy =\sum_{i=1}^d\sum_{j=1}^r\frac{\sigma_j}{\sigma_i^2+\lambda}(\bfu_j^\top\bfy)\,\bfv_i(\bfv_i^\top\bfv_j) =\sum_{i=1}^r\frac{\sigma_i}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\,\bfv_i. \] The second form in the statement is the same sum, with each term multiplied and divided by \(\sigma_i\).
So ridge regression is the pseudoinverse with each term multiplied by a filter factor \(f_i=\sigma_i^2/(\sigma_i^2+\lambda)\in(0,1)\).
| \(\sigma_i^2\gg\lambda\) | \(\sigma_i^2\ll\lambda\) | |
|---|---|---|
| filter factor \(f_i\) | \(\approx1\): keep this direction | \(\approx\sigma_i^2/\lambda\approx0\): suppress it |
| gain \(\sigma_i/(\sigma_i^2+\lambda)\) | \(\approx1/\sigma_i\), as in least squares | \(\approx\sigma_i/\lambda\), small |
Well-measured directions are left alone and poorly measured ones are shrunk toward zero. In particular the gain never exceeds \(1/(2\sqrt{\lambda})\) (by \(\sigma_i^2+\lambda\geq2\sigma_i\sqrt{\lambda}\)), so noise along \(\bfu_i\) can be amplified by at most \(1/(2\sqrt\lambda)\) instead of \(1/\sigma_i\). This is the cure promised in the warning about small singular values.
Two limits are worth noting.
- As \(\lambda\to0^+\), every \(f_i\to1\), so \(\widehat{\bftheta}_\lambda\to\bftheta^\dagger\), even when \(\bfX\) is rank-deficient. This gives another description of the pseudoinverse: \(\bfX^\dagger=\lim_{\lambda\to0^+}(\bfX^\top\bfX+\lambda\bfI_d)^{-1}\bfX^\top\).
- As \(\lambda\to\infty\), every \(f_i\to0\), so \(\widehat{\bftheta}_\lambda\to\mathbf{0}\).
Ridge also repairs the conditioning: by Lemma 1, \(\bfX^\top\bfX+\lambda\bfI_d\) has singular values \(\sigma_i^2+\lambda\) (and \(\lambda\), if \(r<d\)), so its condition number is at most \((\sigma_1^2+\lambda)/\lambda\), however small \(\sigma_d\) is.
- The penalty \(\lambda\|\bftheta\|_2^2\) treats all coefficients alike, so ridge depends on the units of the features: measuring temperature in millikelvin instead of Celsius changes the answer. Standardize each feature (subtract its mean and divide by its standard deviation) before fitting.
- Do not penalize the intercept \(\theta_0\); shrinking it toward \(0\) just biases the fit toward \(y=0\) for no reason. Centering \(\bfX\) and \(\bfy\) and then fitting without an intercept accomplishes this.
- In scikit-learn,
Ridge(alpha=lam)minimizes exactly \(\|\bfy-\bfX\bftheta\|_2^2+\lambda\|\bftheta\|_2^2\) and leaves the intercept unpenalized.
The bias and variance of ridge regression
Ridge regression is linear, \(\widehat{\bftheta}_\lambda=\bfA_\lambda\bfy\) with \(\bfA_\lambda=(\bfX^\top\bfX+\lambda\bfI_d)^{-1}\bfX^\top\), but it is biased, so Gauss–Markov does not apply to it. To keep the formulas clean, assume (A3), so \(r=d\) and \(\bfX=\bfU_d\Sigma_d\bfV_d^\top\) with \(\bfV_d\in\Rbb^{d\times d}\) orthogonal. Write the true coefficients in the basis \(\bfv_1,\ldots,\bfv_d\) as \(\bfalpha=\bfV_d^\top\bftheta\), that is, \(\alpha_i=\bfv_i^\top\bftheta\).
Proof. Let \(\bfD=(\Sigma_d^2+\lambda\bfI_d)^{-1}\), the diagonal matrix with entries \(1/(\sigma_i^2+\lambda)\). By Lemma 1, \((\bfX^\top\bfX+\lambda\bfI_d)^{-1}=\bfV_d\bfD\bfV_d^\top\), so \(\bfA_\lambda=\bfV_d\bfD\Sigma_d\bfU_d^\top\). Using \(\bfU_d^\top\bfU_d=\bfI_d\) and that diagonal matrices commute, \[ \bfA_\lambda\bfX=\bfV_d\bfD\Sigma_d^2\bfV_d^\top, \qquad \bfA_\lambda\bfA_\lambda^\top=\bfV_d\bfD^2\Sigma_d^2\bfV_d^\top. \] By means and covariances under linear maps, \(\Ebb[\widehat{\bftheta}_\lambda]=\bfA_\lambda\bfX\bftheta=\bfV_d\bfD\Sigma_d^2\bfalpha\) and \(\Cov(\widehat{\bftheta}_\lambda)=\sigma^2\bfA_\lambda\bfA_\lambda^\top\), which are the two formulas written out column by column.
For the MSE, \(\bftheta=\bfV_d\bfalpha\) and \(\bfD\Sigma_d^2-\bfI_d=\bfD(\Sigma_d^2-(\Sigma_d^2+\lambda\bfI_d))=-\lambda\bfD\), so the bias is \[ \Ebb[\widehat{\bftheta}_\lambda]-\bftheta=\bfV_d(\bfD\Sigma_d^2-\bfI_d)\bfalpha=-\lambda\bfV_d\bfD\bfalpha. \] Since \(\bfV_d\) is orthogonal, its squared norm is \(\lambda^2\|\bfD\bfalpha\|_2^2\), the bias term. By the cyclic property of the trace, \(\operatorname{tr}\Cov(\widehat{\bftheta}_\lambda)=\sigma^2\operatorname{tr}(\bfD^2\Sigma_d^2\bfV_d^\top\bfV_d)=\sigma^2\operatorname{tr}(\bfD^2\Sigma_d^2)\), the variance term. Now apply the bias–variance decomposition.
(Without (A3) the same computation goes through with the full SVD, and the bias picks up an extra term \(\sum_{i>r}\alpha_i^2\). That is the part of \(\bftheta\) in \(\Ncal(\bfX)\): the data never see it, so no estimator built from \(\bfX\bftheta\) can recover it.)
Setting \(\lambda=0\) recovers least squares: zero bias and variance \(\sigma^2\sum_i\sigma_i^{-2}\), as in where the least squares error lives. As \(\lambda\) grows, every bias term \(\lambda^2\alpha_i^2/(\sigma_i^2+\lambda)^2\) increases and every variance term \(\sigma^2\sigma_i^2/(\sigma_i^2+\lambda)^2\) decreases. This is the bias–variance tradeoff, and we can make it pay.
Proof. By Proposition 1 the MSE is a sum of \(d\) terms \[ g_i(\lambda)=\frac{\lambda^2\alpha_i^2+\sigma^2\sigma_i^2}{(\sigma_i^2+\lambda)^2}, \qquad g_i(0)=\frac{\sigma^2}{\sigma_i^2}, \] and \(\sum_ig_i(0)=\operatorname{MSE}(\widehat{\bftheta}_{\mathrm{LS}})\). By the quotient rule (check this!), \[ g_i'(\lambda)=\frac{2\sigma_i^2(\lambda\alpha_i^2-\sigma^2)}{(\sigma_i^2+\lambda)^3}. \] This is negative whenever \(\lambda\alpha_i^2<\sigma^2\), so each \(g_i\) is strictly decreasing on \([0,\sigma^2/\alpha_i^2]\) (on all of \([0,\infty)\) if \(\alpha_i=0\)). For \(\lambda\) in the stated range, every \(g_i\) is strictly decreasing on \([0,\lambda]\), so \(g_i(\lambda)<g_i(0)\) for every \(i\). Sum over \(i\).
So even the champion of Gauss–Markov can be beaten, by a biased linear estimator. The improvement is largest exactly when least squares is worst: a small \(\sigma_i\) makes \(g_i(0)=\sigma^2/\sigma_i^2\) huge, and ridge cuts it down. This is not a contradiction, since ridge is not unbiased; it simply shows that unbiasedness is an expensive requirement.
Remark (Stein’s paradox). Something even stranger is true. For Gaussian noise and \(d\geq3\), Charles Stein showed in 1956 that there is an estimator (a shrinkage estimator, now named after James and Stein) whose MSE is strictly smaller than that of least squares for every value of \(\bftheta\), not just for small \(\lambda\) chosen with \(\bftheta\) in mind.
As with the shrunken mean, the good range of \(\lambda\) depends on the unknown \(\bftheta\). Each \(g_i\) is minimized at \(\lambda=\sigma^2/\alpha_i^2\), a noise-to-signal ratio. If all \(\alpha_i^2\) were equal to some \(\tau^2\), the best choice would be \(\lambda^\ast=\sigma^2/\tau^2\), and then the filter factor would be \[ f_i=\frac{\sigma_i^2}{\sigma_i^2+\lambda^\ast}=\frac{\mathrm{SNR}_i}{1+\mathrm{SNR}_i}, \qquad \mathrm{SNR}_i=\frac{\sigma_i^2\tau^2}{\sigma^2}, \] the same form as \(c^\ast\) in shrinking the sample mean. Here \(\mathrm{SNR}_i\) compares the size of the signal \(\sigma_i\alpha_i\) in the \(i\)-th observation direction to the noise level \(\sigma\).
Remark (A Bayesian reading (optional)). The textbook (Section 4.6) makes this precise. If we model \(\bftheta\) itself as random, with independent \(\mathcal{N}(0,\tau^2)\) coordinates, and the noise as \(\mathcal{N}(\mathbf{0},\sigma^2\bfI_n)\), then the most probable \(\bftheta\) given the data is exactly \(\widehat{\bftheta}_\lambda\) with \(\lambda=\sigma^2/\tau^2\). In the SVD coordinates the problem splits into \(r\) independent scalar problems “signal plus noise”, and in each one ridge applies the classical Wiener filter gain \(\mathrm{SNR}_i/(1+\mathrm{SNR}_i)\): trust the data where the signal-to-noise ratio is high, fall back on the prior guess \(0\) where it is low.
Choosing \(\lambda\)
We cannot minimize \(\operatorname{MSE}(\widehat{\bftheta}_\lambda)\) directly, since it depends on the unknown \(\bftheta\) and \(\sigma^2\). In practice \(\lambda\) is chosen from the data, usually by cross-validation; see the supplementary notes on choosing hyperparameters.
Replacing \(\lambda\|\bftheta\|_2^2\) with \(\lambda\|\bftheta\|_1\) gives the lasso. It has no closed-form solution, but it tends to set many coefficients exactly to zero, which selects a small set of features. We will see why later in the course.
Summary
| Least squares | Ridge regression | |
|---|---|---|
| Objective | \(\lVert\bfy-\bfX\bftheta\rVert_2^2\) | \(\lVert\bfy-\bfX\bftheta\rVert_2^2+\lambda\lVert\bftheta\rVert_2^2\) |
| Solution | \(\bfX^\dagger\bfy=\sum_i\frac{1}{\sigma_i}(\bfu_i^\top\bfy)\bfv_i\) | \(\sum_i\frac{\sigma_i}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\bfv_i\) |
| Unique solution | only when \(\operatorname{rank}(\bfX)=d\) | always |
| Bias | \(\mathbf{0}\) | nonzero, grows with \(\lambda\) |
| Variance | \(\sigma^2\sum_i\sigma_i^{-2}\) | \(\sigma^2\sum_i\sigma_i^2/(\sigma_i^2+\lambda)^2\), shrinks with \(\lambda\) |
| Optimality | best unbiased linear (Gauss–Markov) | smaller MSE for suitable \(\lambda\) (Theorem 3) |
Small singular values make least squares unstable. Ridge regression trades a little bias for a large drop in variance by shrinking the poorly measured directions, and cross-validation chooses how much to shrink.