1.2 Model validation and ridge regression

Motivation

Least squares answers a question about a fixed feature matrix: which coefficients minimize the residual sum of squares? Choosing the features is a further question. A polynomial of sufficiently high degree can interpolate observations at distinct inputs, yet interpolation alone does not determine how accurately it predicts another response. Training residuals describe agreement with observations that helped choose the coefficients.

Model validation examines predictions at observations withheld from fitting. It makes the distinction between estimating coefficients and selecting a model operational: the training set determines a fitted function, while a validation set compares candidate feature choices or penalty strengths. After selection, the chosen procedure can be fitted again using all development data. A separate final test set evaluates that selected procedure without participating in the choice.

Ridge regression introduces another choice. It balances squared residuals against the squared size of the coefficients. Understanding this balance requires stating the loss scale, which coefficients are penalized, and how the features are measured. Cross-validation then compares these choices through prediction error rather than through the penalized training objective itself.

There is a simple reason why training error cannot settle the model choice. Suppose a new feature matrix contains every column of an old one together with additional columns. Every old coefficient vector extends to a new vector by assigning zero to the added coefficients, producing exactly the old predictions. Minimizing over the larger collection can therefore only decrease or preserve training SSE. The improvement could reflect a useful pattern, or it could absorb fluctuations specific to these observations. The training optimization has no withheld responses with which to distinguish those possibilities. Validation supplies a separate prediction comparison after the optimization has finished.

Definitions

Let $T$ and $V$ be disjoint, nonempty sets of training and validation indices. A model fitted using $T$ is denoted $\widehat f_T$. For squared error, define

$$L_T=\frac{1}{|T|}\sum_{i\in T}(y_i-\widehat f_T(x_i))^2, \qquad L_V=\frac{1}{|V|}\sum_{i\in V}(y_i-\widehat f_T(x_i))^2.$$

The same fitted function occurs in both expressions. Computing validation loss does not refit its coefficients using the validation responses. A candidate model includes the feature construction, preprocessing, and fitting rule, so these parts must also be learned within each training split.

Measuring prediction error

For an evaluation set of size $m$, write $e_i=y_i-\widehat y_i$. Three common metrics are

$$\operatorname{MSE}=\frac1m\sum_i e_i^2, \qquad \operatorname{MAE}=\frac1m\sum_i|e_i|, \qquad \operatorname{RMSE}=\sqrt{\frac1m\sum_i e_i^2}.$$

MSE has squared response units; MAE and RMSE have response units. Squaring gives large residuals more influence, whereas absolute loss grows linearly with their magnitude. RMSE and MSE rank predictions identically on the same evaluation set because the square root is increasing. MAE may give a different ranking. Choose the metric to express the prediction task before comparing candidates.

With the evaluation-set mean $\overline y$, define

$$\operatorname{SST}=\sum_i(y_i-\overline y)^2, \qquad R^2=1-\frac{\sum_i e_i^2}{\operatorname{SST}}, \qquad \operatorname{SST}>0.$$

Here the baseline always predicts that evaluation-set mean. Negative $R^2$ means the candidate has larger squared error than this baseline. On training data, ordinary least squares in a model containing a free constant has $0\le R^2\le1$: the constant predictor is among its competitors. That comparison does not apply to held-out predictions. When all evaluation responses coincide, SST is zero and this definition has no value; MSE remains available. In particular, a one-observation validation fold cannot supply this $R^2$.

To justify the training bound, evaluate the least-squares objective at the response mean. Its SSE equals SST, so the minimizing SSE lies between zero and SST. Division by positive SST gives the bound. This argument uses the same observations for fitting and computing the two sums. A validation response mean is a comparison baseline, not a coefficient estimated from the training subset. Also, averaging fold RMSEs and taking the square root of pooled MSE generally differ: the square root is applied at different stages of the calculation.

Median and percentage errors

Two further summaries use the middle of the error distribution:

$$\operatorname{MdSE}=\operatorname{median}_i(e_i^2), \qquad \operatorname{MdAE}=\operatorname{median}_i|e_i|.$$

