Ridge Shrinkage in Ill-Conditioned Design Matrices

Correlated columns make least squares produce coefficients of large magnitude and opposite sign that cancel in the prediction. The fit is fine on the training data and unstable everywhere else. Ridge regression replaces the unconstrained problem with a penalised one and, in the process, converts a numerically fragile inversion into a stable one. Both effects come from a single algebraic operation: adding $\lambda$ to every singular value of the design before inverting.

The penalty rewritten in the SVD basis

The ridge objective is

$$\lVert y - X\beta \rVert^2 + \lambda \lVert \beta \rVert^2,$$

with a stationary condition $(X^\top X + \lambda I)\beta = X^\top y$. Diagonalising the design by its SVD $X = U\Sigma V^\top$ turns the Gram matrix into $V\Sigma^2 V^\top$, so the solution is

$$\hat\beta = V \operatorname{diag}!\left(\frac{s_i}{s_i^2 + \lambda}\right) U^\top y .$$

Compare with the least-squares solution, which has the coefficients $1/s_i$ where ridge has $s_i/(s_i^2+\lambda)$. In the orthonormal basis of singular vectors, ridge multiplies the $i$-th component of the least-squares solution by the factor

$$d_i = \frac{s_i^2}{s_i^2 + \lambda},$$

a number between 0 and 1 that approaches 1 when $s_i^2 \gg \lambda$ and 0 when $s_i^2 \ll \lambda$. Large-variance directions are left nearly untouched; low-variance directions, meaning the ones where the data barely constrains the answer and the least-squares coefficient is dominated by noise, are squashed toward zero. That is the geometric content of the penalty: it is not a uniform shrink toward the origin but a selective one, applied in decreasing order of how well the data determines each direction.

The conditioning follows. The smallest singular value of the penalised matrix is at least $\sqrt{\lambda}$ in the sense that $(X^\top X + \lambda I)$ has eigenvalue $s_{\min}^2 + \lambda$ in the worst direction, so the condition number of the ridge problem is bounded by $(s_{\max}^2 + \lambda)/\lambda$ regardless of how small $s_{\min}$ is. With $\lambda$ fixed and positive, the problem is invertible even when $X$ has fewer rows than columns. That is the property that makes ridge work in the $p > n$ regime, where ordinary least squares has no unique solution at all: an infinite family of coefficient vectors interpolates the data, and ridge picks the one of smallest norm, the same member the pseudo-inverse would select, but by a smoother route.

Bias, variance, and the effective degrees of freedom

Because the estimator is linear in $y$, $\hat\beta_\lambda = (X^\top X + \lambda I)^{-1}X^\top y$, it has a closed-form risk. Writing the true coefficient as $\beta$ and the noise variance as $\sigma^2$, the mean squared error of the fitted vector decomposes into

$$\text{bias}^2 + \sigma^2\sum_i d_i^2 + \sigma^2\sum_i (1-d_i)^2,$$

where the middle term is the variance and the last is the squared bias induced by shrinking the true signal. Ordinary least squares corresponds to $\lambda = 0$ and $d_i = 1$, which leaves the variance term at $\sigma^2 p$ and the bias at zero. Ridge reduces the variance sum monotonically as $\lambda$ grows and increases the bias. There is always a $\lambda$ at which the total is smaller than the least-squares risk, and in the orthogonal case where all $d_i$ are equal, the optimal amount has a closed form proportional to $p\sigma^2/\lVert\beta\rVert^2$, the shrinkage is stronger when there are many parameters, when the noise is large, and when the signal is weak.

The quantity $\sum_i d_i = \operatorname{tr}!\big(X(X^\top X+\lambda I)^{-1}X^\top\big)$ is called the effective degrees of freedom. It interpolates between 0 as $\lambda \to \infty$ and $\operatorname{rank}(X)$ at $\lambda = 0$. For a design with $p$ columns, a ridge fit with a moderate penalty might carry an effective dimension of 3.7 rather than 8, and that number, not $p$, is what belongs in a corrected information criterion, since the penalty term in AIC-type criteria is a function of the model’s flexibility rather than its nominal parameter count.

Why the intercept stays out of the penalty

Shrinking the intercept toward zero means shrinking toward a prediction of zero on the linear scale, which is a different model from shrinking the slopes toward zero. For standardised features, penalising the intercept shifts predictions away from the response mean in direct proportion to the penalty. The usual fix is to centre both $X$ and $y$ before fitting, drop the intercept from the penalty, and recover it as $\bar y$ when the centred slopes are zero. Centring changes nothing about the slopes when the penalty is applied to them, because the penalty is on the deviation from the origin in a centred coordinate system.

