Ridge Regression

In the multiple linear regression entry we saw least squares get confused: when predictors carry almost the same information, the data no longer decide how to share the effect between them, and the coefficients swing wildly from one dataset to the next.

Ridge regression is the simplest and oldest cure. Its idea fits in one sentence: among all the coefficient vectors that fit the data well, prefer the ones that are small. That modest preference introduces a little bias and, in exchange, removes a lot of instability. It is the first example of regularization, an idea that runs through modern chemometrics and machine learning.

In 1970, Arthur E. Hoerl and Robert W. Kennard, two statisticians at DuPont in Wilmington, Delaware, published a pair of papers in Technometrics. The first, Ridge Regression: Biased Estimation for Nonorthogonal Problems [1], proposed adding a small positive constant to the diagonal of ${\mathbf{X}}^T\mathbf{X}$ before inverting it. The matrix becomes better conditioned and the coefficients become biased but far more stable. The second, Ridge Regression: Applications to Nonorthogonal Problems [11], showed the method at work on real data.

The name came from Hoerl's earlier work on ridge analysis, a technique for exploring response surfaces that he introduced in 1959 [12] and applied to regression problems in the early 1960s [2]. Ridge analysis follows the best attainable response on spheres of increasing radius around a centre point. Ridge regression does the same thing with the residual sum of squares: for every sphere of a given size around the origin, it finds the best-fitting coefficients on it. We will see that picture in a moment.

Hoerl and Kennard were not the first to have the core idea. In the Soviet Union, Andrey N. Tikhonov had studied the stability of inverse problems since the 1940s [3] and formalized his regularization method in 1963 [4]. In physics and applied mathematics the same technique is called Tikhonov regularization. Two communities, two names, one idea: whenever an inverse problem is ill-conditioned, add a penalty that keeps the solution reasonable.

Their most important message was almost philosophical: a biased estimator can beat an unbiased one when the goal is to be close to the truth on average. That idea, the bias-variance trade-off, is one of the pillars of statistical learning.

Suppose you predict a concentration from the absorbances at three neighbouring wavelengths, say 500, 501 and 502 nm. After centring, the matrix ${\mathbf{X}}^T\mathbf{X}$ might look like this:

Its eigenvalues are about 299.4, 0.5 and 0.1. One direction in coefficient space (all three coefficients moving together) is very well determined by the data; the other two are barely determined at all. The ratio between the largest and smallest eigenvalue, about 3000, is the condition number of ${\mathbf{X}}^T\mathbf{X}$. When the matrix is inverted, the small eigenvalues become large factors ($1/0.1=10$, compared with $1/299.4\approx 0.003$), and the variance of the coefficients along those directions is about 3000 times larger than along the well-determined one.

The result is the familiar symptom: coefficients such as $+41$, $-73$ and $+43$ that almost cancel each other. Add a single sample and they may become $+25$, $-41$ and $+27$. Their sum, the well-determined part, stays at 11, so the fitted values barely change, but the model has learned nothing trustworthy about the individual wavelengths, and predictions for samples that do not follow the calibration pattern exactly can be badly off.

Instead of minimizing only the sum of squared residuals, ridge regression minimizes:

The first term asks the model to fit the data. The second, the L2 penalty (also called the ridge or Tikhonov penalty, and weight decay in neural networks), asks the coefficients to stay small. The tuning parameter $\lambda \ge 0$ sets the balance between the two wishes.

Taking the gradient and setting it to zero, exactly as for the normal equations:

which rearranges into a slightly modified set of normal equations:

The only change from ordinary least squares is the term $\lambda \mathbf{I}$, a constant added to the diagonal. Adding $\lambda$ to the diagonal raises every eigenvalue by $\lambda$. In our example, $\lambda =1$ turns the eigenvalues into 300.4, 1.5 and 1.1, and the condition number drops from about 3000 to about 270. And since every eigenvalue is now at least $\lambda$, the matrix can be inverted even when ${\mathbf{X}}^T\mathbf{X}$ is singular, for example when there are more wavelengths than samples.

At the two extremes: $\lambda =0$ gives back ordinary least squares, and $\lambda \to \infty$ shrinks every slope to zero, so the model predicts $\bar{y}$ for every sample. The art is in choosing something in between.

There is an equivalent way to write the ridge problem that makes its behaviour visible. Minimizing the penalized sum is the same as minimizing the ordinary sum of squares subject to a budget on the size of the coefficients:

Each $\lambda$ corresponds to some budget $t$: a larger penalty means a smaller budget. With two coefficients, the budget is a disc centred at the origin. The sum of squares, as a function of the coefficients, forms the bowl we met in gradient descent; its contour lines are ellipses centred at the least squares estimate. The ridge solution is the point where the smallest ellipse touches the disc.

Play with it. When the least squares point already lies inside the disc, the constraint does nothing. When it lies outside, the solution slides to the edge of the disc, closer to the origin. Increase the correlation and the ellipses become long and thin: that is collinearity seen from coefficient space, a long valley along which the fit hardly changes. The disc cuts that valley short. Switch to the diamond, or tick "Show both shapes", to get a preview of the lasso.

The singular value decomposition, $\mathbf{X}=\mathbf{UD}{\mathbf{V}}^T$, shows precisely what the penalty does [7]. Here the columns ${\mathbf{v}}_j$ of $\mathbf{V}$ are orthogonal directions in coefficient space, the $d_j$ are the singular values (how much the data spread along each direction), and the ${\mathbf{u}}_j$ are the corresponding directions in sample space. The least squares solution can be written as:

A tiny $d_j$ means a huge factor $1/d_j$: whatever noise sits in that direction is multiplied enormously. Ridge regression replaces the factor:

The shrinkage factor is close to 1 for directions where the data are rich ($d_j^2\gg \lambda$) and close to 0 where they are poor ($d_j^2\ll \lambda$). Ridge regression therefore leaves the well-determined directions almost untouched and quietly switches off the poorly determined ones, which are exactly the directions where least squares was amplifying noise.

The sum of the shrinkage factors is a useful number too. It is called the effective degrees of freedom of the fit:

When $\lambda =0$ it equals the rank of $\mathbf{X}$ (the number of non-zero singular values, which is $p$ when there are more samples than variables), and it decreases smoothly towards 0 as $\lambda$ grows: ridge regression behaves like a model with a non-integer, adjustable number of parameters. Keep this formula in mind when you read about principal component regression, which uses factors that are exactly 1 or exactly 0.

The Gauss-Markov theorem says that ordinary least squares has the smallest variance among unbiased linear estimators. Ridge regression is biased: on average, its coefficients are pulled towards zero. Why would anyone accept that?

Because the error we care about has two parts. The mean squared error of an estimator is its squared bias plus its variance:

Least squares has zero bias but, with collinear predictors, an enormous variance. A small penalty adds a little bias and removes a lot of variance. Hoerl and Kennard proved that, for any problem, there is always some $\lambda >0$ for which ridge regression has a smaller mean squared error than least squares [1]. The catch is that the best $\lambda$ depends on the unknown true coefficients, so it has to be estimated from the data.

Hoerl and Kennard suggested looking at the whole family of solutions at once: plot every coefficient as a function of $\lambda$. This plot is the ridge trace. The demo below computes it for a simulated near-infrared calibration with ten wavelengths, together with the calibration error and the cross-validated error for every value of $\lambda$.

Read it from left to right. With almost no penalty you see the least squares coefficients: large, with alternating signs, the typical fingerprint of collinearity. As $\lambda$ increases they settle down, and eventually they all slide towards zero. Notice also the two error curves: the calibration error can only grow with the penalty, because ridge deliberately fits the calibration data a little worse. The cross-validated error first decreases and then increases. The bottom of that curve is the best compromise for prediction. The hollow circle marks a slightly stronger penalty, chosen by the one-standard-error rule described in the next section.

The most widely used approach, and the one in the demo:

Remember to standardize inside each fold, using only the training samples, so that the left-out samples really are unseen.

Golub, Heath and Wahba proposed a fast approximation to leave-one-out cross-validation that needs a single fit per $\lambda$ [5]:

${\mathbf{H}}_{\lambda }$ is the ridge hat matrix (it turns $\mathbf{y}$ into $\hat{\mathbf{y}}$), and its trace is precisely the effective degrees of freedom $\text{df}(\lambda )$ introduced above. With an intercept fitted, add 1 to the trace, since the intercept uses one more degree of freedom.

Plot the size of the solution, $\log \Vert {\hat{\mathbf{\beta }}}_{\lambda }\Vert$, against the size of the residual, $\log \Vert \mathbf{y}-\mathbf{X}{\hat{\mathbf{\beta }}}_{\lambda }\Vert$, for many values of $\lambda$. The curve usually looks like an L. Its corner marks the point where shrinking the coefficients further starts to cost a lot of fit. The L-curve is popular in inverse problems and spectroscopic deconvolution [6].