For an even sample size, average the two middle ordered values. MdSE has squared response units; MdAE has response units. Means reflect the average error across all observations, while medians describe its middle magnitude. The two summaries can therefore distinguish different features of a residual distribution.

For $y=(1,2,3,7)^T$ and predictions $(1,2,0,0)^T$, residuals are $(0,0,3,7)^T$. MSE is $14.5$ and MdSE is $4.5$; MAE is $2.5$ and MdAE is $1.5$. The mean absolute percentage error is

$$\operatorname{MAPE}=\frac{100}{m}\sum_{i=1}^m\left|\frac{e_i}{y_i}\right|, \qquad y_i\ne0\quad(1\le i\le m).$$

MAPE expresses relative error as a percentage. All four observed responses above are nonzero, and their absolute relative errors are $0,0,1,1$, giving $50\%$. This definition requires every observed response to be nonzero.

Learning a feature scale

For a nonconstant training feature, its training mean and population standard deviation give

$$\mu_T=\frac1{|T|}\sum_{i\in T}x_i, \qquad s_T=\sqrt{\frac1{|T|}\sum_{i\in T}(x_i-\mu_T)^2}, \qquad z=\frac{x-\mu_T}{s_T}.$$

Use the same $\mu_T$ and $s_T$ when transforming validation inputs. Training min-max scaling similarly uses training extrema and requires a positive training range. A constant feature can be removed before standardization; the intercept already represents a constant contribution. A future transformed value can lie outside the training interval. Its location records a new input, rather than an error in the transformation.

For $A\in\mathbb R^{n\times p}$ and $y\in\mathbb R^n$, the fully penalized ridge objective in this section is

$$J_\lambda(\beta)=\|y-A\beta\|_2^2+\lambda\|\beta\|_2^2, \qquad \lambda>0.$$

This uses the unnormalized sum of squared errors. An objective $\|y-A\beta\|_2^2/n+\rho\|\beta\|_2^2$ has the same minimizer when $\lambda=n\rho$. If the penalty is written $\alpha^2\|\beta\|_2^2$, then $\lambda=\alpha^2$. These parameters are numerical descriptions of the chosen objective, so their normalization belongs with the reported value.

Theorem/proposition

Proposition

A unique ridge fit at every positive penalty

For every real $A$ and $y$ and every $\lambda>0$, $A^TA+\lambda I$ is positive definite. The unique minimizer of $J_\lambda$ is characterized by

$$(A^TA+\lambda I)\widehat\beta_\lambda=A^Ty.$$

This conclusion includes rank-deficient design matrices. Positive definiteness provides both existence and uniqueness of the coefficient vector.

Proposition

The training tradeoff as the penalty increases

Fix $A$ and $y$, and let $0<\lambda_1<\lambda_2$. Define

$$S_j=\|y-A\widehat\beta_{\lambda_j}\|_2^2, \qquad P_j=\|\widehat\beta_{\lambda_j}\|_2^2.$$

Then $P_2\le P_1$ and $S_1\le S_2$. A larger penalty reduces or preserves the coefficient norm and increases or preserves the training residual sum of squares. The inequalities make no ordering claim about validation error.

Proof sketch or proof idea

For any nonzero vector $h$,

$$h^T(A^TA+\lambda I)h=\|Ah\|_2^2+\lambda\|h\|_2^2>0.$$

Differentiation of $J_\lambda$ gives $2(A^TA+\lambda I)\beta-2A^Ty$. Its unique zero solves the displayed linear system. More directly, expansion around this solution yields

$$J_\lambda(\widehat\beta_\lambda+h) =J_\lambda(\widehat\beta_\lambda) +h^T(A^TA+\lambda I)h.$$

Every nonzero change therefore increases the objective. This argument proves a global minimum, including when ordinary least squares has several coefficient solutions.

The same ridge problem is an ordinary least-squares problem on augmented observations:

$$\widetilde A=\begin{bmatrix}A\\\sqrt\lambda I\end{bmatrix}, \qquad \widetilde y=\begin{bmatrix}y\\0\end{bmatrix}, \qquad J_\lambda(\beta)=\|\widetilde y-\widetilde A\beta\|_2^2.$$