Centring also changes the numerical problem for the better. Raw features measured in different units can make $\lambda$ meaningless as a single number: a penalty large enough to matter for a feature measured in dollars is negligible for one measured in millions of dollars, so the shrinkage on each coefficient depends on its arbitrary scale. Standardising to unit variance makes $\lambda$ comparable across coefficients, and combined with centring reduces the correlation between the intercept column and the rest of the design to zero, which removes the largest source of conditioning trouble in raw count and price data.

import numpy as np

def ridge_standardised(X, y, lam):
    xm, ym = X.mean(axis=0), y.mean()
    Xc, yc = X - xm, y - ym
    scale = Xc.std(axis=0, ddof=0)
    scale[scale == 0.0] = 1.0            # constant column: avoid divide by zero
    Z = Xc / scale
    # penalise in standardised coordinates, then map back
    beta_std = np.linalg.solve(Z.T @ Z + lam * np.eye(Z.shape[1]), Z.T @ yc)
    beta = beta_std / scale
    intercept = ym - xm @ beta
    return intercept, beta

Mapping the coefficients back through the scaling is mandatory; a coefficient reported in standardised units is not the coefficient on the original feature, and the two differ by the feature’s standard deviation. Dividing by scale after solving is the correct direction, since a unit change in the original feature is a $1/\sigma_j$ change in the standardised one.

The path is monotone but the coefficients are not

The ridge solution varies continuously with $\lambda$, and the ridge trace (each coefficient plotted against $\lambda$) has a shape worth understanding. Every coefficient converges to zero as $\lambda \to \infty$, but not necessarily in the same direction: a coefficient that is small and positive at $\lambda = 0$ can grow in magnitude before shrinking, because the shrinkage redistributes weight among correlated columns. This is why the trace is used as a diagnostic rather than a control: coefficients that cross or change sign along the path indicate that the corresponding features are nearly collinear and their individual values carry little information.

Cross-validation over $\lambda$ is the usual selection rule. The leave-one-out error for a linear smoother has a closed form,

$$\text{LOO} = \frac{1}{n}\sum_i \left(\frac{y_i - \hat y_i}{1 - h_{ii}}\right)^2,$$

with $h_{ii}$ the leverage of the ridge fit, so the whole path can be evaluated without refitting $n$ times. Generalised cross-validation replaces $\sum h_{ii}$ by its trace, which costs one matrix trace rather than $n$ diagonal entries. The selection is noisy when $n$ is small: the standard error of the cross-validation curve is often wide enough that a broad plateau of $\lambda$ values is indistinguishable, and any point in the plateau behaves similarly.

import numpy as np

def loo_gcv(X, y, lams):
    U, s, Vt = np.linalg.svd(X, full_matrices=False)
    Uty = U.T @ y
    out = []
    for lam in lams:
        d = s**2 / (s**2 + lam)                 # shrinkage factors
        fit = U @ (d * Uty)
        h = (U**2) @ d                          # leverage per observation
        loo = np.mean(((y - fit) / (1 - h))**2)
        gcv = np.mean((y - fit)**2) / (1 - d.sum() / len(y))**2
        out.append((lam, loo, gcv))
    return out

Both criteria need $h_{ii} < 1$; if the penalty is so small that an observation is interpolated, the denominator vanishes and the error estimate is unbounded. That degenerate case only arises when a row has no counterpart in the design, which is a data problem rather than a tuning problem.

What ridge cannot do

Ridge never sets a coefficient exactly to zero, so it cannot select features. Every column stays in the model with a nonzero weight, and for a model with hundreds of thousands of sparse categorical indicators this means dense storage and dense arithmetic throughout. The $\ell_1$ penalty achieves selection by making the objective non-differentiable at the origin, which changes the whole numerical apparatus: no closed form, no single SVD, and a solution path that is piecewise linear in $\lambda$ rather than smooth.

Ridge also assumes the penalty is appropriate to the feature scaling, which is why standardisation is part of the method rather than a preprocessing preference. And it does not repair the estimand: what ends up being estimated is a biased version of $\beta$, so coefficients are not interpretable as effects. They are predictions-oriented weights, useful for forecasting and for stabilising a design, not for reading off marginal associations. The same trade appears in the logistic setting, where the penalty is also what makes a separable problem well posed, as described in logistic regression as conditional likelihood.

Finally, the penalty interacts with the loss. It shrinks toward zero, which is meaningful on a linear predictor scale but arbitrary on a probability scale, so the fitted coefficients of a penalised logistic model are not shrunken toward the base rate but toward even odds. For a GLM with a log link the penalty therefore shrinks toward a relative effect of one, treated in link functions and Poisson regression, and the choice of link silently decides what “small coefficient” means.