Ridge regression also has a neat probabilistic interpretation. Suppose, before seeing the data, you believe that each coefficient is probably small, and you express that belief as a Gaussian prior centred at zero:

With Gaussian measurement noise of variance ${\sigma }^2$, the most probable coefficients after seeing the data (the posterior mode, which here is also the posterior mean) are exactly the ridge solution:

So $\lambda$ is the ratio between the noise variance and the prior variance of the coefficients. A large $\lambda$ means "I strongly believe the coefficients are small"; a small $\lambda$ lets the data speak. When collinearity leaves the data ambiguous about some directions, the prior fills the gap. This extends the maximum likelihood reasoning we followed in your first step into chemometrics with one extra ingredient, the prior. The result is no longer a maximum likelihood estimate but a maximum a posteriori (MAP) estimate: the coefficients that are most probable once the data and the prior belief are combined. Replace the Gaussian prior by a Laplace prior and you get the lasso.

import numpy as np
from sklearn.linear_model import Ridge, RidgeCV, LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import KFold, cross_val_score
import matplotlib.pyplot as plt

# --- Simulate spectra: analyte band + overlapping interferent + noise ---
rng = np.random.default_rng(42)
n_samples, n_wavelengths = 40, 30
wavelengths = np.linspace(450, 700, n_wavelengths)
analyte = np.exp(-((wavelengths - 550) ** 2) / 3000)
interferent = np.exp(-((wavelengths - 600) ** 2) / 4000)
conc = rng.uniform(0.5, 10.0, n_samples)            # analyte concentration
c_int = rng.uniform(0.5, 10.0, n_samples)           # interferent, unrelated to y
X = (np.outer(conc, analyte) + np.outer(c_int, interferent)
     + rng.normal(0, 0.05, (n_samples, n_wavelengths)))
y = conc + rng.normal(0, 0.2, n_samples)            # reference method error

# Standardize predictors (the penalty must treat all variables equally)
X_scaled = StandardScaler().fit_transform(X)

# --- Ridge trace: coefficients versus lambda ---
# scikit-learn's Ridge minimizes ||y - Xb||^2 + alpha ||b||^2, so alpha is our lambda
lambdas = np.logspace(-3, 4, 200)
coefs = np.array([Ridge(alpha=a).fit(X_scaled, y).coef_ for a in lambdas])

plt.figure(figsize=(10, 5))
plt.semilogx(lambdas, coefs, linewidth=0.8)
plt.axhline(0, color="black", linewidth=0.5)
plt.xlabel("lambda")
plt.ylabel("Coefficient")
plt.title("Ridge trace")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# --- Choose lambda by 10-fold cross-validation ---
# The scaler sits inside the pipeline, so each fold is standardized with its own training data
model = make_pipeline(StandardScaler(),
                      RidgeCV(alphas=np.logspace(-3, 4, 100), cv=10,
                              scoring="neg_mean_squared_error"))
model.fit(X, y)
best_lambda = model[-1].alpha_
print(f"Optimal lambda: {best_lambda:.4f}")

# --- Compare prediction errors on the same folds ---
folds = KFold(n_splits=10, shuffle=True, random_state=0)
for name, reg in [("Least squares", LinearRegression()),
                  ("Ridge", Ridge(alpha=best_lambda))]:
    mse = -cross_val_score(make_pipeline(StandardScaler(), reg), X, y, cv=folds,
                           scoring="neg_mean_squared_error").mean()
    print(f"{name:13s} RMSECV = {np.sqrt(mse):.3f}")

# --- Compare least squares and ridge coefficients ---
ols = LinearRegression().fit(X_scaled, y)
ridge = Ridge(alpha=best_lambda).fit(X_scaled, y)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].bar(wavelengths, ols.coef_, width=6, color="tomato")
axes[0].set_title("Least squares coefficients")
axes[1].bar(wavelengths, ridge.coef_, width=6, color="steelblue")
axes[1].set_title(f"Ridge coefficients (lambda = {best_lambda:.3f})")
for ax in axes:
    ax.set_xlabel("Wavelength (nm)")
    ax.grid(True, alpha=0.3)
axes[0].set_ylabel("Coefficient")
plt.tight_layout()
plt.show()