The additional rows contribute exactly the coefficient penalty. More generally, replacing $\sqrt\lambda I$ by a matrix $\Gamma$ gives a Tikhonov penalty $\|\Gamma\beta\|_2^2$. Its minimizer is unique precisely when the only vector annihilated by both $A$ and $\Gamma$ is zero. Indeed, the corresponding quadratic form is $\|Ah\|_2^2+\|\Gamma h\|_2^2$. This criterion explains why leaving an intercept unpenalized can still give a unique fit: the data matrix observes the constant direction even though the penalty does not constrain it.

For the second proposition, optimality at the two penalty values gives

$$S_1+\lambda_1P_1\le S_2+\lambda_1P_2, \qquad S_2+\lambda_2P_2\le S_1+\lambda_2P_1.$$

Adding yields $(\lambda_2-\lambda_1)(P_2-P_1)\le0$, so $P_2\le P_1$. The first inequality then implies $S_1-S_2\le\lambda_1(P_2-P_1)\le0$.

The spectrum explains which directions shrink. If $A^TA v_j=d_jv_j$ for orthonormal eigenvectors and $d_j>0$, then

$$v_j^T\widehat\beta_\lambda =\frac{v_j^TA^Ty}{d_j+\lambda} =\frac{d_j}{d_j+\lambda}\,v_j^T\widehat\beta_0^{+}.$$

Here $\widehat\beta_0^{+}$ is the minimum-norm ordinary least-squares solution. A zero-eigenvalue direction has $Av_j=0$, hence $v_j^TA^Ty=0$ and ridge coefficient zero in that direction. Small positive $d_j$ receive stronger relative shrinkage. This is a statement about orthogonal spectral coordinates; an individual coefficient in the original feature coordinates need not move monotonically toward zero.

Ridge with a coefficient-norm budget

The same fitting preference can be expressed through a constraint. Write $f(\beta)=\|y-A\beta\|_2^2$ and consider

$$\min_{\beta:\,\|\beta\|_2^2\le\gamma}f(\beta), \qquad \gamma\ge0.$$

This is a ball of radius $\sqrt\gamma$. For a positive penalty, let $\gamma_\lambda=\|\widehat\beta_\lambda\|_2^2$. Any vector $b$ satisfying that budget obeys

$$f(\widehat\beta_\lambda)+\lambda\gamma_\lambda \le f(b)+\lambda\|b\|_2^2 \le f(b)+\lambda\gamma_\lambda.$$

Canceling the penalty bound shows that the ridge vector minimizes the constrained loss. This construction links the selected penalty to its budget through the fitted coefficient norm.

Conversely, take $\gamma>0$. The zero vector is strictly feasible, so the convex differentiable problem satisfies Slater's condition. Its Karush–Kuhn–Tucker (KKT) conditions are necessary and sufficient:

$$\begin{aligned} A^T(A\beta_*-y)+\lambda\beta_*&=0,\\ \|\beta_*\|_2^2&\le\gamma,\\ \lambda&\ge0,\\ \lambda(\|\beta_*\|_2^2-\gamma)&=0. \end{aligned}$$

Here the Lagrangian is $f(\beta)+\lambda(\|\beta\|_2^2-\gamma)$. Dividing its stationarity equation by two gives the displayed first line, with the same SSE penalty convention used above. Convexity makes a stationary point minimize the Lagrangian. Feasibility and complementary slackness then certify the constrained minimum. A positive multiplier gives the unique ridge vector; at a zero multiplier the stationarity equation is ordinary least squares, whose coefficients need not be unique.

A slack norm constraint forces $\lambda=0$. A constraint holding with equality can also have zero multiplier when an ordinary least-squares minimizer lies on the boundary. Positive multipliers therefore do not describe every possible budget. At $\gamma=0$, strict feasibility fails and the only feasible vector is zero; a finite ridge penalty need not reach it. If the intercept is unpenalized, constrain only the slope vector and leave the intercept free.

For a scalar example, take $A=[1]$ and $y=2$. Then

$$\widehat\beta_\lambda=\frac2{1+\lambda},\qquad \gamma_\lambda=\frac4{(1+\lambda)^2},\qquad \beta_*^{(\gamma)}=\min(2,\sqrt\gamma).$$

