Multicollinearity, Variance Inflation, Regularized Shrinkage (Ridge & LASSO) & Dimension Reduction

Let $\mathbf{X} = [\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_p] \in \mathbb{R}^{n \times p}$ denote the design matrix of $p$ predictors (or let $\mathbf{X}$ denote the $n \times k$ matrix of centered and standardized non-constant predictors after sweeping out the intercept via the Frisch–Waugh–Lovell theorem, so that $\mathbf{X}^\top \mathbf{X} = \mathbf{R}_{XX}$ is the $k \times k$ sample correlation matrix of the explanatory variables).

Exact (Perfect) Multicollinearity and Non-Identifiability

Exact multicollinearity occurs when there exists a non-zero vector $\mathbf{c} = (c_1, c_2, \dots, c_p)^\top \neq \mathbf{0}$ such that an exact linear combination of the columns of $\mathbf{X}$ vanishes across all $n$ observations in the sample:

Equivalently, $\operatorname{rank}(\mathbf{X}) < p$, which implies $\operatorname{rank}(\mathbf{X}^\top \mathbf{X}) < p$ and $\det(\mathbf{X}^\top \mathbf{X}) = 0$. What happens to the normal equations $(\mathbf{X}^\top \mathbf{X})\hat{\boldsymbol{\beta}} = \mathbf{X}^\top \mathbf{Y}$ when $\mathbf{X}\mathbf{c} = \mathbf{0}$?

1. Uniqueness of the Fitted Projection $\hat{\mathbf{Y}}$: Because the column space $\mathcal{C}(\mathbf{X}) \subset \mathbb{R}^n$ is a well-defined linear subspace of dimension $r = \operatorname{rank}(\mathbf{X}) < p$, the orthogonal projection $\hat{\mathbf{Y}} = \mathbf{P}_{\mathcal{C}(\mathbf{X})}\mathbf{Y}$ of $\mathbf{Y}$ onto $\mathcal{C}(\mathbf{X})$ (and hence the residual vector $\mathbf{e} = \mathbf{Y} - \hat{\mathbf{Y}}$ and $\text{SSE} = \mathbf{e}^\top \mathbf{e}$) remains completely unique! Specifically, for any generalized inverse $(\mathbf{X}^\top \mathbf{X})^-$ satisfying $(\mathbf{X}^\top \mathbf{X})(\mathbf{X}^\top \mathbf{X})^-(\mathbf{X}^\top \mathbf{X}) = \mathbf{X}^\top \mathbf{X}$, the projection matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}^\top \mathbf{X})^-\mathbf{X}^\top$ is invariant to the choice of generalized inverse.

2. Non-Identifiability of the Parameter Vector $\boldsymbol{\beta}$: If $\hat{\boldsymbol{\beta}}_0 = (\mathbf{X}^\top \mathbf{X})^-\mathbf{X}^\top \mathbf{Y}$ is any particular solution to the normal equations, then for every scalar $\alpha \in \mathbb{R}$, the shifted vector $\hat{\boldsymbol{\beta}}(\alpha) = \hat{\boldsymbol{\beta}}_0 + \alpha \mathbf{c}$ generates the exact same fitted values:

Thus, there exists an infinite affine subspace of solutions $\hat{\boldsymbol{\beta}}_0 + \operatorname{Null}(\mathbf{X})$ that achieve the exact same minimum sum of squared residuals $\text{SSE}$. Individual parameters $\beta_j$ associated with nonzero entries of $\mathbf{c}$ are statistically unidentifiable: the data cannot distinguish the separate partial effect of $X_j$ holding other variables constant because $X_j$ cannot vary independently of the other predictors in the sample. (Only linear functions $\boldsymbol{\ell}^\top \boldsymbol{\beta}$ where $\boldsymbol{\ell} \in \mathcal{R}(\mathbf{X}^\top \mathbf{X}) = \operatorname{Null}(\mathbf{X})^\perp$ are estimable).

Near (Approximate) Multicollinearity and Geometric Instability

In practice, exact linear dependence usually arises from deterministic specification errors—such as including all $m$ dummy variables for an $m$-category factor alongside an intercept (the Dummy Variable Trap), including both nominal income, real income, and the price deflator in logarithmic form ($\ln Y_{\text{real}} = \ln Y_{\text{nom}} - \ln P$), or having fewer observations than parameters ($n < p$). Far more pervasive is near multicollinearity, where no exact vector $\mathbf{c} \neq \mathbf{0}$ satisfies $\mathbf{X}\mathbf{c} = \mathbf{0}$, but there exists a unit vector $\mathbf{c}$ ($\|\mathbf{c}\|_2 = 1$) such that:

Geometrically, in $n$-dimensional observation space $\mathbb{R}^n$, the column vectors $\mathbf{x}_1, \dots, \mathbf{x}_p$ span a $p$-dimensional hyperplane $\mathcal{C}(\mathbf{X})$, but two or more vectors point in nearly identical (or nearly opposite) directions, separated by a tiny acute angle $\theta \approx 0$. Just as attempting to balance a flat tabletop on two support legs placed only a millimeter apart results in wild tilting along the perpendicular axis whenever a speck of dust lands on the surface, projecting $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ onto two nearly collinear vectors $\mathbf{x}_1$ and $\mathbf{x}_2$ requires huge, oppositely signed coordinates $\hat{\beta}_1 \gg 0$ and $\hat{\beta}_2 \ll 0$ to represent small orthogonal components of the noise vector $\boldsymbol{\varepsilon}$.

Practical Symptoms of Severe Multicollinearity in Empirical Regression

When near multicollinearity is present, standard OLS output exhibits a characteristic syndrome of diagnostic anomalies:

1. The High-$R^2$ / Few-Significant-$t$-Ratios Paradox: The overall regression achieves a very high coefficient of determination $R^2$ and a highly significant omnibus ANOVA $F$-test (rejecting $H_0: \beta_1 = \dots = \beta_k = 0$ with a tiny $p$-value), yet the individual partial $t$-statistics $t_j = \hat{\beta}_j / \operatorname{SE}(\hat{\beta}_j)$ for the collinear predictors are statistically insignificant! Because the predictors move together, the data confirm with high certainty that the *group* of variables jointly predicts $\mathbf{Y}$, but cannot apportion credit among individual members of the group.

2. Strong Negative Correlation Between Estimated Coefficients: If two standardized predictors $X_1$ and $X_2$ have a strong positive correlation $r_{12} > 0$, their OLS slope estimators $\hat{\beta}_1$ and $\hat{\beta}_2$ exhibit a strong negative sampling correlation $\operatorname{Corr}(\hat{\beta}_1, \hat{\beta}_2) = -r_{12} < 0$. An overestimate of $\beta_1$ in a given sample is automatically compensated by a corresponding underestimate of $\beta_2$ so that the sum $\hat{\beta}_1 X_{i1} + \hat{\beta}_2 X_{i2}$ remains close to $Y_i$.

3. Extreme Sensitivity to Sample Perturbations: Deleting or adding a single observation, or perturbing $Y_i$ in the third decimal place, can cause $\hat{\boldsymbol{\beta}}_{\text{OLS}}$ to swing dramatically in magnitude and even flip algebraic sign.

---

To quantify the exact factor by which multicollinearity inflates the sampling variance of the $j$-th regression coefficient $\hat{\beta}_j$, let us center all variables around their sample means and standardize each predictor so that $\mathbf{z}_j = \frac{\mathbf{x}_j - \bar{x}_j \mathbf{1}_n}{\sqrt{S_{x_j x_j}}}$ has unit length ($\mathbf{z}_j^\top \mathbf{z}_j = 1$). Let $\mathbf{Z} = [\mathbf{z}_1, \dots, \mathbf{z}_k] \in \mathbb{R}^{n \times k}$ be the standardized design matrix, so that $\mathbf{Z}^\top \mathbf{Z} = \mathbf{R}_{XX}$ is the $k \times k$ correlation matrix of the explanatory variables.

Exact Derivation via the Frisch–Waugh–Lovell Theorem and Block Matrix Inversion

Without loss of generality, focus on the $j$-th original predictor $\mathbf{x}_j$ and partition the centered design matrix as $\tilde{\mathbf{X}} = [\tilde{\mathbf{x}}_j, \; \tilde{\mathbf{X}}_{-j}]$, where $\tilde{\mathbf{x}}_j = \mathbf{x}_j - \bar{x}_j \mathbf{1}_n$ is the centered $j$-th predictor (with sum of squares $S_{jj} = \tilde{\mathbf{x}}_j^\top \tilde{\mathbf{x}}_j = \sum_{i=1}^n (X_{ij} - \bar{X}_j)^2$) and $\tilde{\mathbf{X}}_{-j} \in \mathbb{R}^{n \times (k-1)}$ contains the remaining $k - 1$ centered predictors. The centered Gram matrix is partitioned into a $2 \times 2$ block matrix:

By the Schur complement formula for the $(1, 1)$ entry of the inverse of a partitioned symmetric matrix (or equivalently by the Frisch–Waugh–Lovell theorem), the $(j, j)$ diagonal element of $(\tilde{\mathbf{X}}^\top \tilde{\mathbf{X}})^{-1}$ is the reciprocal of the Schur complement of $\tilde{\mathbf{X}}_{-j}^\top \tilde{\mathbf{X}}_{-j}$:

where $\mathbf{H}_{-j} = \tilde{\mathbf{X}}_{-j} (\tilde{\mathbf{X}}_{-j}^\top \tilde{\mathbf{X}}_{-j})^{-1} \tilde{\mathbf{X}}_{-j}^\top$ is the orthogonal projection matrix onto the column space of the other $k - 1$ centered predictors! Now recognize what the quadratic form in the denominator represents:

• In an auxiliary OLS regression of $X_j$ on an intercept and all other $k - 1$ explanatory variables $X_1, \dots, X_{j-1}, X_{j+1}, \dots, X_k$:

the total sum of squares of $X_j$ is $\text{SST}_j = \tilde{\mathbf{x}}_j^\top \tilde{\mathbf{x}}_j = S_{jj}$, the explained sum of squares is $\text{SSR}_j = \tilde{\mathbf{x}}_j^\top \mathbf{H}_{-j} \tilde{\mathbf{x}}_j$, and the residual sum of squares is $\text{SSE}_j = \tilde{\mathbf{x}}_j^\top (\mathbf{I}_n - \mathbf{H}_{-j}) \tilde{\mathbf{x}}_j$.

• Dividing $\text{SSE}_j$ by $\text{SST}_j = S_{jj}$ invokes the coefficient of determination $R_j^2 = \frac{\text{SSR}_j}{\text{SST}_j} = 1 - \frac{\text{SSE}_j}{\text{SST}_j}$ from regressing $X_j$ on the remaining $k - 1$ predictors:

Substituting this exact identity into $\operatorname{Var}(\hat{\beta}_j) = \sigma^2 \left[(\tilde{\mathbf{X}}^\top \tilde{\mathbf{X}})^{-1}\right]_{jj}$ yields the fundamental variance decomposition of multiple regression:

Definition and Interpretation of $\text{VIF}_j$ and Tolerance

The multiplier that inflates $\operatorname{Var}(\hat{\beta}_j)$ relative to an ideal orthogonal design (where $X_j$ is uncorrelated with all other predictors, so $R_j^2 = 0$) is defined by Donald Marquardt (1970) as the Variance Inflation Factor ($\text{VIF}_j$):

Notice two critical mathematical facts:

1. Diagonal of the Inverse Correlation Matrix: $\text{VIF}_j$ is identically equal to the $j$-th diagonal element of the inverse of the $k \times k$ predictor correlation matrix $\mathbf{R}_{XX}^{-1}$.