% Requires the Statistics and Machine Learning Toolbox (cvpartition)

% --- Simulate spectra: analyte band + overlapping interferent + noise ---
rng(42);
n_samples = 40;
n_wavelengths = 30;
wavelengths = linspace(450, 700, n_wavelengths);
analyte = exp(-((wavelengths - 550).^2) / 3000);
interferent = exp(-((wavelengths - 600).^2) / 4000);
conc = 0.5 + (10.0 - 0.5) * rand(n_samples, 1);     % analyte concentration
c_int = 0.5 + (10.0 - 0.5) * rand(n_samples, 1);    % interferent, unrelated to y
X = conc * analyte + c_int * interferent + 0.05 * randn(n_samples, n_wavelengths);
y = conc + 0.2 * randn(n_samples, 1);               % reference method error

% Standardize predictors and centre the response
Xs = (X - mean(X)) ./ std(X);
yc = y - mean(y);

% --- Ridge trace ---
lambdas = logspace(-3, 4, 200);
coefs = zeros(numel(lambdas), n_wavelengths);
for i = 1:numel(lambdas)
    coefs(i, :) = ((Xs' * Xs + lambdas(i) * eye(n_wavelengths)) \ (Xs' * yc))';
end
figure;
semilogx(lambdas, coefs, 'LineWidth', 0.8);
xlabel('\lambda'); ylabel('Coefficient'); title('Ridge trace');
yline(0, 'k'); grid on;

% --- 10-fold cross-validation (standardizing inside each fold) ---
cvp = cvpartition(n_samples, 'KFold', 10);
cv_mse = zeros(numel(lambdas), 1);
for i = 1:numel(lambdas)
    sq_err = 0;
    for f = 1:cvp.NumTestSets
        tr = training(cvp, f);  te = test(cvp, f);
        mu = mean(X(tr, :));  sd = std(X(tr, :));
        Xtr = (X(tr, :) - mu) ./ sd;  Xte = (X(te, :) - mu) ./ sd;
        ybar = mean(y(tr));
        b = (Xtr' * Xtr + lambdas(i) * eye(n_wavelengths)) \ (Xtr' * (y(tr) - ybar));
        sq_err = sq_err + sum((y(te) - ybar - Xte * b).^2);
    end
    cv_mse(i) = sq_err / n_samples;
end
[~, best] = min(cv_mse);
best_lambda = lambdas(best);
fprintf('Optimal lambda: %.4f (RMSECV = %.3f)\n', best_lambda, sqrt(cv_mse(best)));
fprintf('Near least squares (lambda = %.3f): RMSECV = %.3f\n', lambdas(1), sqrt(cv_mse(1)));

% --- Compare least squares and ridge coefficients ---
beta_ols = Xs \ yc;
beta_ridge = (Xs' * Xs + best_lambda * eye(n_wavelengths)) \ (Xs' * yc);
figure;
subplot(1, 2, 1); bar(wavelengths, beta_ols, 'FaceColor', [0.9 0.4 0.3]);
title('Least squares coefficients'); xlabel('Wavelength (nm)'); grid on;
subplot(1, 2, 2); bar(wavelengths, beta_ridge, 'FaceColor', [0.3 0.5 0.8]);
title(sprintf('Ridge coefficients (\\lambda = %.3f)', best_lambda));
xlabel('Wavelength (nm)'); grid on;

library(glmnet)

# --- Simulate spectra: analyte band + overlapping interferent + noise ---
set.seed(42)
n_samples <- 40
n_wavelengths <- 30
wavelengths <- seq(450, 700, length.out = n_wavelengths)
analyte <- exp(-((wavelengths - 550)^2) / 3000)
interferent <- exp(-((wavelengths - 600)^2) / 4000)
conc <- runif(n_samples, 0.5, 10.0)       # analyte concentration
c_int <- runif(n_samples, 0.5, 10.0)      # interferent, unrelated to y
X <- outer(conc, analyte) + outer(c_int, interferent) +
  matrix(rnorm(n_samples * n_wavelengths, 0, 0.05), n_samples, n_wavelengths)
y <- conc + rnorm(n_samples, 0, 0.2)      # reference method error

# --- Ridge trace (alpha = 0 selects the pure L2 penalty) ---
# Note: glmnet scales the loss by 1/(2n) and standardizes X internally,
# so its lambda values are not on the same scale as the formula in the text.
ridge_fit <- glmnet(X, y, alpha = 0)
plot(ridge_fit, xvar = "lambda", main = "Ridge trace")