The constrained solution projects $2$ onto $[-\sqrt\gamma,\sqrt\gamma]$. Budget $\gamma=1$ gives $\beta_*=1$ and $\lambda=1$. Budget $\gamma=4$ gives $\beta_*=2$ on the boundary with $\lambda=0$; budget $\gamma=9$ gives the same coefficient strictly inside the ball, again with zero multiplier. Budget $\gamma=0$ forces $\beta_*=0$, while $2/(1+\lambda)>0$ for every finite $\lambda\ge0$. In this example, the zero budget is attained only in the limit $\lambda\to\infty$.

Worked examples

Worked example

Example 1: transform a withheld input with the training scale

Suppose one training feature takes values $0$ and $2$, and a validation input is $4$. The training mean is $1$ and population standard deviation is $1$, so the training standardized values are $-1$ and $1$; the validation standardized value is $3$.

Training min-max scaling instead gives $0$, $1$, and validation value $2$. Both procedures reuse their training parameters. Recomputing the mean or range after inserting the validation input changes the transformation used to define the fitted model. It lets the held-out observations influence a stage of fitting before their prediction errors are measured.

If another cross-validation fold has different training observations, recompute its scale from those observations. The preprocessing belongs to the fitting procedure, rather than to a single preliminary pass over the complete data.

Worked example

Example 2: a constant model under leave-one-out validation

Fit a constant to responses $0$, $2$, and $4$. Using all three gives $\widehat c=2$, residual sum of squares $8$, and training MSE $8/3$.

For leave-one-out cross-validation, refit the constant three times. Holding out $0$ leaves training mean $3$; holding out $2$ leaves mean $2$; holding out $4$ leaves mean $1$. Thus

$$\widehat y^{(-1)}=3,\qquad \widehat y^{(-2)}=2,\qquad \widehat y^{(-3)}=1, \qquad \operatorname{CV}_{\mathrm{LOO}}=\frac{9+0+9}{3}=6.$$

Each held-out prediction was obtained without that response. Using the full-data mean for all three predictions would reproduce the training MSE, not leave-one-out error. The calculation compares a full-data fit with three fits using two observations; their difference is visible even for a model with just one coefficient.

Worked example

Example 3: ridge separates residual error from coefficient size

Take $A=\operatorname{diag}(1,2)$, $y=(2,2)^T$, and $\lambda=1$, with both coefficients penalized. Then

$$A^TA+I=\begin{bmatrix}2&0\\0&5\end{bmatrix}, \qquad A^Ty=\begin{bmatrix}2\\4\end{bmatrix}, \qquad \widehat\beta_1=\begin{bmatrix}1\\4/5\end{bmatrix}.$$

The predictions are $(1,8/5)^T$ and residuals $(1,2/5)^T$. Therefore

$$\|y-A\widehat\beta_1\|_2^2=\frac{29}{25}, \qquad \|\widehat\beta_1\|_2^2=\frac{41}{25}, \qquad J_1(\widehat\beta_1)=\frac{14}{5}.$$

Ordinary least squares gives $(2,1)^T$, with zero residual error and squared coefficient norm $5$. Its value in the same $J_1$ objective is $5$, exceeding $14/5$. Ridge accepts residual error to reduce the combined objective. The two spectral shrinkage factors are $1/2$ and $4/5$, matching the two coefficients above.

These two computations can be reproduced directly:

import numpy as np

y = np.array([0.0, 2.0, 4.0])
held_predictions = np.array([
    np.delete(y, i).mean() for i in range(len(y))
])
train_mse = np.mean((y - y.mean()) ** 2)
loo_mse = np.mean((y - held_predictions) ** 2)

A = np.diag([1.0, 2.0])
target = np.array([2.0, 2.0])
beta = np.linalg.solve(A.T @ A + np.eye(2), A.T @ target)
residual = target - A @ beta
objective = residual @ residual + beta @ beta
print(train_mse, loo_mse, beta, objective)

The outputs are approximately $2.666667$, $6$, $(1,0.8)$, and $2.8$. Solving the positive definite system avoids explicitly computing its inverse.

Choosing a model with cross-validation