2. Standard Error Magnification: Because standard errors are square roots of variances, $\operatorname{SE}(\hat{\beta}_j)$ and the width of the $100(1-\alpha)\%$ confidence interval for $\beta_j$ are multiplied by $\sqrt{\text{VIF}_j}$:

• If $R_j^2 = 0.00$ (complete orthogonality): $\text{TOL}_j = 1.00$, $\text{VIF}_j = 1.0$, and $\sqrt{\text{VIF}_j} = 1.00$.
• If $R_j^2 = 0.75$: $\text{TOL}_j = 0.25$, $\text{VIF}_j = 4.0$, and standard errors are doubled ($\sqrt{\text{VIF}_j} = 2.00$).
• If $R_j^2 = 0.80$: $\text{TOL}_j = 0.20$, $\text{VIF}_j = 5.0$, and standard errors are inflated by $\sqrt{5} \approx 2.236$ (a common moderate-collinearity warning threshold).
• If $R_j^2 = 0.90$: $\text{TOL}_j = 0.10$, $\text{VIF}_j = 10.0$, and standard errors are inflated by $\sqrt{10} \approx 3.162$ (the standard severe-collinearity cutoff).
• If $R_j^2 = 0.99$: $\text{TOL}_j = 0.01$, $\text{VIF}_j = 100.0$, and standard errors explode $10$-fold!

---

While Variance Inflation Factors clearly indicate which variables have inflated variances, $\text{VIF}_j$ suffers from two diagnostic limitations: first, it cannot distinguish whether a single near-linear dependency or several distinct coexisting linear dependencies are operating among the predictors; second, if a near-dependency involves the intercept column (for instance, when a continuous predictor such as absolute Kelvin temperature has a tiny coefficient of variation across observations, $\operatorname{SD}(X_j) / \bar{X}_j \ll 1$), centered correlation-based $\text{VIF}_j$ values completely miss the ill-conditioning of the uncentered matrix $\mathbf{X}$.

To resolve both limitations, David Belsley, Edwin Kuh, and Roy Welsch (1980) developed a comprehensive spectral diagnostic framework based on the Singular Value Decomposition (SVD) of the column-equilibrated $n \times p$ design matrix $\mathbf{X}$ (where each column, including the intercept, is scaled to have unit Euclidean norm $\|\mathbf{x}_j\|_2 = 1$).

Singular Value Decomposition and the Spectral Condition Number

Any $n \times p$ matrix $\mathbf{X}$ ($n \ge p$) admits the thin Singular Value Decomposition:

where:

• $\mathbf{U} = [\mathbf{u}_1, \dots, \mathbf{u}_p] \in \mathbb{R}^{n \times p}$ has orthonormal columns ($\mathbf{U}^\top \mathbf{U} = \mathbf{I}_p$) spanning $\mathcal{C}(\mathbf{X})$.
• $\mathbf{V} = [\mathbf{v}_1, \dots, \mathbf{v}_p] \in \mathbb{R}^{p \times p}$ is an orthogonal matrix ($\mathbf{V}^\top \mathbf{V} = \mathbf{V}\mathbf{V}^\top = \mathbf{I}_p$) whose columns $\mathbf{v}_\ell$ are the orthonormal right singular vectors (principal axes / eigenvectors of $\mathbf{X}^\top \mathbf{X}$).
• $\mathbf{D} = \operatorname{diag}(\mu_1, \mu_2, \dots, \mu_p)$ contains the ordered non-negative singular values $\mu_1 \ge \mu_2 \ge \dots \ge \mu_p \ge 0$.

Cross-multiplying $\mathbf{X}^\top \mathbf{X}$ immediately yields its orthogonal eigendecomposition:

Thus, the eigenvalues of $\mathbf{X}^\top \mathbf{X}$ are the squares of the singular values of $\mathbf{X}$ ($\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p > 0$), and the determinant is $\det(\mathbf{X}^\top \mathbf{X}) = \prod_{\ell=1}^p \lambda_\ell = \prod_{\ell=1}^p \mu_\ell^2$. Whenever one or more singular values $\mu_\ell$ approach zero, the corresponding eigenvectors $\mathbf{v}_\ell$ identify the exact linear combinations of columns that are nearly collinear ($\|\mathbf{X}\mathbf{v}_\ell\|_2 = \mu_\ell \approx 0$).

We define the $\ell$-th Condition Index ($\eta_\ell$) and the overall Spectral Condition Number ($\kappa(\mathbf{X})$) as:

By construction, $\eta_\ell \ge 1$, with $\eta_1 = 1$. If the columns of $\mathbf{X}$ are mutually orthonormal, all singular values equal $1$ and $\kappa(\mathbf{X}) = 1$. As collinearity intensifies, $\mu_p \to 0$ and $\kappa(\mathbf{X}) \to \infty$:

• $10 \le \eta_\ell < 30$: Indicates moderate to strong linear dependence associated with the $\ell$-th singular value.
• $\eta_\ell \ge 30$ (equivalently $\lambda_1 / \lambda_\ell \ge 900$): Signals severe multicollinearity; the number of condition indices exceeding $30$ reveals the exact number of distinct near-linear dependencies present in $\mathbf{X}$.

Belsley–Kuh–Welsch Variance Decomposition Proportions ($\pi_{\ell j}$)

Crucially, because $(\mathbf{X}^\top \mathbf{X})^{-1} = \mathbf{V} \mathbf{D}^{-2} \mathbf{V}^\top = \sum_{\ell=1}^p \frac{1}{\mu_\ell^2} \mathbf{v}_\ell \mathbf{v}_\ell^\top$, the sampling variance of the $j$-th regression coefficient $\hat{\beta}_j$ decomposes additively across the $p$ orthogonal principal components:

where $v_{j\ell}$ is the $(j, \ell)$ entry of $\mathbf{V}$. Define $\phi_{j\ell} = \frac{v_{j\ell}^2}{\mu_\ell^2}$ and $\phi_j = \sum_{\ell=1}^p \phi_{j\ell}$. The Variance Decomposition Proportion ($\pi_{\ell j}$) is the fraction of the sampling variance of $\hat{\beta}_j$ attributable to the $\ell$-th singular value $\mu_\ell$:

The Belsley–Kuh–Welsch Diagnostic Rule: A degrading collinearity relation is diagnosed whenever (i) a singular value has a high condition index ($\eta_\ell \ge 30$) and (ii) two or more regression coefficients $\hat{\beta}_j, \hat{\beta}_m$ simultaneously have high variance decomposition proportions ($\pi_{\ell j} > 0.50$ and $\pi_{\ell m} > 0.50$) associated with that same singular value $\mu_\ell$.

---

When the eigenvalues $\lambda_\ell$ of $\mathbf{X}^\top \mathbf{X}$ are near zero, the total mean squared error (expected squared Euclidean distance) of the OLS estimator explodes:

If $\lambda_p = 10^{-4}$, then $1/\lambda_p = 10{,}000$, causing $\hat{\boldsymbol{\beta}}_{\text{OLS}}$ to lie on average very far from the true parameter vector $\boldsymbol{\beta}$, with an inflated squared Euclidean norm $\mathbb{E}[\|\hat{\boldsymbol{\beta}}_{\text{OLS}}\|_2^2] = \|\boldsymbol{\beta}\|_2^2 + \sigma^2 \sum_{\ell=1}^p \frac{1}{\lambda_\ell} \gg \|\boldsymbol{\beta}\|_2^2$.

Formulation of Hoerl–Kennard Ridge Regression (1970)

Drawing on Andrey Tikhonov's (1943) regularization theory for ill-posed inverse problems, Arthur Hoerl and Robert Kennard (1970) proposed adding a strictly positive scalar multiple $\lambda \mathbf{I}_p$ ($\lambda > 0$, often denoted $k > 0$) to the ill-conditioned Gram matrix $\mathbf{X}^\top \mathbf{X}$ before inversion. Equivalently, after centering $\mathbf{Y}$ and standardizing the $p$ predictors so that the intercept $\hat{\beta}_0 = \bar{Y}$ is estimated separately without penalization, Ridge Regression minimizes the $\ell_2$-penalized residual sum of squares:

By Lagrangian duality, this penalized unconstrained optimization is mathematically equivalent to the constrained optimization problem $\min_{\boldsymbol{\beta}} \|\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}\|_2^2$ subject to the spherical Euclidean ball constraint $\|\boldsymbol{\beta}\|_2^2 \le t$ (where $t > 0$ has a one-to-one inverse correspondence with $\lambda > 0$). It also admits a Bayesian Conjugate Normal Interpretation: if the prior distribution on $\boldsymbol{\beta}$ is isotropic Gaussian $\boldsymbol{\beta} \sim \mathcal{N}_p(\mathbf{0}, \tau^2 \mathbf{I}_p)$ and the likelihood is $\mathbf{Y} \mid \boldsymbol{\beta} \sim \mathcal{N}_n(\mathbf{X}\boldsymbol{\beta}, \sigma^2 \mathbf{I}_n)$, then the posterior distribution $\boldsymbol{\beta} \mid \mathbf{Y}$ is Gaussian with posterior mean (and MAP mode) equal to $\hat{\boldsymbol{\beta}}_{\text{ridge}}(\lambda)$ for $\lambda = \sigma^2 / \tau^2$.

Differentiating the Ridge objective $Q_\lambda(\boldsymbol{\beta}) = (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta})^\top (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}) + \lambda \boldsymbol{\beta}^\top \boldsymbol{\beta}$ with respect to $\boldsymbol{\beta}$ and setting the gradient to zero yields the Ridge Normal Equations:

Because the eigenvalues of $\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}_p$ are $\lambda_\ell + \lambda \ge \lambda > 0$, the matrix $\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}_p$ is always strictly positive definite and invertible for every $\lambda > 0$, even if $\operatorname{rank}(\mathbf{X}) < p$ or $p > n$! Its spectral condition number is stabilized to $\kappa\left(\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}_p\right) = \frac{\lambda_1 + \lambda}{\lambda_p + \lambda} \ll \frac{\lambda_1}{\lambda_p}$. Solving for $\hat{\boldsymbol{\beta}}_{\text{ridge}}(\lambda)$ gives:

where $\mathbf{W}_\lambda = (\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}_p)^{-1}(\mathbf{X}^\top \mathbf{X}) = \left(\mathbf{I}_p + \lambda (\mathbf{X}^\top \mathbf{X})^{-1}\right)^{-1}$ is the linear shrinkage matrix.

Spectral Filter Representation, Bias Vector, and Covariance Matrix

Substituting the Singular Value Decomposition $\mathbf{X} = \mathbf{U}\mathbf{D}\mathbf{V}^\top$ (with $\lambda_\ell = \mu_\ell^2$) into the Ridge estimator reveals how Ridge operates along each principal axis $\mathbf{v}_\ell$:

Compared to OLS, which weights each principal direction by $1$, Ridge multiplies the $\ell$-th canonical component by the spectral shrinkage factor $f_\ell(\lambda) = \frac{\lambda_\ell}{\lambda_\ell + \lambda} \in (0, 1)$. Directions with large eigenvalues ($\lambda_\ell \gg \lambda$, high information) experience negligible shrinkage ($f_\ell(\lambda) \approx 1$), whereas collinear directions with tiny eigenvalues ($\lambda_\ell \ll \lambda$, high noise) are heavily damped toward zero! The Effective Degrees of Freedom of the Ridge fit is the trace of the Ridge hat matrix $\mathbf{H}_\lambda = \mathbf{X}(\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}_p)^{-1}\mathbf{X}^\top$:

Now let us compute the exact Bias Vector and Covariance Matrix of $\hat{\boldsymbol{\beta}}_{\text{ridge}}(\lambda) = \mathbf{W}_\lambda \hat{\boldsymbol{\beta}}_{\text{OLS}}$:

1. Expectation and Bias:

2. Variance-Covariance Matrix:

3. Total Mean Squared Error Decomposition: Transforming to canonical coordinates $\boldsymbol{\alpha} = \mathbf{V}^\top \boldsymbol{\beta}$ (where $\|\boldsymbol{\alpha}\|_2^2 = \|\boldsymbol{\beta}\|_2^2$), the total scalar MSE is the sum of the trace of the covariance matrix and the squared Euclidean norm of the bias vector:

Theorem and Proof: The Hoerl–Kennard Existence Theorem

[!IMPORTANT] Theorem (Hoerl and Kennard, 1970): For any full-rank linear regression model with $\sigma^2 > 0$ and bounded true parameter vector $\|\boldsymbol{\beta}\|_2 < \infty$, there always exists a strictly positive regularization parameter $\lambda^* > 0$ such that for all $\lambda \in \left(0, \, \frac{2\sigma^2}{\max_\ell \alpha_\ell^2}\right)$ (and specifically for all $\lambda \in \left(0, \, \frac{\sigma^2}{\max_\ell \alpha_\ell^2}\right)$ in the matrix Loewner ordering), the Ridge estimator achieves a strictly smaller total Mean Squared Error than the OLS estimator:

Proof: Differentiate the variance term $V(\lambda)$ and the squared bias term $B^2(\lambda)$ with respect to $\lambda$ for $\lambda \ge 0$:

Now evaluate both derivatives at the OLS boundary $\lambda = 0^+$:

Adding the two derivatives yields the initial slope of $\operatorname{MSE}(\lambda)$ at $\lambda = 0^+$:

Why does this happen? At $\lambda = 0$, the squared bias $B^2(\lambda) = O(\lambda^2)$ has a zero first derivative (horizontal tangent), whereas the variance $V(\lambda)$ drops with a strictly negative linear slope proportional to $\sum \lambda_\ell^{-2}$! Consequently, moving $\lambda$ slightly above zero trades a second-order $O(\lambda^2)$ increase in squared bias for a first-order $O(\lambda)$ reduction in variance. More specifically, $\operatorname{MSE}'(\lambda) = 2\sum_{\ell=1}^p \frac{\lambda_\ell (\lambda \alpha_\ell^2 - \sigma^2)}{(\lambda_\ell + \lambda)^3} < 0$ is guaranteed for every $\lambda \in \left(0, \, \frac{\sigma^2}{\max_{1 \le \ell \le p} \alpha_\ell^2}\right)$. $\blacksquare$

In empirical practice, visual selection of $\lambda$ is aided by the Ridge Trace Plot, which graphs the $p$ coordinates of $\hat{\boldsymbol{\beta}}_{\text{ridge}}(\lambda)$ against $\lambda \ge 0$: as $\lambda$ increases from $0$, wild oscillations and sign reversals among collinear coefficients rapidly stabilize at a modest value of $\lambda$.

---

Although Ridge Regression dramatically improves estimation accuracy under multicollinearity, the smooth quadratic penalty $\|\boldsymbol{\beta}\|_2^2 = \sum_{j=1}^p \beta_j^2$ shrinks coefficients asymptotically toward zero without ever setting any coefficient exactly equal to zero for finite $\lambda > 0$. Consequently, Ridge retains all $p$ variables in the model and does not perform automatic subset selection.

Formulation of the LASSO and Subgradient KKT Optimality

To achieve simultaneous continuous shrinkage and exact variable selection (sparsity), Robert Tibshirani (1996) introduced the Least Absolute Shrinkage and Selection Operator (LASSO), which replaces the $\ell_2$ penalty with the $\ell_1$ norm $\|\boldsymbol{\beta}\|_1 = \sum_{j=1}^p |\beta_j|$:

Equivalently, LASSO minimizes $\|\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}\|_2^2$ subject to the $\ell_1$ diamond/polytope constraint $\sum_{j=1}^p |\beta_j| \le t$. From a Bayesian standpoint, the $\ell_1$ penalty corresponds to an independent Laplace (Double Exponential) prior $p(\beta_j) = \frac{\lambda}{2\sigma^2}\exp\left(-\frac{\lambda}{\sigma^2}|\beta_j|\right)$, whose sharp cusp at the origin concentrates posterior mode mass directly on sparse coordinate axes.

Geometrically, why does LASSO produce exact zeros while Ridge does not? The residual sum of squares $\|\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}\|_2^2$ has elliptical level contours centered at the unconstrained OLS estimate $\hat{\boldsymbol{\beta}}_{\text{OLS}}$.

• In Ridge Regression, the feasible region $\|\boldsymbol{\beta}\|_2^2 \le t$ is a smooth hypersphere with no corners; the expanding RSS ellipses almost surely touch the sphere at an interior point where all coordinates $\hat{\beta}_j \neq 0$.
• In LASSO, the feasible region $\|\boldsymbol{\beta}\|_1 \le t$ is a cross-polytope (a diamond in $\mathbb{R}^2$, an octahedron in $\mathbb{R}^3$) with sharp vertices and edges lying directly on the coordinate axes where one or more $\beta_j = 0$. Expanding elliptical contours frequently intersect a corner of the diamond first, setting the corresponding coordinate(s) identically to zero!

Closed-Form Derivation of the Soft-Thresholding Operator via Subgradient Calculus

Because the absolute value function $|\beta_j|$ is non-differentiable at $\beta_j = 0$, standard derivative calculus is replaced by convex subdifferential calculus. The subdifferential $\partial |\beta_j|$ of $|\beta_j|$ is the set-valued mapping:

By the Karush–Kuhn–Tucker (KKT) optimality conditions for convex programs, $\hat{\boldsymbol{\beta}}$ minimizes $\frac{1}{2}\|\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}\|_2^2 + \lambda \|\boldsymbol{\beta}\|_1$ if and only if $\mathbf{0} \in -\mathbf{X}^\top(\mathbf{Y} - \mathbf{X}\hat{\boldsymbol{\beta}}) + \lambda \partial \|\hat{\boldsymbol{\beta}}\|_1$, which states coordinate-wise for each $j = 1, \dots, p$:

Suppose the columns of $\mathbf{X}$ are orthonormal ($\mathbf{X}^\top \mathbf{X} = \mathbf{I}_p$), so that the unconstrained OLS estimate is $\hat{\beta}_{j,\text{OLS}} = \mathbf{x}_j^\top \mathbf{Y}$. Then $\mathbf{x}_j^\top(\mathbf{Y} - \mathbf{X}\hat{\boldsymbol{\beta}}) = \hat{\beta}_{j,\text{OLS}} - \hat{\beta}_j$, and the KKT stationarity condition becomes:

We analyze the three exhaustive cases:

1. Case 1 ($\hat{\beta}_j > 0 \implies s_j = +1$): Then $\hat{\beta}_j = \hat{\beta}_{j,\text{OLS}} - \lambda > 0$, which holds if and only if $\hat{\beta}_{j,\text{OLS}} > \lambda$.

2. Case 2 ($\hat{\beta}_j < 0 \implies s_j = -1$): Then $\hat{\beta}_j = \hat{\beta}_{j,\text{OLS}} + \lambda < 0$, which holds if and only if $\hat{\beta}_{j,\text{OLS}} < -\lambda$.

3. Case 3 ($\hat{\beta}_j = 0 \implies s_j \in [-1, +1]$): Then $\hat{\beta}_{j,\text{OLS}} = \lambda s_j \in [-\lambda, +\lambda]$, which holds if and only if $|\hat{\beta}_{j,\text{OLS}}| \le \lambda$.

Combining all three cases yields Donoho and Johnstone's Soft-Thresholding Operator $S_\lambda(\cdot)$:

For general non-orthogonal designs $\mathbf{X}$, Cyclic Coordinate Descent (Friedman, Hastie, and Tibshirani, 2010) updates each coordinate $\hat{\beta}_j$ iteratively by computing the partial residual $\mathbf{r}^{(-j)} = \mathbf{Y} - \sum_{\ell \neq j} \mathbf{x}_\ell \hat{\beta}_\ell$ and applying soft-thresholding: $\hat{\beta}_j \leftarrow \frac{1}{\mathbf{x}_j^\top \mathbf{x}_j} S_\lambda\left(\mathbf{x}_j^\top \mathbf{r}^{(-j)}\right)$. Furthermore, Bradley Efron, Trevor Hastie, Iain Johnstone, and Robert Tibshirani (2004) proved that the entire continuous LASSO regularization path $\{\hat{\boldsymbol{\beta}}_{\text{LASSO}}(\lambda) : \lambda \ge 0\}$ is piecewise linear as a function of $\lambda$! Their Least Angle Regression (LARS) algorithm computes the exact entire LASSO solution path with the same order of computational complexity as a single Ordinary Least Squares fit ($O(np^2)$ operations) by moving along equiangular directions that maintain equal absolute correlation with all currently active predictors until a new variable enters the active set or an existing coefficient crosses zero.

Despite its remarkable geometric elegance, standard LASSO applies the exact same penalty rate $\lambda$ to every coefficient, which introduces a constant shrinkage bias $\pm \lambda$ even for large, genuinely non-zero true parameters $\beta_j$. Moreover, Zhao and Yu (2006) demonstrated that standard LASSO achieves consistent model selection (recovering the exact true support $S_0 = \{j : \beta_j \neq 0\}$ with probability approaching one) if and only if the design matrix satisfies the stringent Irrepresentable Condition, which bounds the regression coefficients of irrelevant predictors on relevant predictors strictly below unity in $\ell_1$ norm—a condition frequently violated under strong multicollinearity. To overcome this limitation and eliminate the constant bias on large coefficients, Hui Zou (2006) introduced the Adaptive LASSO, which assigns data-dependent inverse weights $w_j = 1 / |\hat{\beta}_{j,\text{init}}|^\gamma$ ($\gamma > 0$, where $\hat{\boldsymbol{\beta}}_{\text{init}}$ is a root-$n$ consistent initial estimator such as OLS or Ridge) to the $\ell_1$ penalty $\lambda \sum_{j=1}^p w_j |\beta_j|$. Adaptive LASSO enjoys the celebrated Oracle Property: asymptotically, it selects the exact true subset of non-zero variables with probability tending to one while estimating their non-zero coefficients with the exact same asymptotic normal covariance matrix as if the true sparse model had been known in advance.

Zou and Hastie's Elastic Net (2005)

While LASSO excels at sparse selection, it exhibits three limitations under severe multicollinearity: (a) when $p > n$, LASSO can select at most $n$ variables before saturating; (b) if a group of predictors exhibits very high pairwise correlation, LASSO tends to arbitrarily select one variable from the group and drop the others; and (c) when $n > p$ and correlations are high, Ridge empirically outperforms LASSO in prediction accuracy.

Hui Zou and Trevor Hastie (2005) resolved these issues via the Elastic Net, which combines the $\ell_1$ LASSO penalty and the $\ell_2$ Ridge penalty:

Often parameterized with $\lambda > 0$ and mixing parameter $\alpha \in [0, 1]$ as $\lambda\left[\alpha \|\boldsymbol{\beta}\|_1 + \frac{1-\alpha}{2}\|\boldsymbol{\beta}\|_2^2\right]$, the strictly convex $\ell_2$ term $(\lambda_2/2)\|\boldsymbol{\beta}\|_2^2$ induces a grouping effect: any two standardized predictors with sample correlation $r_{j\ell} \to 1$ are guaranteed to receive nearly identical estimated coefficients, bounded explicitly by the inequality $|\hat{\beta}_j - \hat{\beta}_\ell| \le \frac{\|\mathbf{Y}\|_2}{\lambda_2}\sqrt{2(1 - r_{j\ell})}$, while the sharp corners of the $\ell_1$ term $\lambda_1 \|\boldsymbol{\beta}\|_1$ preserve exact sparsity. Under orthonormal columns ($\mathbf{X}^\top \mathbf{X} = \mathbf{I}_p$), the Elastic Net solution combines soft-thresholding with Ridge scaling:

---

An alternative to continuous penalty shrinkage is **discrete spectral dimension reduction**, which projects the $p$ correlated predictors onto an $m$-dimensional orthogonal subspace ($m < p$) before applying OLS. ### 1. Principal Component Regression (PCR) Using the Singular Value Decomposition $\mathbf{X} = \mathbf{U}\mathbf{D}\mathbf{V}^\top$ of the centered and standardized design matrix (so $\mathbf{X}^\top \mathbf{X} = \mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^\top$ with eigenvalues $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p$), define the $n \times p$ matrix of **Principal Component Scores**: $$\mathbf{Z} = \mathbf{X}\mathbf{V} = [\mathbf{z}_1, \mathbf{z}_2, \dots, \mathbf{z}_p], \qquad \text{where } \mathbf{z}_\ell = \mathbf{X}\mathbf{v}_\ell = \mu_\ell \mathbf{u}_\ell$$ Notice that the principal components $\mathbf{z}_1, \dots, \mathbf{z}_p$ are **mutually orthogonal**, with sums of squares equal to the eigenvalues of $\mathbf{X}^\top \mathbf{X}$: $$\mathbf{Z}^\top \mathbf{Z} = \mathbf{V}^\top (\mathbf{X}^\top \mathbf{X}) \mathbf{V} = \boldsymbol{\Lambda} = \operatorname{diag}(\lambda_1, \lambda_2, \dots, \lambda_p)$$ Rewriting the linear model in terms of the orthogonal scores gives $\mathbf{Y} = \mathbf{X}\mathbf{V}\mathbf{V}^\top \boldsymbol{\beta} + \boldsymbol{\varepsilon} = \mathbf{Z}\boldsymbol{\alpha} + \boldsymbol{\varepsilon}$, where $\boldsymbol{\alpha} = \mathbf{V}^\top \boldsymbol{\beta}$. Because $\mathbf{Z}^\top \mathbf{Z}$ is diagonal, the OLS estimator of the $\ell$-th canonical parameter $\alpha_\ell$ is simply $\hat{\alpha}_\ell = \frac{\mathbf{z}_\ell^\top \mathbf{Y}}{\mathbf{z}_\ell^\top \mathbf{z}_\ell} = \frac{\mathbf{v}_\ell^\top \mathbf{X}^\top \mathbf{Y}}{\lambda_\ell}$ with variance $\operatorname{Var}(\hat{\alpha}_\ell) = \frac{\sigma^2}{\lambda_\ell}$. In **Principal Component Regression (PCR)**, we retain only the first $m < p$ principal components $\mathbf{Z}_m = [\mathbf{z}_1, \dots, \mathbf{z}_m]$ corresponding to the $m$ largest eigenvalues $\lambda_1 \ge \dots \ge \lambda_m$, and hard-threshold the remaining $p - m$ high-variance, low-eigenvalue components to zero ($\hat{\alpha}_{m+1} = \dots = \hat{\alpha}_p = 0$). Transforming back to the original coordinate basis via $\mathbf{V}_m = [\mathbf{v}_1, \dots, \mathbf{v}_m]$ yields the **PCR Estimator**: $$\hat{\boldsymbol{\beta}}_{\text{PCR}}(m) = \mathbf{V}_m \hat{\boldsymbol{\alpha}}_m = \sum_{\ell=1}^m \left(\frac{\mathbf{u}_\ell^\top \mathbf{Y}}{\mu_\ell}\right) \mathbf{v}_\ell = \sum_{\ell=1}^m \frac{1}{\lambda_\ell} \mathbf{v}_\ell \mathbf{v}_\ell^\top \mathbf{X}^\top \mathbf{Y}$$ with covariance matrix $\operatorname{Var}\left(\hat{\boldsymbol{\beta}}_{\text{PCR}}(m)\right) = \sigma^2 \sum_{\ell=1}^m \frac{1}{\lambda_\ell} \mathbf{v}_\ell \mathbf{v}_\ell^\top$ and bias $\operatorname{Bias}\left(\hat{\boldsymbol{\beta}}_{\text{PCR}}(m)\right) = -\sum_{\ell=m+1}^p \alpha_\ell \mathbf{v}_\ell$. While Ridge applies smooth weights $\frac{\lambda_\ell}{\lambda_\ell + \lambda} \in (0, 1)$, PCR applies a discrete step-function filter ($1$ for $\ell \le m$, $0$ for $\ell > m$). Because orthogonal vectors $\mathbf{v}_\ell$ preserve Euclidean norms, the total Mean Squared Error of the $m$-component PCR estimator decomposes cleanly into retained variance plus omitted signal squared: $$\operatorname{MSE}\left(\hat{\boldsymbol{\beta}}_{\text{PCR}}(m)\right) = \sigma^2 \sum_{\ell=1}^m \frac{1}{\lambda_\ell} + \sum_{\ell=m+1}^p \alpha_\ell^2$$ Consequently, discarding the $\ell$-th principal component strictly reduces total Mean Squared Error whenever its omitted signal squared $\alpha_\ell^2$ is smaller than its estimation noise variance $\sigma^2 / \lambda_\ell$. ### 2. Partial Least Squares (PLS) Regression A potential vulnerability of PCR is that the principal directions $\mathbf{v}_1, \dots, \mathbf{v}_m$ are chosen in an **unsupervised** manner solely to maximize predictor variance $\operatorname{Var}(\mathbf{X}\mathbf{v}_\ell) = \lambda_\ell$, without consulting the response $\mathbf{Y}$. As pointed out by Ian Jolliffe (1982), there is no theoretical guarantee that the high-variance directions of $\mathbf{X}$ are the directions most strongly correlated with $\mathbf{Y}$; occasionally, a low-eigenvalue principal component carries substantial predictive power for $\mathbf{Y}$. Herman Wold (1966) developed **Partial Least Squares (PLS)** to construct orthogonal latent components $\mathbf{t}_1, \dots, \mathbf{t}_m$ in a **supervised** manner that balances large predictor variance with high predictor-response correlation. Specifically, at step $\ell$, PLS seeks a unit weight vector $\mathbf{w}_\ell$ maximizing the squared sample covariance between the linear combination $\mathbf{X}^{(\ell-1)}\mathbf{w}_\ell$ and $\mathbf{Y}$, where $\operatorname{Cov}^2(\mathbf{X}\mathbf{w}, \mathbf{Y}) = \operatorname{Var}(\mathbf{X}\mathbf{w}) \cdot \operatorname{Var}(\mathbf{Y}) \cdot \operatorname{Corr}^2(\mathbf{X}\mathbf{w}, \mathbf{Y})$: - For the first PLS direction $\mathbf{w}_1$, solve: $$\mathbf{w}_1 = \arg\max_{\|\mathbf{w}\|_2 = 1} \operatorname{Cov}^2(\mathbf{X}\mathbf{w}, \mathbf{Y}) \quad \Longrightarrow \quad \mathbf{w}_1 = \frac{\mathbf{X}^\top \mathbf{Y}}{\|\mathbf{X}^\top \mathbf{Y}\|_2}$$ - Compute the first score vector $\mathbf{t}_1 = \mathbf{X}\mathbf{w}_1$, regress $\mathbf{Y}$ on $\mathbf{t}_1$ to obtain coefficient $\hat{\theta}_1 = \frac{\mathbf{t}_1^\top \mathbf{Y}}{\mathbf{t}_1^\top \mathbf{t}_1}$, orthogonalize (deflate) each column of $\mathbf{X}$ with respect to $\mathbf{t}_1$ via $\mathbf{X}^{(1)} = \mathbf{X} - \mathbf{t}_1 (\mathbf{t}_1^\top \mathbf{t}_1)^{-1}\mathbf{t}_1^\top \mathbf{X}$, and iterate to extract $m$ mutually orthogonal supervised components spanning the $m$-dimensional **Krylov subspace** $\mathcal{K}_m(\mathbf{X}^\top \mathbf{X}, \mathbf{X}^\top \mathbf{Y}) = \operatorname{span}\left\{\mathbf{X}^\top \mathbf{Y}, \, (\mathbf{X}^\top \mathbf{X})\mathbf{X}^\top \mathbf{Y}, \, \dots, \, (\mathbf{X}^\top \mathbf{X})^{m-1}\mathbf{X}^\top \mathbf{Y}\right\}$. Because PLS explicitly targets $\mathbf{X}^\top \mathbf{Y}$, it typically achieves comparable in-sample fit to PCR using fewer components ($m_{\text{PLS}} \le m_{\text{PCR}}$), though its non-linear dependence on $\mathbf{Y}$ can slightly increase variance relative to PCR. ### 3. Selecting the Regularization Parameter: $K$-Fold CV, Exact LOO Shortcut, and GCV Because $\text{SSE}$ on the training data is minimized at $\lambda = 0$ (OLS), the regularization parameter $\lambda$ (or $m$) must be tuned using out-of-sample predictive error criteria: - **$K$-Fold Cross-Validation**: Partition the $n$ observations into $K$ disjoint folds $F_1, \dots, F_K$ (typically $K = 5$ or $10$), fit the regularized model on the remaining $n - |F_k|$ observations for a grid of $\lambda$ values, and minimize $\operatorname{CV}(\lambda) = \frac{1}{n}\sum_{k=1}^K \sum_{i \in F_k} \left(Y_i - \mathbf{x}_i^\top \hat{\boldsymbol{\beta}}^{(-k)}(\lambda)\right)^2$. - **Exact Leave-One-Out (PRESS) Shortcut for Ridge**: For any linear smoother $\hat{\mathbf{Y}}(\lambda) = \mathbf{H}_\lambda \mathbf{Y}$, the $n$-fold Leave-One-Out cross-validation error can be computed from a **single** fit without refitting the model $n$ times: $$\operatorname{LOOCV}(\lambda) = \frac{1}{n}\sum_{i=1}^n \left( Y_i - \mathbf{x}_i^\top \hat{\boldsymbol{\beta}}^{(-i)}(\lambda) \right)^2 = \frac{1}{n}\sum_{i=1}^n \left( \frac{Y_i - \hat{Y}_i(\lambda)}{1 - (\mathbf{H}_\lambda)_{ii}} \right)^2$$ - **Golub–Heath–Wahba (1979) Generalized Cross-Validation (GCV)**: Replaces the individual leverage diagonals $(\mathbf{H}_\lambda)_{ii}$ with their rotation-invariant average $\frac{1}{n}\operatorname{tr}(\mathbf{H}_\lambda) = \frac{1}{n}\operatorname{df}(\lambda)$: $$\operatorname{GCV}(\lambda) = \frac{\frac{1}{n}\left\|\mathbf{Y} - \mathbf{H}_\lambda \mathbf{Y}\right\|_2^2}{\left(1 - \frac{1}{n}\operatorname{tr}(\mathbf{H}_\lambda)\right)^2} = \frac{n \cdot \text{SSE}(\lambda)}{\left(n - \sum_{\ell=1}^p \frac{\lambda_\ell}{\lambda_\ell + \lambda}\right)^2}$$ ---