# --- Cross-validation to choose lambda ---
cv_ridge <- cv.glmnet(X, y, alpha = 0, nfolds = 10)
plot(cv_ridge)
cat("lambda.min:", cv_ridge$lambda.min, "\n")
cat("lambda.1se:", cv_ridge$lambda.1se, "\n")

# --- Compare least squares and ridge coefficients ---
beta_ols <- coef(lm(y ~ X))[-1]
beta_ridge <- as.numeric(coef(cv_ridge, s = "lambda.min"))[-1]

par(mfrow = c(1, 2))
barplot(beta_ols, names.arg = round(wavelengths), col = "tomato",
        main = "Least squares coefficients", xlab = "Wavelength (nm)")
barplot(beta_ridge, names.arg = round(wavelengths), col = "steelblue",
        main = "Ridge coefficients", xlab = "Wavelength (nm)")

The simulated spectra contain an analyte band that overlaps an interferent band, plus some error in the reference values, so least squares has something to get wrong. In the Python run, cross-validation picks a small but non-zero $\lambda$ (about 0.01), and the cross-validated error drops from 0.40 for least squares to 0.23 for ridge. The least squares coefficients swing up to about ±10 with alternating signs, while the ridge coefficients stay below 2 in absolute value. The MATLAB and R versions draw different random numbers, so their optimal $\lambda$ differs (about 1 in MATLAB, with the error falling from 0.35 to 0.23), but the picture is the same. Remember that the ridge error is slightly optimistic, because $\lambda$ was tuned on the same data.

Judging ridge by the calibration fit. Ridge regression deliberately fits the calibration data a little worse. Its calibration $R^2$ will always be lower than that of least squares, which tells you nothing. Compare models on validation data or by cross-validation.

Tuning and reporting on the same folds. If you choose $\lambda$ by cross-validation and then report that same cross-validated error as your expected performance, the number is optimistic, because $\lambda$ was picked to make it small. An independent test set, or a nested cross-validation, gives an honest estimate.

Over-interpreting individual coefficients. Ridge coefficients are shrunk, and their sizes depend on $\lambda$. Use them to understand the model, not as measurements of the "true" effect of each wavelength.

Expecting variable selection. Ridge shrinks every coefficient but sets none exactly to zero. If you need a short list of wavelengths, look at the lasso.

Ridge regression helps when:

• the predictors are strongly correlated, as in spectra;
• there are more predictors than samples ($p>n$): ${\mathbf{X}}^T\mathbf{X}+\lambda \mathbf{I}$ is always invertible for $\lambda >0$, so the solution is unique;
• the predictors are noisy and least squares would amplify that noise.

It does little or nothing when:

• the predictors are nearly orthogonal, for example from a designed experiment, where least squares is already stable;
• you need a sparse model with only a few variables;
• the model is misspecified: a penalty cannot fix a linear model fitted to a strongly nonlinear relationship.

In spectroscopic calibration, ridge regression competes directly with PCR and PLS. In an influential comparison, Frank and Friedman found that the three methods often perform similarly, with ridge frequently slightly ahead in their simulations, though by a small margin and with PLS and PCR close behind [8]. In chemometric practice PLS is more common, mostly because its latent variables are easier to inspect and interpret.

Ridge regression is least squares with a gentle preference for small coefficients. Mathematically it adds $\lambda$ to the diagonal of ${\mathbf{X}}^T\mathbf{X}$; geometrically it confines the coefficients to a sphere; in terms of the SVD it smoothly switches off the directions the data cannot determine; and from a Bayesian point of view it is a Gaussian prior centred at zero. All four descriptions say the same thing: trade a little bias for a large reduction in variance.

The same idea reappears in other forms:

• Lasso regression replaces the squared penalty by absolute values, which sets some coefficients exactly to zero.
• Elastic net mixes both penalties to combine the stability of ridge with the sparsity of the lasso [9].
• Principal component regression applies a hard version of the same shrinkage, keeping some directions completely and discarding the rest.
• Partial least squares also shrinks most of the poorly determined directions (though, unlike ridge and PCR, some of its factors can exceed 1), and it chooses them with the response in mind [8].

For a gentle textbook treatment, see James et al. [10]; for the full statistical picture, Hastie, Tibshirani and Friedman [7].