For $n\ge2$, leave-one-out validation trains on all but one observation and averages the $n$ squared held-out errors. Leave-$p$-out, with $1\le p<n$, considers every validation subset of size $p$, requiring $\binom np$ fits per candidate. In $k$-fold validation, with $2\le k\le n$, partition observations into $k$ disjoint nonempty blocks; each block is withheld once. When $k=n$, this is leave-one-out validation.

For unequal fold sizes, pool squared errors and divide by the total number of held-out observations to weight observations equally. Averaging fold MSEs equally instead weights folds equally. Repeated random holdout, or Monte Carlo cross-validation, samples several training-validation partitions. Within each split the sets are disjoint, while an observation may be validated repeatedly across splits or never be selected. Use the same splits for all candidate penalties to make their comparisons refer to the same prediction tasks.

The code below uses three fixed folds, one nonconstant feature, and an unpenalized intercept. Its penalty matrix $P=\operatorname{diag}(0,1)$ gives objective $\|y-A\beta\|_2^2+\lambda\beta_1^2$. Every training subset contains varying feature values, so each computed standard deviation is positive.

import numpy as np

X = np.arange(6.0)[:, None]
y = np.array([1.0, 2.0, 2.0, 4.0, 4.0, 5.0])
folds = [np.array([0, 1]), np.array([2, 3]), np.array([4, 5])]
indices = np.arange(len(y))
penalty = np.diag([0.0, 1.0])

for lam in (0.0, 1.0, 4.0):
    squared_error = 0.0
    count = 0
    for valid in folds:
        train = np.setdiff1d(indices, valid)
        mu = X[train].mean(axis=0)
        scale = X[train].std(axis=0, ddof=0)
        Z_train = (X[train] - mu) / scale
        Z_valid = (X[valid] - mu) / scale
        A = np.column_stack((np.ones(len(train)), Z_train))
        B = np.column_stack((np.ones(len(valid)), Z_valid))
        beta = np.linalg.solve(A.T @ A + lam * penalty, A.T @ y[train])
        error = y[valid] - B @ beta
        squared_error += error @ error
        count += len(valid)
    print(lam, squared_error / count)

For $\lambda=0,1,4$, the pooled MSEs are approximately $0.263829$, $0.179066$, and $0.808482$, respectively. Thus $\lambda=1$ wins among these three candidates. Select the candidate with the smallest pooled validation MSE, then learn a fresh scale and fit using all development observations. The intercept is free to represent the response level. With positive $\lambda$, the only unpenalized direction here is the intercept, whose nonzero constant column ensures uniqueness. The earlier $\lambda I$ theorem penalizes every coefficient; changing the penalty matrix changes that stated problem.

If the candidate is instead defined by a mean-loss penalty strength $\rho$, use $\lambda_T=|T|\rho$ inside each training split. Training subsets of different sizes then implement the same mean-loss convention. Keeping the numerical SSE parameter fixed implements a different comparison. Once candidates have been compared repeatedly against validation observations, those observations have helped select the procedure. Refitting the selected procedure on the full development set uses that information constructively; a final test set answers the subsequent evaluation question. A fixed split makes a numerical comparison reproducible, while several chosen splits reveal how strongly the comparison depends on which observations were withheld.

Feature scaling also changes the penalty. If $z_j=(x_j-\mu_j)/s_j$ and the standardized coefficient is $\theta_j$, the raw slope is $\beta_j=\theta_j/s_j$. Penalizing $\sum_j\theta_j^2$ therefore penalizes $\sum_j s_j^2\beta_j^2$ in raw coordinates. Standardization gives a chosen, consistent unit of comparison; it is part of the model definition. Sensitivity can then be explored by perturbing a feature or penalty and recomputing predictions while keeping other choices fixed.

Sensitivity to input changes

Fixed-model sensitivity changes inputs while keeping fitted parameters fixed. Refitting after changing training observations also changes the model itself; its resulting predictions measure a different form of stability.

One-at-a-time (OAT) analysis changes one input while fixing the others. For the fixed interaction model $f(x_1,x_2)=x_1x_2$ at $(0,0)$,

