How does ridge regression reduce the variance of the estimates of $\beta$

Question

In the scikit-learn library, Ridge class, there is a note that reads: "Regularization improves the conditioning of the problem and reduces the variance of the estimates."

Given the expression of beta for ridge regression:

$\beta = (X^T X + \lambda I_d)^{-1} X^T y$

How does the addition of $\lambda I_d$ to the moments matrix improve the condition number of the inverse? And how exactly does it reduce the variance of the estimates (mathematically)?

Answer 1

The additional term $\lambda I_d$ can be interpreted as adding some bias and it is the bias that reduces the variance.

1. The addition of $\lambda I_d$ is equivalent to adding a penalty term to the least squares optimization problem $$\min_{\beta}: \overbrace{ (y-X\beta )^T(y-X\beta ) }^{\text{least squares}}+ \overbrace{\lambda (\beta^T\beta)}^{\text{additional penalty}}$$

2. The penalty term is equivalent to constraining the parameters to a particular maximum size.

   $$\min_{\beta}: { (y-X\beta )^T(y-X\beta ) } \qquad \text{and} \qquad \beta^T\beta < c$$

3. Constraining the parameters will reduce the variance of the parameters.

Answer 2

$\DeclareMathOperator{\rid}{\hat{\boldsymbol\beta}_\text{ridge}}\DeclareMathOperator{\ols}{\hat{\boldsymbol\beta}}\DeclareMathOperator{\bias}{\hat{\boldsymbol\beta^\ast}} \DeclareMathOperator{\tr} {trace}\DeclareMathOperator{\xx}{\mathbf X^\mathsf T\mathbf X}$

Condition number (cf. $\rm [I],$ chapter $4,$ p. $189$) $\kappa:\mathcal M\to \mathbb R^{\geq 0}$ is defined as

$$\kappa(\mathbf A) := \Vert \mathbf A \Vert \left\Vert \mathbf A^{-1}\right\Vert,\tag 1$$ where $\Vert \mathbf A\Vert := \displaystyle\max_{\Vert \mathbf x\Vert = 1} \Vert\mathbf{Ax} \Vert.$ Corresponding to $\Vert \mathbf x\Vert^2 =\sum_{i=1}^n |x_i|^2, $ the matrix norm subordinate to it can be shown (cf.$\rm [II],$ chapter $1,$ p. $60$) to be $\Vert \mathbf A \Vert= \sqrt{\lambda_1}, $ the largest eigenvalue of $\mathbf A^\mathsf T\mathbf A; $ if $\bf A$ is symmteric, then $$ \kappa(\mathbf A) =\frac{\lambda_\max(\mathbf A) }{\lambda_\min(\mathbf A)},\tag{1.I}$$

Eigenvalues of $\xx+ k\mathbf I$ are of the form $\lambda_i+ k, ~\lambda_i$ being an eigenvalue of $\xx;$ then from $\rm(1.I), $ $$\kappa(\mathbf A) =\frac{\lambda_\max + k }{\lambda_\min + k}.\tag 2$$ Since $k> 0,$ one is able to circumvent the scenario of dividing the numerator by a very small number $\lambda_\min, $ i.e. at least there would be $k:$ this would be alleviating from what would have been a very large condition number.

The problem with $\ols$ is that, albeit being unbiased, the variance is large and the point estimate is unstable. One tradeoff could be that of allowing biased estimators, say $\bias,$ in that, as echoed in Sextus Empiricus' post, the variance of $\bias$ can be lowered such that the mean square error of $\bias$ is smaller than the variance of $\ols$ (cf.$\rm[III]$, chapter $9,$ p. $305$).

If $L_1^2:= \left(\ols-\boldsymbol\beta\right) ^\mathsf T\left(\ols-\boldsymbol\beta\right), $ then since $\mathbb E[\boldsymbol\varepsilon ]=\mathbf 0,~\mathbb E\left[\boldsymbol\varepsilon\boldsymbol\varepsilon^\mathsf T\right]=\sigma^2\mathbf I, $ \begin{align}\mathbb E\left[ L^2_1\right]&= \sigma^2\tr\left(\xx\right)^{-1}\\&= \sigma^2\sum_{i=1}^p\frac1{\lambda_i},\tag 3\end{align} $\lambda_i$ being the eigenvalues of $\xx.$

Now, $$\rid= \underbrace{\left[\xx + k\mathbf I\right]^{-1}}_{:=\mathbf W}\mathbf X^\mathsf T\mathbf y; \tag{4.I}$$ equivalently $$\rid =\underbrace{\left[\mathbf I +k\left(\xx\right)^{-1}\right]^{-1}}_{:=\mathbf Z}\ols.\tag{4.II}$$ As $\mathbf Z=\mathbf W\xx, $ $$\mathbf Z= \mathbf I-k\mathbf W. \tag 5$$