$$f(1/5,0)=f(0,1/5)=0, \qquad f(1/5,1/5)=1/25.$$

Neither individual perturbation changes the output, but the joint perturbation does. Unchanged OAT outputs therefore do not establish joint insensitivity. Multivariable analysis examines interactions; scenario analysis specifies joint changes, while Monte Carlo sensitivity samples joint inputs from a specified distribution.

Common mistakes

Common mistake

Compare prediction error after the complete fitting procedure

Training SSE generally favors more flexible models and weaker penalties. Validation compares predictions after fitting and preprocessing on each training subset. The ridge objective includes a coefficient penalty, so its numerical value is not a validation MSE. Increasing $\lambda$ guarantees the training tradeoff proved above, while the location of the best validation error depends on the data.

Summary

Define the evaluation metric, split observations, learn preprocessing on the training subset, and fit each candidate there. Pool held-out errors to compare candidates. For ridge, record the loss normalization and intercept convention with the penalty value; positive full-coefficient regularization gives a unique solution and shrinks spectral directions according to their eigenvalues.

Exercises

Checkpoint

Exercise 1: training feature values are $0$ and $2$. Standardize a validation input $4$ using the population standard deviation. What are the three transformed values?

Compute the mean and standard deviation using the two training values.

Checkpoint

Exercise 2: evaluation responses are $0$ and $2$, with predictions $3$ and $3$. Find MSE, MAE, RMSE, and $R^2$.

Use the evaluation-set mean in SST and retain the sign of $R^2$.

Checkpoint

Exercise 3: let $A=\begin{bmatrix}1&1\\1&1\end{bmatrix}$, $y=(2,2)^T$, and $\lambda=2$. Penalize both coefficients. Find the unique ridge coefficient vector.

Form $A^TA+2I$ even though the two columns of $A$ coincide.

Checkpoint

Exercise 4: a training set has $n=5$ and objective $\|y-A\beta\|_2^2/5+0.2\|\beta\|_2^2$. What is the equivalent $\lambda$ in the unnormalized SSE convention?

Multiply the entire objective by $5$.

Checkpoint

Exercise 5: For $A=[1]$, $y=2$, and budget $\gamma=1/4$, find the ball radius, constrained minimizer, and corresponding SSE ridge parameter.

Project the unconstrained coefficient onto the ball, then solve for the penalty.

Solutions

Solution · Exercise 1 solution

The training mean is $1$ and population standard deviation is $1$. Thus $(0-1)/1=-1$, $(2-1)/1=1$, and $(4-1)/1=3$. The validation input does not change the scale parameters.

Solution · Exercise 2 solution

The residuals are $-3$ and $-1$, so SSE is $10$, MSE is $5$, MAE is $2$, and RMSE is $\sqrt5$. The response mean is $1$ and SST is $2$, giving $R^2=1-10/2=-4$. The negative value expresses performance below the constant-mean baseline on these observations.

Solution · Exercise 3 solution

The system is

$$\begin{bmatrix}4&2\\2&4\end{bmatrix}\beta =\begin{bmatrix}4\\4\end{bmatrix}.$$

Subtracting the two equations gives $\beta_1=\beta_2$, and either equation then gives $\beta_1=\beta_2=2/3$. The matrix has eigenvalues $6$ and $2$, both positive. The duplicate original columns do not prevent the ridge solution from being unique.

Solution · Exercise 4 solution

Multiplying by $5$ gives $\|y-A\beta\|_2^2+\|\beta\|_2^2$, so $\lambda=1$. Multiplication by a positive constant preserves minimizers; rescaling only the residual term would change their tradeoff.

Solution · Exercise 5 solution

The radius is $\sqrt\gamma=1/2$, so $\beta_*=1/2$. Solving $2/(1+\lambda)=1/2$ gives $\lambda=3$. Stationarity is $\beta_*-2+\lambda\beta_*=1/2-2+3/2=0$, and the positive multiplier matches the binding budget.

Related notes

Least squares and residuals develops the fitting equations used here. Lasso and matching pursuit studies penalties and algorithms that select sparse coefficients. SVD and principal components explains the spectral directions behind shrinkage. The course reading list collects further reading.