If $L^2(k):= \left(\rid-\boldsymbol\beta\right)^\mathsf T \left(\rid-\boldsymbol\beta\right), $ \begin{align}\mathbb E\left[L^2(k)\right]&= \mathbb E\left[\left(\ols-\boldsymbol\beta\right)^\mathsf T \mathbf Z^\mathsf T\mathbf Z\left(\ols-\boldsymbol\beta\right)\right]+ \left(\mathbf Z{\boldsymbol\beta}-\boldsymbol\beta\right)^\mathsf T \left(\mathbf Z{\boldsymbol\beta}-\boldsymbol\beta\right)\\ &= \mathbb E\left[{\boldsymbol\varepsilon}^\mathsf T\mathbf X\left(\xx\right)^{-1}\mathbf Z^\mathsf T\mathbf Z\left(\xx\right)^{-1}\mathbf X^\mathsf T\boldsymbol\varepsilon\right]+ \left(\mathbf Z{\boldsymbol\beta}-\boldsymbol\beta\right)^\mathsf T \left(\mathbf Z{\boldsymbol\beta}-\boldsymbol\beta\right)\\&= \sigma^2\tr\left[\left(\xx\right)^{-1}\mathbf Z^\mathsf T\mathbf Z\right]+ \boldsymbol\beta^\mathsf T\left(\mathbf Z-\mathbf I\right)^\mathsf T \left(\mathbf Z-\mathbf I\right) {\boldsymbol\beta}\\ &\overset{(5)}{=} \sigma^2\left[\tr\mathbf W-k\tr\mathbf W^2\right]+ k^2\boldsymbol\beta^\mathsf T\mathbf W^{2}\boldsymbol\beta\\ &= \underbrace{\sigma^2\sum_{i=1}^p \frac{\lambda_i}{(\lambda_i + k)^2}}_{\gamma_1(k) }+ \underbrace{k^2\sum_{i=1}^p \frac{\alpha_i^2}{(\lambda_i + k)^2}}_{\gamma_2(k)};\tag 6 \end{align} where $\boldsymbol\alpha = \mathbf P\boldsymbol\beta,~\mathbf P $ being the orthogonal matrix such that $\xx = \mathbf{ P\Lambda P}^\mathsf T, ~\mathbf\Lambda :=\operatorname{diag}(\lambda_i).$

Observation $1.$ $\gamma_1(k) ,$ the variance, is decreasing function of $k$ and $\gamma^\prime_1(k) \to -\infty$ as $k\to +0, ~\lambda_\min\to 0.$

\begin{align}\lim_{k\to +0}\gamma^\prime_1(k)&= \lim_{k\to +0}-2\sigma^2\sum_{i=1}^p\frac{\lambda_i}{(\lambda_i + k)^3} \\ &= -2\sigma^2\sum_{i=1}^p\frac{ 1}{\lambda_i^2}.\tag 7\end{align}

Observation $2.$ $\gamma_2(k) ,$ the squared bias is increasing function of $k$ and $\gamma^\prime_2(k) \to 0$ as $k\to +0.$

\begin{align}\lim_{k\to +0}\gamma^\prime_2(k)&= \lim_{k\to +0}2k\sum_{i=1}^p\frac{\lambda_i\alpha_i^2}{(\lambda_i + k)^3} \\ &= 0.\tag 8\end{align}

These two observations indicate there are "admissible" values or $k$ for which $\operatorname{MSE}(\rid) $ is lesser than the variance of $\ols.$

Theorem 1. (cf. $\rm[IV]$) There exists $k > 0$ such that $\mathbb E\left[L_1^2(k)\right] < \mathbb E\left[L_1^2(0)\right].$

Proof. Note $\gamma_1(0)= \sigma^2\sum_{i=1}^p\frac1{\lambda_i}, ~\gamma_2(0) = 0.$

Now, \begin{align}\frac{\mathrm d}{\mathrm dk}\mathbb E\left[L_1^2(k)\right] &= -2\sigma^2\sum_{i=1}^p\frac{\lambda_i}{(\lambda_i + k)^3}+ 2k\sum_{i=1}^p\frac{\lambda_i\alpha_i^2}{(\lambda_i + k)^3}; \tag 9\end{align} based on both the observations, it suffices to show that there exists $k>0$ such that $\frac{\mathrm d}{\mathrm dk}\mathbb E\left[L_1^2(k)\right] <0.$ From $(9), $ the value of $k$ should be $$k < \frac{\sigma^2}{\alpha^2_\max}.\tag{10}$$ $\square$

It must be reiterated again for the sake of gravity, also as mentioned in Sextus Empiricus' comment, the agenda is not to solely decrease the variance but rather lower the mean square error of $\rid$ than the variance of the least square estimator $\ols.$ This is done by introducing a little bias and substantially reducing the variance for certain $k>0$ as shown above.

References:

$[\rm I]$ Algebraic Eigenvalue Problem, J. H. Wilkinson, Oxford University Press, $1965.$

$[\rm II]$ Computational Methods of Linear Algebra, V. N. Faddeeva, Dover Publications, $1959.$

$\rm[III]$ Introduction to Linear Regression Analysis, Douglas C. Montgomery, Elizabeth A. Peck, G. Geoffrey Vining, John Wiley & Sons, $2012.$

$\rm [IV]$ Ridge Regression: Biased Estimation for Nonorthogonal Problems, Arthur E. Hoerl, Robert W. Kennard, Technometrics $42,$ no. $1~ (2000): ~80–86. $ https://doi.org/10.2307/1271436.