Holding out rows

Choosing Hyperparameters

Most methods in this course come with a knob that the data do not set by themselves: the regularization parameter $\lambda$ in ridge regression, the number of principal components to keep, the degree of a polynomial model, the number of clusters. These notes are about how to set such a knob. The main tool is cross-validation, and we watch it choose the regularization parameter of ridge regression on a deblurring problem.

Parameters and hyperparameters

Definition 1 (Hyperparameter) A hyperparameter is a number that must be fixed before a method can be fit to data, as opposed to a parameter, which the fitting itself determines.

In ridge regression, \[ \widehat{\bftheta}_\lambda=\argmin_{\bftheta\in\Rbb^d}\Big\{\|\bfy-\bfX\bftheta\|_2^2+\lambda\|\bftheta\|_2^2\Big\}, \] the coefficients $\bftheta$ are parameters and $\lambda$ is a hyperparameter: once $\lambda$ is fixed, the ridge solution determines $\widehat{\bftheta}_\lambda$, but nothing in the objective tells us which $\lambda$ to use.

Lecture 10 tells us what the best $\lambda$ would be. Under the standard assumptions, the mean squared error is \[ \operatorname{MSE}(\widehat{\bftheta}_\lambda) =\sum_{i=1}^d\frac{\lambda^2\alpha_i^2}{(\sigma_i^2+\lambda)^2} +\sigma^2\sum_{i=1}^d\frac{\sigma_i^2}{(\sigma_i^2+\lambda)^2}, \qquad \alpha_i=\bfv_i^\top\bftheta, \] and we would choose $\lambda$ to make this as small as possible. But it depends on the true $\bftheta$ and the noise level $\sigma^2$, which we do not know. We need a yardstick that can be computed from the data alone.

Training error is the wrong yardstick

The first idea is to choose the $\lambda$ whose fit matches the data best. This never works. Throughout, $\bfX\in\Rbb^{n\times d}$ has rank $r\geq1$ and compact SVD $\bfX=\sum_{i=1}^r\sigma_i\bfu_i\bfv_i^\top$.

Proposition 1 (Training error always prefers $\lambda=0$) The training residual of ridge regression satisfies \[ \|\bfy-\bfX\widehat{\bftheta}_\lambda\|_2^2 =\sum_{i=1}^r\Big(\frac{\lambda}{\sigma_i^2+\lambda}\Big)^2(\bfu_i^\top\bfy)^2+\|\bfy-\bfP\bfy\|_2^2, \] where $\bfP=\bfU_r\bfU_r^\top$ is the orthogonal projection onto $\Ccal(\bfX)$. In particular, it is an increasing function of $\lambda$.

Proof

Proof. By ridge regression filters the singular values and $\bfX\bfv_i=\sigma_i\bfu_i$, the fitted values are $\bfX\widehat{\bftheta}_\lambda=\sum_{i=1}^r\frac{\sigma_i^2}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\bfu_i$. Since $\bfy=\sum_{i=1}^r(\bfu_i^\top\bfy)\bfu_i+(\bfy-\bfP\bfy)$, \[ \bfy-\bfX\widehat{\bftheta}_\lambda =\sum_{i=1}^r\frac{\lambda}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\bfu_i+(\bfy-\bfP\bfy). \] The vectors $\bfu_1,\ldots,\bfu_r$ are orthonormal and $\bfy-\bfP\bfy\perp\Ccal(\bfX)$, so the Pythagorean theorem gives the formula. Each $\lambda/(\sigma_i^2+\lambda)$ increases with $\lambda$.

Of course: least squares is defined as the minimizer of the training error, so anything else does worse on the training data. The same trap appears for every hyperparameter that controls how flexible a model is. More principal components, a higher polynomial degree, or more clusters always fit the training data at least as well, so the training error always asks for the most flexible model, including its ability to fit the noise. To judge a hyperparameter honestly, we must test the fit on data that was not used to compute it.

Validation and cross-validation

The simplest option is to hold out a validation set: fit on most of the data, measure the error on the rest, and pick the $\lambda$ with the smallest validation error. This wastes data, and the answer depends on which points happened to be held out. Cross-validation fixes both problems by rotating the role of validation set through the data.

Algorithm 1 ($K$-fold cross-validation) Input: data $(\bfX,\bfy)$, a method that fits an estimate $\widehat{\bftheta}_\lambda$ to data for each value of a hyperparameter $\lambda$, a grid of candidate values $\lambda_1,\ldots,\lambda_m$, and a number of folds $K$ (typically $5$ or $10$).

Randomly shuffle the indices $\{1,\ldots,n\}$ and split them into $K$ disjoint folds $\Fcal_1,\ldots,\Fcal_K$ of nearly equal size.
For each candidate $\lambda$ and each fold $k=1,\ldots,K$:
1. Fit $\widehat{\bftheta}_\lambda^{(-k)}$ using only the observations not in $\Fcal_k$.
2. Compute its error on the held-out fold, $\displaystyle\mathrm{Err}_k(\lambda)=\frac{1}{|\Fcal_k|}\sum_{i\in\Fcal_k}\big(y_i-\bfx_i^\top\widehat{\bftheta}_\lambda^{(-k)}\big)^2$.
Average over the folds: $\displaystyle\mathrm{CV}(\lambda)=\frac1K\sum_{k=1}^K\mathrm{Err}_k(\lambda)$.
Choose $\lambda^\ast=\argmin_{\lambda_j}\mathrm{CV}(\lambda_j)$, and refit on all the data with $\lambda^\ast$.

Every observation is used for validation exactly once and for fitting $K-1$ times. The choice $K=n$ is called leave-one-out cross-validation. Nothing in the algorithm needs a formula for the estimate, only a way to fit it to part of the data and to predict the rest, so it works for any method. It can even choose between methods, if we treat the method itself as one more hyperparameter. For ridge regression the grid of $\lambda$ values is usually logarithmically spaced, say $10^{-3},\ldots,10^{3}$, since $\lambda$ only matters relative to the $\sigma_i^2$.

$\mathrm{CV}(\lambda)$ is itself an average of $K$ noisy numbers, so its minimum is not exact. The spread of the fold errors measures how uncertain it is, and suggests a more cautious choice (Hastie et al. 2009, sec. 7.10).

Remark (The one-standard-error rule)

Remark (The one-standard-error rule). Let $\mathrm{SE}(\lambda)$ be the sample standard deviation of $\mathrm{Err}_1(\lambda),\ldots,\mathrm{Err}_K(\lambda)$, divided by $\sqrt{K}$, and let $\lambda^\ast$ minimize $\mathrm{CV}$. Values of $\lambda$ whose $\mathrm{CV}(\lambda)$ is within $\mathrm{SE}(\lambda^\ast)$ of the minimum cannot really be told apart from $\lambda^\ast$. Among them, choose the one that regularizes most, which for ridge regression is the largest.

Keep the held-out data held out

Anything learned from the data counts as fitting, including the means and standard deviations used to standardize the features. Compute them from the training folds only, then apply them to the held-out fold. And if you want an honest estimate of the final model’s error, you need a separate test set that played no role in choosing $\lambda$ at all.

Cross-validation for ridge regression in NumPy

Using the same folds for every $\lambda$ makes the comparison fair. Solving with $\bfX^\top\bfX+\lambda\bfI_d$ directly is acceptable here, because $\lambda>0$ keeps its condition number under control.

import numpy as np

def ridge(X, y, lam):
    d = X.shape[1]
    return np.linalg.solve(X.T @ X + lam * np.eye(d), X.T @ y)

def cv_error(X, y, lam, K=5, seed=0):
    rng = np.random.default_rng(seed)
    folds = np.array_split(rng.permutation(len(y)), K)
    errors = []
    for k in range(K):
        test = folds[k]
        train = np.concatenate([folds[j] for j in range(K) if j != k])
        theta = ridge(X[train], y[train], lam)
        errors.append(np.mean((y[test] - X[test] @ theta) ** 2))
    return np.mean(errors)

lambdas = np.logspace(-3, 3, 50)
best_lam = min(lambdas, key=lambda lam: cv_error(X, y, lam))

In scikit-learn, RidgeCV(alphas=lambdas, cv=5) does the same search, and GridSearchCV does it for any model and any hyperparameters.

What cross-validation measures

Cross-validation scores a fit by how well it predicts new observations. That is not the same as how close $\widehat{\bftheta}$ is to $\bftheta$, and the difference is worth making precise.

Proposition 2 (Prediction error) Let $\widehat{\bftheta}$ be any estimate of $\bftheta$, and let $y=\bfx^\top\bftheta+\varepsilon$ be a new observation at a fixed $\bfx\in\Rbb^d$, whose noise $\varepsilon$ has mean $0$ and variance $\sigma^2$ and is independent of $\widehat{\bftheta}$. Then \[ \Ebb\big[(y-\bfx^\top\widehat{\bftheta})^2\big] =\sigma^2+\Ebb\big[\big(\bfx^\top(\widehat{\bftheta}-\bftheta)\big)^2\big]. \]

Proof

Proof. Write $y-\bfx^\top\widehat{\bftheta}=\varepsilon-\bfx^\top(\widehat{\bftheta}-\bftheta)$ and expand the square: \[ \Ebb\big[(y-\bfx^\top\widehat{\bftheta})^2\big] =\Ebb[\varepsilon^2]-2\,\Ebb\big[\varepsilon\,\bfx^\top(\widehat{\bftheta}-\bftheta)\big]+\Ebb\big[\big(\bfx^\top(\widehat{\bftheta}-\bftheta)\big)^2\big]. \] The first term is $\Var(\varepsilon)=\sigma^2$. By independence, the middle expectation is $\Ebb[\varepsilon]\,\Ebb[\bfx^\top(\widehat{\bftheta}-\bftheta)]=0$.

In $K$-fold cross-validation, the noise in the held-out fold is independent of the fit $\widehat{\bftheta}_\lambda^{(-k)}$, which never saw it. So $\mathrm{CV}(\lambda)$ estimates $\sigma^2$, which no choice of $\lambda$ can change, plus the average prediction error $(\bfx_i^\top(\widehat{\bftheta}-\bftheta))^2$ over the held-out rows. In particular, the cross-validation error cannot go much below $\sigma^2$.

The prediction error weighs the error in each direction by how well the data see that direction. Averaging over all the rows of $\bfX$ and writing $\widehat{\bftheta}-\bftheta$ in the basis $\bfv_1,\ldots,\bfv_d$ (assuming (A3), so $r=d$), \[ \frac1n\|\bfX(\widehat{\bftheta}-\bftheta)\|_2^2 =\frac1n\sum_{i=1}^d\sigma_i^2\big(\bfv_i^\top(\widehat{\bftheta}-\bftheta)\big)^2, \qquad\text{while}\qquad \|\widehat{\bftheta}-\bftheta\|_2^2 =\sum_{i=1}^d\big(\bfv_i^\top(\widehat{\bftheta}-\bftheta)\big)^2. \] An error along a poorly measured direction, with $\sigma_i$ small, is nearly invisible to cross-validation, even when it dominates $\|\widehat{\bftheta}-\bftheta\|_2^2$. When the goal is prediction, this is exactly right: an error the data cannot see will not hurt future predictions either. When the goal is $\bftheta$ itself, as in the example below, cross-validation tends to be too lenient toward those errors.

Example: deblurring a photograph

We return to the out-of-focus camera from the variance of least squares, this time pointed at a cat. A greyscale photograph $\bfTheta$ of Chelsea, the cat that ships with scikit-image, is blurred by the same symmetric matrix $\bfX$, whose row $j$ is a narrow Gaussian bump centred on pixel $j$, and the camera adds a little noise. We observe \[ \bfY=\bfX\bfTheta\bfX^\top+\bfE, \qquad E_{pq}\ \text{independent},\ \mathcal N(0,\sigma^2), \] where the noise level $\sigma$ is small compared to the range from black to white. The smallest singular values of the blur are many orders of magnitude below the largest, so least squares turns the noise into static, and we must regularize. The question is how much.

This example has one luxury that real problems do not: we have the sharp photograph. So for every $\lambda$ we can compute the true error, the mean squared difference per pixel between the estimate and the sharp photograph, and grade the choice that cross-validation makes without ever showing it the answer. The code is folded throughout; the figures are the point.

Set up the blurred photograph

import base64
import io

import numpy as np
import plotly.graph_objects as go
from matplotlib import colormaps
from PIL import Image
from plotly.subplots import make_subplots
from skimage import color, data
from threadpoolctl import threadpool_limits

threadpool_limits(1)  # this page does many small matrix products, which run fastest on one thread
rng = np.random.default_rng(123)

# The blur: row j of X is a Gaussian bump of width w centred on pixel j.
d = 100
w = 2.0
pixels = np.arange(d)
X = np.exp(-(pixels[:, None] - pixels[None, :])**2 / (2 * w**2)) / (np.sqrt(2 * np.pi) * w)
sigma = 0.01  # the sensor noise: 1% of the range from black to white

# scikit-image's photograph of Chelsea the cat in greyscale, cropped from 300 x 451 to a
# 300 x 300 square around her face, then averaged down to 100 x 100 to keep the page fast.
photo = color.rgb2gray(data.chelsea())[:, 90:390].reshape(d, 3, d, 3).mean(axis=(1, 3))

# Blur every column, then every row, then add sensor noise to every pixel.
blurred = X @ photo @ X.T + sigma * rng.standard_normal((d, d))

U, s, Vt = np.linalg.svd(X)
V = Vt.T
XV = X @ V  # equals U diag(s)

Figure style and image helpers

# The same muted palette as the examples, which reads correctly against
# both the light and the dark version of this site.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
BLUE, ORANGE = "#2a78d6", "#eb6834"
GAP = 6  # transparent pixels between images placed side by side


def style(fig, height=420, title=None, showlegend=True):
    """Apply the shared layout to a plot with axes."""
    fig.update_layout(
        template="simple_white",
        paper_bgcolor="rgba(0,0,0,0)",
        plot_bgcolor="rgba(0,0,0,0)",
        font=dict(color=INK, size=13),
        title=dict(text=title, x=0, xanchor="left", font=dict(size=15)),
        margin=dict(l=60, r=30, t=60 if title else 40, b=50),
        height=height,
        showlegend=showlegend,
        hovermode="closest",
        legend=dict(orientation="h", y=-0.18, x=0),
    )
    fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID, showgrid=True)
    fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID, showgrid=True)
    return fig


def png_uri(array):
    """Encode a greyscale array as a PNG data URI, clipping brightness to [0, 1].

    NaN pixels become transparent. A PNG is far smaller in the page than a heatmap.
    """
    grey = (np.clip(np.nan_to_num(array), 0, 1) * 255).astype(np.uint8)
    alpha = np.where(np.isnan(array), 0, 255).astype(np.uint8)
    buffer = io.BytesIO()
    Image.fromarray(np.stack([grey, alpha], axis=-1), "LA").save(buffer, "PNG")
    return "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode()


def residual_image(R):
    """Show residuals as brightness: mid-grey is 0, white is +4 sigma, black is -4 sigma."""
    return 0.5 + R / (8 * sigma)


def side_by_side(images):
    """Place d x d images in a row, separated by transparent pixels."""
    out = np.full((d, len(images) * (d + GAP) - GAP), np.nan)
    for j, image in enumerate(images):
        out[:, j * (d + GAP):j * (d + GAP) + d] = image
    return out


def label_images(fig, labels):
    """Put one label under each image of a side_by_side picture and hide the axes."""
    fig.update_xaxes(tickvals=[j * (d + GAP) + (d - 1) / 2 for j in range(len(labels))],
                     ticktext=labels, showgrid=False, zeroline=False, showline=False, ticks="")
    fig.update_yaxes(showticklabels=False, showgrid=False, zeroline=False, showline=False,
                     ticks="")
    return fig


def image_row(images, labels, height=300):
    """A static row of images with a label under each."""
    fig = go.Figure(go.Image(source=png_uri(side_by_side(images)), hoverinfo="skip"))
    style(fig, height=height, showlegend=False)
    label_images(fig, labels)
    fig.update_layout(margin=dict(l=20, r=20, t=20, b=40))
    return fig


def error_figure(lams, fold_errors, true_error, train_error, height=460):
    """The error on each fold (thin), their average CV(lambda) (thick), the true error,
    and the training error.

    Large dots mark the minimum of each averaged curve.
    """
    fig = go.Figure()
    for k, errors in enumerate(fold_errors):
        fig.add_trace(go.Scatter(
            x=lams, y=errors, mode="lines", name="error on each held-out fold",
            legendgroup="folds", showlegend=k == 0,
            line=dict(color=BLUE, width=1), opacity=0.3, hoverinfo="skip",
        ))
    curves = [(train_error, "training error", INK, "dash"),
              (fold_errors.mean(axis=0), "cross-validation error CV(λ)", BLUE, "solid"),
              (true_error, "true error (needs the sharp photo)", ORANGE, "solid")]
    for values, name, color, dash in curves:
        fig.add_trace(go.Scatter(
            x=lams, y=values, mode="lines", name=name, line=dict(color=color, width=2.5, dash=dash),
            hovertemplate="λ = %{x:.2e}<br>%{y:.3e}<extra></extra>",
        ))
        k = int(np.argmin(values))
        fig.add_trace(go.Scatter(
            x=[lams[k]], y=[values[k]], mode="markers", showlegend=False, hoverinfo="skip",
            marker=dict(color=color, size=13, line=dict(color="white", width=2)),
        ))
    fig.add_trace(go.Scatter(
        x=[lams[0], lams[-1]], y=[sigma**2, sigma**2], mode="lines", name="noise variance σ²",
        line=dict(color=INK, width=1, dash="dot"), hoverinfo="skip",
    ))
    style(fig, height=height)
    fig.update_xaxes(title="λ", type="log", exponentformat="power")
    fig.update_yaxes(title="mean squared error per pixel", type="log", exponentformat="power")
    return fig


def slider_figure(lams, steps, estimate, residuals, cv, true_error):
    """Sharp photo, the estimate, and its cross-validated residuals, with a slider for lambda.

    `estimate(j)` is the estimate fit to all the data with lambda = lams[j].
    """
    cv_best, best = int(np.argmin(cv)), int(np.argmin(true_error))

    def title(j):
        note = " (chosen by cross-validation)" if j == cv_best else (
            " (best possible)" if j == best else "")
        return (f"CV error {cv[j]:.3e}, RMS error of the estimate "
                f"{np.sqrt(true_error[j]):.3f}{note}")

    # The sharp photograph is drawn once; each step shows an estimate and its residuals.
    fig = go.Figure(go.Image(source=png_uri(photo), hoverinfo="skip"))
    for j in steps:
        fig.add_trace(go.Image(source=png_uri(estimate(j)), x0=d + GAP,
                               visible=j == cv_best, hoverinfo="skip"))
        fig.add_trace(go.Image(source=png_uri(residual_image(residuals[j])), x0=2 * (d + GAP),
                               visible=j == cv_best, hoverinfo="skip"))
    fig.update_layout(sliders=[dict(
        active=steps.index(cv_best), currentvalue=dict(prefix="λ = ", font=dict(color=INK)),
        font=dict(color="rgba(0,0,0,0)"),  # hide the crowded tick labels
        pad=dict(t=10), len=0.9, x=0.05,
        steps=[
            dict(method="update", label=f"{lams[j]:.1e}",
                 args=[{"visible": [True] + [i == j for i in steps for _ in range(2)]},
                       {"title.text": title(j)}])
            for j in steps
        ],
    )])
    style(fig, height=440, showlegend=False, title=title(cv_best))
    label_images(fig, ["sharp", "estimate", "cross-validated residuals"])
    fig.update_xaxes(range=[-0.5, 3 * d + 2 * GAP - 0.5])
    fig.update_layout(margin=dict(l=20, r=20, t=50, b=90))
    return fig

Our observations are the $d^2$ pixels of the blurred photograph. We hold out whole rows of it rather than scattered pixels, because then the training data still have the form “blur, then add noise”, and ridge regression can be fit to them unchanged. Figure 1 shows one of the $K$ folds.

Why whole rows?

Deleting the rows in $\Fcal_k$ from $\bfY$ deletes the same rows from the column blur, so the remaining data are \[ \bfY_{-k}=\bfX_{-k}\bfTheta\bfX^\top+\bfE_{-k}, \] where $\bfX_{-k}$ is $\bfX$ without the rows in $\Fcal_k$. This is the same kind of model with a shorter blur matrix, and it is fit the same way. The fit $\widehat{\bfTheta}^{(-k)}$ then predicts the held-out rows as $\bfX_k\widehat{\bfTheta}^{(-k)}\bfX^\top$, where $\bfX_k$ keeps only the rows in $\Fcal_k$, and the held-out residuals are $\bfY_k-\bfX_k\widehat{\bfTheta}^{(-k)}\bfX^\top$. Holding out scattered pixels instead would break the structure $\bfTheta\mapsto\bfX\bfTheta\bfX^\top$ that gives ridge regression its formula.

Split the rows into folds and plot one fold

K = 10
folds = np.array_split(rng.permutation(d), K)  # K random sets of rows of the photograph

held_out = np.isin(pixels, folds[0])[:, None] & np.ones(d, bool)
image_row([np.where(held_out, np.nan, blurred), np.where(held_out, blurred, np.nan)],
          ["training data", "held-out rows"]).show()

Figure 1: One of the folds. Left: the training data, the blurred photograph with the held-out rows removed. Right: the held-out rows, which the fit must predict.

For each $\lambda$, we fit to the training data, predict the held-out rows, and record the residuals. Doing this for all $K$ folds gives a residual at every pixel of the photograph, each one computed by a fit that never saw that pixel. We call this the cross-validated residual image; $\mathrm{CV}(\lambda)$ is the average over the folds of its mean square. Proposition 2 says what a good residual image looks like: when the fit predicts everything that can be predicted, what is left is the noise, so the residuals look like static with standard deviation $\sigma$.

Ridge regression

Flattening $\bfTheta$ into a vector in $\Rbb^{d^2}$ turns $\bfTheta\mapsto\bfX\bfTheta\bfX^\top$ into a $d^2\times d^2$ matrix with singular vectors $\bfv_i\bfv_j^\top$ and singular values $\sigma_i\sigma_j$. So by ridge regression filters the singular values, ridge regression for the photograph is \[ \widehat{\bfTheta}_\lambda=\sum_{i,j=1}^d\frac{(\sigma_i\sigma_j)^2}{(\sigma_i\sigma_j)^2+\lambda}\cdot\frac{\bfu_i^\top\bfY\bfu_j}{\sigma_i\sigma_j}\,\bfv_i\bfv_j^\top. \] Since this is a formula, ridge regression is cheap to fit, and we can afford a fine, logarithmically spaced grid of $\lambda$ values.

Cross-validate ridge regression

lams = np.logspace(-9, -1, 65)  # eight values of lambda per factor of ten
S2 = np.outer(s, s)             # entry (i, j) is the singular value sigma_i sigma_j


def ridge(Y, lam):
    """Ridge regression for the whole photograph Y, by the SVD formula."""
    return V @ ((U.T @ Y @ U) * S2 / (S2**2 + lam)) @ V.T


def ridge_true_error(Y):
    """The true error of ridge regression for every lambda (needs the sharp photo).

    V is orthogonal, so the error can be measured in the coordinates v_i v_j^T.
    """
    coeffs, target = U.T @ Y @ U, V.T @ photo @ V
    return np.array([np.mean((coeffs * S2 / (S2**2 + lam) - target)**2) for lam in lams])


def ridge_cv(Y, folds):
    """Cross-validate ridge regression on the blurred photograph Y.

    Returns the error on each fold for every lambda, a K x len(lams) array, and the
    cross-validated residual image for every lambda.
    """
    fold_errors = np.zeros((K, len(lams)))
    residuals = np.zeros((len(lams), d, d))
    for k, test in enumerate(folds):
        train = np.setdiff1d(pixels, test)
        # Ridge regression on the rows `train` is the same formula with the SVD of X[train]:
        # the estimate is V1 Z V^T, which predicts the held-out rows as (X[test] V1) Z (X V)^T.
        U1, s1, V1t = np.linalg.svd(X[train], full_matrices=False)
        S = np.outer(s1, s)
        coeffs = U1.T @ Y[train] @ U
        to_held_out = X[test] @ V1t.T
        for j, lam in enumerate(lams):
            Z = coeffs * S / (S**2 + lam)
            residuals[j, test] = Y[test] - to_held_out @ Z @ XV.T
            fold_errors[k, j] = np.mean(residuals[j, test]**2)
    return fold_errors, residuals


fold_errors, ridge_residuals = ridge_cv(blurred, folds)
cv = fold_errors.mean(axis=0)
true_error = ridge_true_error(blurred)
train_error = np.array([np.mean((blurred - X @ ridge(blurred, lam) @ X.T)**2) for lam in lams])
cv_best, best = int(np.argmin(cv)), int(np.argmin(true_error))

print(f"lambda chosen by cross-validation  = {lams[cv_best]:.1e}, "
      f"RMS error {np.sqrt(true_error[cv_best]):.3f} per pixel")
print(f"best lambda, using the sharp photo = {lams[best]:.1e}, "
      f"RMS error {np.sqrt(true_error[best]):.3f} per pixel")
print(f"for comparison, the blurred photo itself has RMS error "
      f"{np.sqrt(np.mean((blurred - photo)**2)):.3f}")

lambda chosen by cross-validation  = 5.6e-04, RMS error 0.068 per pixel
best lambda, using the sharp photo = 5.6e-03, RMS error 0.043 per pixel
for comparison, the blurred photo itself has RMS error 0.074

Figure 2 shows every curve that goes into the choice.

Plot the errors

error_figure(lams, fold_errors, true_error, train_error).show()

Figure 2: Judging lambda for ridge regression, per pixel, on a logarithmic scale. Grey, dashed: the training error. Thin blue: the error on each held-out fold. Thick blue: their average, the cross-validation error. Orange: the true error, which needs the sharp photograph. Large dots mark the minima. The dotted line is the noise variance sigma^2.

Three things stand out.

The training error does what Proposition 1 says: its minimum is at the left end of the grid, which is essentially least squares, whose true error is in the thousands. It even dips below $\sigma^2$, the sign of a fit that has absorbed some of the noise.
Each fold on its own gives a noisy verdict, but most folds have their minimum in the same neighbourhood, and averaging them makes a smooth curve with a clear minimum.
The minimum of the cross-validation error is in the right region, but to the left of the minimum of the true error.

To see what cross-validation is looking at, drag the slider in Figure 3. The middle image is the estimate, and the right image is its cross-validated residual image, with mid-grey for $0$, lighter for positive residuals, and darker for negative ones.

Build the ridge slider

steps = sorted(set(range(0, len(lams), 2)) | {cv_best, best})  # four values per factor of ten
slider_figure(lams, steps, lambda j: ridge(blurred, lams[j]), ridge_residuals, cv,
              true_error).show()

Figure 3: Left: the sharp photograph. Middle: the ridge estimate. Right: its cross-validated residual image, where mid-grey is 0 and lighter and darker pixels are positive and negative residuals. Drag the slider to change lambda; it starts at the value chosen by cross-validation.

Read the residuals from left to right.

$\lambda$ too small. The estimate is static. It predicts the held-out rows badly, because it fit the noise in the training rows, so the residuals are larger than the noise.
$\lambda$ chosen by cross-validation. The residuals look like pure static with no trace of the photograph: everything predictable has been predicted, and only the noise is left.
$\lambda$ too large. The estimate is too smooth and too dark, and the cat’s eyes and nose appear in the residuals. That is the bias: there is predictable structure in the data that the fit failed to capture.

Cross-validation picks the $\lambda$ that leaves the least behind.

Why it stops a little short

Move the slider from the cross-validation choice to the best $\lambda$, a few steps to the right. The estimate loses its mottled texture and gets visibly closer to the sharp photograph. But in the residuals, the only change is that a faint outline of the cat’s eyes and nose starts to appear. Cross-validation sees the bias arriving, but not the variance leaving. Proposition 2 explains why. Cross-validation measures how well we predict blurred data, and an error in $\widehat{\bfTheta}$ along $\bfv_i\bfv_j^\top$ shows up in the predicted data shrunk by $\sigma_i\sigma_j$. The mottled texture is made of exactly those poorly measured directions, so it costs the estimate a lot and the predictions almost nothing. That is also why the blue curves in Figure 2 are so much flatter than the orange one on the left.

Is this an accident of one noise draw and one split into folds? To find out, we repeat the whole experiment many times, each time with fresh noise and a fresh random split into folds. Each time we record the $\lambda$ that minimizes $\mathrm{CV}$, the $\lambda$ chosen by the one-standard-error rule, and the $\lambda$ that was actually best.

Repeat the experiment with fresh noise and fresh folds

def one_standard_error(fold_errors):
    """The largest lambda whose CV error is within one standard error of the minimum."""
    cv = fold_errors.mean(axis=0)
    se = fold_errors.std(axis=0, ddof=1) / np.sqrt(len(fold_errors))
    j = int(np.argmin(cv))
    return int(np.flatnonzero(cv <= cv[j] + se[j]).max())


rules = ["minimum of CV", "one-standard-error rule", "best possible"]
picks = {rule: [] for rule in rules}  # the index of the chosen lambda, in each repetition
rms_of = {rule: [] for rule in rules}  # its RMS error
repeat_rng = np.random.default_rng(2026)
for _ in range(50):
    Y = X @ photo @ X.T + sigma * repeat_rng.standard_normal((d, d))
    errors, _ = ridge_cv(Y, np.array_split(repeat_rng.permutation(d), K))
    trial_true = ridge_true_error(Y)
    for rule, j in zip(rules, [np.argmin(errors.mean(axis=0)), one_standard_error(errors),
                               np.argmin(trial_true)]):
        picks[rule].append(int(j))
        rms_of[rule].append(np.sqrt(trial_true[j]))
rms_of = {rule: np.array(values) for rule, values in rms_of.items()}

for rule in rules:
    print(f"{rule:24s} average RMS error {rms_of[rule].mean():.4f}, worst {rms_of[rule].max():.4f}")

minimum of CV            average RMS error 0.0572, worst 0.1158
one-standard-error rule  average RMS error 0.0479, worst 0.0579
best possible            average RMS error 0.0424, worst 0.0430

Plot the repetitions

colors = {"minimum of CV": BLUE, "one-standard-error rule": INK, "best possible": ORANGE}
fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.14,
                    subplot_titles=["which λ is chosen", "RMS error per pixel"])
for rule in rules:
    counts = np.bincount(picks[rule], minlength=len(lams))
    shown = counts > 0
    fig.add_trace(go.Bar(
        x=np.log10(lams[shown]), y=counts[shown], name=rule, marker_color=colors[rule],
        width=0.1, customdata=lams[shown], legendgroup=rule,
        hovertemplate="λ = %{customdata:.2e}<br>%{y} of 50<extra></extra>",
    ), row=1, col=1)
for rule in rules[:2]:
    fig.add_trace(go.Scatter(
        x=rms_of["best possible"], y=rms_of[rule], mode="markers", showlegend=False,
        legendgroup=rule,
        marker=dict(color=colors[rule], size=8, opacity=0.75, line=dict(color="white", width=1)),
        hovertemplate=f"{rule}<br>best λ: %{{x:.4f}}<br>this rule: %{{y:.4f}}<extra></extra>",
    ), row=1, col=2)
low = 0.95 * rms_of["best possible"].min()
high = 1.05 * max(rms_of[rule].max() for rule in rules)
fig.add_trace(go.Scatter(
    x=[low, high], y=[low, high], mode="lines", showlegend=False,
    line=dict(color=INK, width=1, dash="dash"), hoverinfo="skip",
), row=1, col=2)
style(fig, height=420)
fig.update_layout(barmode="group")
decades = np.arange(-4, 0)  # the bars are drawn at log10(lambda), labelled as powers of ten
fig.update_xaxes(title="λ", range=[-4, -1], tickvals=decades,
                 ticktext=[f"10<sup>{k}</sup>" for k in decades], row=1, col=1)
fig.update_yaxes(title="number of repetitions", row=1, col=1)
fig.update_xaxes(title="with the best λ", range=[low, high], row=1, col=2)
fig.update_yaxes(title="with the λ chosen from the data", range=[low, high], row=1, col=2)
fig.update_annotations(font=dict(color=INK, size=14))
fig.show()

Figure 4: Repetitions of the ridge experiment, each with fresh noise and fresh folds. Left: how often each lambda was chosen by the minimum of the cross-validation error (blue), by the one-standard-error rule (grey), and how often it was the best (orange). Right: for each repetition, the RMS error with the best lambda against the RMS error with the lambda chosen by each rule; the dashed line is where the two would be equal.

Most of the time the minimum of $\mathrm{CV}$ lands well to the left of the best $\lambda$, by the mechanism above. But the split into folds matters more than one might hope. In a few repetitions an unlucky split moves the minimum far to the left, and the estimate fills with static, sometimes ending up worse than the blurred photograph. The one-standard-error rule guards against exactly this. It never goes badly wrong, it moves the choice toward the best $\lambda$, and it gives a smaller RMS error on average and a much smaller one in the worst case. In every repetition, though, the best possible $\lambda$ stays out of reach: nothing computed from the blurred data alone can see the errors that the blur hides.

Summary

A hyperparameter is set before fitting, and the training error cannot choose it, because the training error always prefers the most flexible model.
$K$-fold cross-validation (alg. 1) judges each candidate by how well fits on part of the data predict the rest, and picks the best.
Cross-validation estimates the prediction error plus the noise variance (Proposition 2), not the error in $\bftheta$ itself. Errors along poorly measured directions are nearly invisible to it, so for an inverse problem like deblurring it tends to regularize too little.
It works for any method that can be fit to part of the data and used to predict the rest, and it can compare different methods as well as different values of $\lambda$.
The fold errors also say how uncertain $\mathrm{CV}(\lambda)$ is; the one-standard-error rule uses this to pick the most regularized $\lambda$ that is statistically tied with the minimum, which is more reliable from one data set to the next.
At a good choice of $\lambda$, the cross-validated residuals look like pure noise; structure left in them is bias, and residuals larger than the noise are variance.
Keep everything that depends on the data, including standardization, inside the folds, and keep a separate test set if you want an honest estimate of the final error.

References

Hastie, Trevor, Robert Tibshirani, and J. H. Friedman. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2nd ed. Springer Series in Statistics. Springer.

--- jupyter: python3 number-sections: false format: html: code-fold: true code-line-numbers: true code-tools: true out-width: "100%" execute: echo: true draft: false --- # Choosing Hyperparameters Most methods in this course come with a knob that the data do not set by themselves: the regularization parameter $\lambda$ in [ridge regression](../lectures/lecture-10.qmd#def-ridge), the number of principal components to keep, the degree of a polynomial model, the number of clusters. These notes are about how to set such a knob. The main tool is *cross-validation*, and we watch it choose the regularization parameter of ridge regression on a deblurring problem. ## Parameters and hyperparameters ::: {#def-hyperparameter} ## Hyperparameter A *hyperparameter* is a number that must be fixed before a method can be fit to data, as opposed to a *parameter*, which the fitting itself determines. ::: In ridge regression, $$ \widehat{\bftheta}_\lambda=\argmin_{\bftheta\in\Rbb^d}\Big\{\|\bfy-\bfX\bftheta\|_2^2+\lambda\|\bftheta\|_2^2\Big\}, $$ the coefficients $\bftheta$ are parameters and $\lambda$ is a hyperparameter: once $\lambda$ is fixed, [the ridge solution](../lectures/lecture-10.qmd#thm-ridge-solution) determines $\widehat{\bftheta}_\lambda$, but nothing in the objective tells us which $\lambda$ to use. [Lecture 10](../lectures/lecture-10.qmd#prp-ridge-bias-variance) tells us what the best $\lambda$ would be. Under the [standard assumptions](../lectures/lecture-08.qmd#def-standard-assumptions), the mean squared error is $$ \operatorname{MSE}(\widehat{\bftheta}_\lambda) =\sum_{i=1}^d\frac{\lambda^2\alpha_i^2}{(\sigma_i^2+\lambda)^2} +\sigma^2\sum_{i=1}^d\frac{\sigma_i^2}{(\sigma_i^2+\lambda)^2}, \qquad \alpha_i=\bfv_i^\top\bftheta, $$ and we would choose $\lambda$ to make this as small as possible. But it depends on the true $\bftheta$ and the noise level $\sigma^2$, which we do not know. We need a yardstick that can be computed from the data alone. ## Training error is the wrong yardstick The first idea is to choose the $\lambda$ whose fit matches the data best. This never works. Throughout, $\bfX\in\Rbb^{n\times d}$ has rank $r\geq1$ and [compact SVD](../lectures/lecture-04.qmd#thm-svd) $\bfX=\sum_{i=1}^r\sigma_i\bfu_i\bfv_i^\top$. ::: {#prp-ridge-training-error} ## Training error always prefers $\lambda=0$ The training residual of ridge regression satisfies $$ \|\bfy-\bfX\widehat{\bftheta}_\lambda\|_2^2 =\sum_{i=1}^r\Big(\frac{\lambda}{\sigma_i^2+\lambda}\Big)^2(\bfu_i^\top\bfy)^2+\|\bfy-\bfP\bfy\|_2^2, $$ where $\bfP=\bfU_r\bfU_r^\top$ is the orthogonal projection onto $\Ccal(\bfX)$. In particular, it is an increasing function of $\lambda$. ::: ::: {.proof} By [ridge regression filters the singular values](../lectures/lecture-10.qmd#thm-ridge-svd) and $\bfX\bfv_i=\sigma_i\bfu_i$, the fitted values are $\bfX\widehat{\bftheta}_\lambda=\sum_{i=1}^r\frac{\sigma_i^2}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\bfu_i$. Since $\bfy=\sum_{i=1}^r(\bfu_i^\top\bfy)\bfu_i+(\bfy-\bfP\bfy)$, $$ \bfy-\bfX\widehat{\bftheta}_\lambda =\sum_{i=1}^r\frac{\lambda}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\bfu_i+(\bfy-\bfP\bfy). $$ The vectors $\bfu_1,\ldots,\bfu_r$ are orthonormal and $\bfy-\bfP\bfy\perp\Ccal(\bfX)$, so the [Pythagorean theorem](../lectures/lecture-05.qmd#thm-pythagoras) gives the formula. Each $\lambda/(\sigma_i^2+\lambda)$ increases with $\lambda$. ::: Of course: least squares is *defined* as the minimizer of the training error, so anything else does worse on the training data. The same trap appears for every hyperparameter that controls how flexible a model is. More principal components, a higher polynomial degree, or more clusters always fit the training data at least as well, so the training error always asks for the most flexible model, including its ability to fit the noise. To judge a hyperparameter honestly, we must test the fit on data that was not used to compute it. ## Validation and cross-validation The simplest option is to hold out a *validation set*: fit on most of the data, measure the error on the rest, and pick the $\lambda$ with the smallest validation error. This wastes data, and the answer depends on which points happened to be held out. Cross-validation fixes both problems by rotating the role of validation set through the data. ::: {#alg-k-fold-cv} ## $K$-fold cross-validation Input: data $(\bfX,\bfy)$, a method that fits an estimate $\widehat{\bftheta}_\lambda$ to data for each value of a hyperparameter $\lambda$, a grid of candidate values $\lambda_1,\ldots,\lambda_m$, and a number of folds $K$ (typically $5$ or $10$). 1. Randomly shuffle the indices $\{1,\ldots,n\}$ and split them into $K$ disjoint *folds* $\Fcal_1,\ldots,\Fcal_K$ of nearly equal size. 2. For each candidate $\lambda$ and each fold $k=1,\ldots,K$: a. Fit $\widehat{\bftheta}_\lambda^{(-k)}$ using only the observations *not* in $\Fcal_k$. b. Compute its error on the held-out fold, $\displaystyle\mathrm{Err}_k(\lambda)=\frac{1}{|\Fcal_k|}\sum_{i\in\Fcal_k}\big(y_i-\bfx_i^\top\widehat{\bftheta}_\lambda^{(-k)}\big)^2$. 3. Average over the folds: $\displaystyle\mathrm{CV}(\lambda)=\frac1K\sum_{k=1}^K\mathrm{Err}_k(\lambda)$. 4. Choose $\lambda^\ast=\argmin_{\lambda_j}\mathrm{CV}(\lambda_j)$, and refit on *all* the data with $\lambda^\ast$. ::: Every observation is used for validation exactly once and for fitting $K-1$ times. The choice $K=n$ is called *leave-one-out* cross-validation. Nothing in the algorithm needs a formula for the estimate, only a way to fit it to part of the data and to predict the rest, so it works for any method. It can even choose between methods, if we treat the method itself as one more hyperparameter. For ridge regression the grid of $\lambda$ values is usually logarithmically spaced, say $10^{-3},\ldots,10^{3}$, since $\lambda$ only matters relative to the $\sigma_i^2$. $\mathrm{CV}(\lambda)$ is itself an average of $K$ noisy numbers, so its minimum is not exact. The spread of the fold errors measures how uncertain it is, and suggests a more cautious choice [@hastie_elements_2009, Section 7.10]. ::: {.remark name="The one-standard-error rule"} Let $\mathrm{SE}(\lambda)$ be the sample standard deviation of $\mathrm{Err}_1(\lambda),\ldots,\mathrm{Err}_K(\lambda)$, divided by $\sqrt{K}$, and let $\lambda^\ast$ minimize $\mathrm{CV}$. Values of $\lambda$ whose $\mathrm{CV}(\lambda)$ is within $\mathrm{SE}(\lambda^\ast)$ of the minimum cannot really be told apart from $\lambda^\ast$. Among them, choose the one that regularizes most, which for ridge regression is the largest. ::: ::: {.callout-warning} ## Keep the held-out data held out Anything learned from the data counts as fitting, including the means and standard deviations used to standardize the features. Compute them from the training folds only, then apply them to the held-out fold. And if you want an honest estimate of the final model's error, you need a separate *test set* that played no role in choosing $\lambda$ at all. ::: ::: {.callout-tip} ## Cross-validation for ridge regression in NumPy Using the same folds for every $\lambda$ makes the comparison fair. Solving with $\bfX^\top\bfX+\lambda\bfI_d$ directly is acceptable here, because $\lambda>0$ keeps its condition number under control. ```python import numpy as np def ridge(X, y, lam): d = X.shape[1] return np.linalg.solve(X.T @ X + lam * np.eye(d), X.T @ y) def cv_error(X, y, lam, K=5, seed=0): rng = np.random.default_rng(seed) folds = np.array_split(rng.permutation(len(y)), K) errors = [] for k in range(K): test = folds[k] train = np.concatenate([folds[j] for j in range(K) if j != k]) theta = ridge(X[train], y[train], lam) errors.append(np.mean((y[test] - X[test] @ theta) ** 2)) return np.mean(errors) lambdas = np.logspace(-3, 3, 50) best_lam = min(lambdas, key=lambda lam: cv_error(X, y, lam)) ``` In scikit-learn, `RidgeCV(alphas=lambdas, cv=5)` does the same search, and `GridSearchCV` does it for any model and any hyperparameters. ::: ## What cross-validation measures Cross-validation scores a fit by how well it predicts new observations. That is not the same as how close $\widehat{\bftheta}$ is to $\bftheta$, and the difference is worth making precise. ::: {#prp-prediction-error} ## Prediction error Let $\widehat{\bftheta}$ be any estimate of $\bftheta$, and let $y=\bfx^\top\bftheta+\varepsilon$ be a new observation at a fixed $\bfx\in\Rbb^d$, whose noise $\varepsilon$ has mean $0$ and variance $\sigma^2$ and is independent of $\widehat{\bftheta}$. Then $$ \Ebb\big[(y-\bfx^\top\widehat{\bftheta})^2\big] =\sigma^2+\Ebb\big[\big(\bfx^\top(\widehat{\bftheta}-\bftheta)\big)^2\big]. $$ ::: ::: {.proof} Write $y-\bfx^\top\widehat{\bftheta}=\varepsilon-\bfx^\top(\widehat{\bftheta}-\bftheta)$ and expand the square: $$ \Ebb\big[(y-\bfx^\top\widehat{\bftheta})^2\big] =\Ebb[\varepsilon^2]-2\,\Ebb\big[\varepsilon\,\bfx^\top(\widehat{\bftheta}-\bftheta)\big]+\Ebb\big[\big(\bfx^\top(\widehat{\bftheta}-\bftheta)\big)^2\big]. $$ The first term is $\Var(\varepsilon)=\sigma^2$. By independence, the middle expectation is $\Ebb[\varepsilon]\,\Ebb[\bfx^\top(\widehat{\bftheta}-\bftheta)]=0$. ::: In $K$-fold cross-validation, the noise in the held-out fold is independent of the fit $\widehat{\bftheta}_\lambda^{(-k)}$, which never saw it. So $\mathrm{CV}(\lambda)$ estimates $\sigma^2$, which no choice of $\lambda$ can change, plus the average *prediction error* $(\bfx_i^\top(\widehat{\bftheta}-\bftheta))^2$ over the held-out rows. In particular, the cross-validation error cannot go much below $\sigma^2$. The prediction error weighs the error in each direction by how well the data see that direction. Averaging over all the rows of $\bfX$ and writing $\widehat{\bftheta}-\bftheta$ in the basis $\bfv_1,\ldots,\bfv_d$ (assuming (A3), so $r=d$), $$ \frac1n\|\bfX(\widehat{\bftheta}-\bftheta)\|_2^2 =\frac1n\sum_{i=1}^d\sigma_i^2\big(\bfv_i^\top(\widehat{\bftheta}-\bftheta)\big)^2, \qquad\text{while}\qquad \|\widehat{\bftheta}-\bftheta\|_2^2 =\sum_{i=1}^d\big(\bfv_i^\top(\widehat{\bftheta}-\bftheta)\big)^2. $$ An error along a poorly measured direction, with $\sigma_i$ small, is nearly invisible to cross-validation, even when it dominates $\|\widehat{\bftheta}-\bftheta\|_2^2$. When the goal is prediction, this is exactly right: an error the data cannot see will not hurt future predictions either. When the goal is $\bftheta$ itself, as in the example below, cross-validation tends to be too lenient toward those errors. ## Example: deblurring a photograph We return to the out-of-focus camera from [the variance of least squares](../examples/least-squares-variance/least-squares-variance.qmd), this time pointed at a cat. A greyscale photograph $\bfTheta$ of Chelsea, the cat that ships with scikit-image, is blurred by the same symmetric matrix $\bfX$, whose row $j$ is a narrow Gaussian bump centred on pixel $j$, and the camera adds a little noise. We observe $$ \bfY=\bfX\bfTheta\bfX^\top+\bfE, \qquad E_{pq}\ \text{independent},\ \mathcal N(0,\sigma^2), $$ where the noise level $\sigma$ is small compared to the range from black to white. The smallest singular values of the blur are many orders of magnitude below the largest, so least squares turns the noise into static, and we must regularize. The question is how much. This example has one luxury that real problems do not: we have the sharp photograph. So for every $\lambda$ we can compute the *true error*, the mean squared difference per pixel between the estimate and the sharp photograph, and grade the choice that cross-validation makes without ever showing it the answer. The code is folded throughout; the figures are the point. ```{python} #| code-summary: "Set up the blurred photograph" import base64 import io import numpy as np import plotly.graph_objects as go from matplotlib import colormaps from PIL import Image from plotly.subplots import make_subplots from skimage import color, data from threadpoolctl import threadpool_limits threadpool_limits(1) # this page does many small matrix products, which run fastest on one thread rng = np.random.default_rng(123) # The blur: row j of X is a Gaussian bump of width w centred on pixel j. d = 100 w = 2.0 pixels = np.arange(d) X = np.exp(-(pixels[:, None] - pixels[None, :])**2 / (2 * w**2)) / (np.sqrt(2 * np.pi) * w) sigma = 0.01 # the sensor noise: 1% of the range from black to white # scikit-image's photograph of Chelsea the cat in greyscale, cropped from 300 x 451 to a # 300 x 300 square around her face, then averaged down to 100 x 100 to keep the page fast. photo = color.rgb2gray(data.chelsea())[:, 90:390].reshape(d, 3, d, 3).mean(axis=(1, 3)) # Blur every column, then every row, then add sensor noise to every pixel. blurred = X @ photo @ X.T + sigma * rng.standard_normal((d, d)) U, s, Vt = np.linalg.svd(X) V = Vt.T XV = X @ V # equals U diag(s) ``` ```{python} #| code-summary: "Figure style and image helpers" # The same muted palette as the examples, which reads correctly against # both the light and the dark version of this site. INK = "#7c7c78" GRID = "rgba(128, 128, 128, 0.28)" BLUE, ORANGE = "#2a78d6", "#eb6834" GAP = 6 # transparent pixels between images placed side by side def style(fig, height=420, title=None, showlegend=True): """Apply the shared layout to a plot with axes.""" fig.update_layout( template="simple_white", paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)", font=dict(color=INK, size=13), title=dict(text=title, x=0, xanchor="left", font=dict(size=15)), margin=dict(l=60, r=30, t=60 if title else 40, b=50), height=height, showlegend=showlegend, hovermode="closest", legend=dict(orientation="h", y=-0.18, x=0), ) fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID, showgrid=True) fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID, showgrid=True) return fig def png_uri(array): """Encode a greyscale array as a PNG data URI, clipping brightness to [0, 1]. NaN pixels become transparent. A PNG is far smaller in the page than a heatmap. """ grey = (np.clip(np.nan_to_num(array), 0, 1) * 255).astype(np.uint8) alpha = np.where(np.isnan(array), 0, 255).astype(np.uint8) buffer = io.BytesIO() Image.fromarray(np.stack([grey, alpha], axis=-1), "LA").save(buffer, "PNG") return "data:image/png;base64," + base64.b64encode(buffer.getvalue()).decode() def residual_image(R): """Show residuals as brightness: mid-grey is 0, white is +4 sigma, black is -4 sigma.""" return 0.5 + R / (8 * sigma) def side_by_side(images): """Place d x d images in a row, separated by transparent pixels.""" out = np.full((d, len(images) * (d + GAP) - GAP), np.nan) for j, image in enumerate(images): out[:, j * (d + GAP):j * (d + GAP) + d] = image return out def label_images(fig, labels): """Put one label under each image of a side_by_side picture and hide the axes.""" fig.update_xaxes(tickvals=[j * (d + GAP) + (d - 1) / 2 for j in range(len(labels))], ticktext=labels, showgrid=False, zeroline=False, showline=False, ticks="") fig.update_yaxes(showticklabels=False, showgrid=False, zeroline=False, showline=False, ticks="") return fig def image_row(images, labels, height=300): """A static row of images with a label under each.""" fig = go.Figure(go.Image(source=png_uri(side_by_side(images)), hoverinfo="skip")) style(fig, height=height, showlegend=False) label_images(fig, labels) fig.update_layout(margin=dict(l=20, r=20, t=20, b=40)) return fig def error_figure(lams, fold_errors, true_error, train_error, height=460): """The error on each fold (thin), their average CV(lambda) (thick), the true error, and the training error. Large dots mark the minimum of each averaged curve. """ fig = go.Figure() for k, errors in enumerate(fold_errors): fig.add_trace(go.Scatter( x=lams, y=errors, mode="lines", name="error on each held-out fold", legendgroup="folds", showlegend=k == 0, line=dict(color=BLUE, width=1), opacity=0.3, hoverinfo="skip", )) curves = [(train_error, "training error", INK, "dash"), (fold_errors.mean(axis=0), "cross-validation error CV(λ)", BLUE, "solid"), (true_error, "true error (needs the sharp photo)", ORANGE, "solid")] for values, name, color, dash in curves: fig.add_trace(go.Scatter( x=lams, y=values, mode="lines", name=name, line=dict(color=color, width=2.5, dash=dash), hovertemplate="λ = %{x:.2e} %{y:.3e}<extra></extra>", )) k = int(np.argmin(values)) fig.add_trace(go.Scatter( x=[lams[k]], y=[values[k]], mode="markers", showlegend=False, hoverinfo="skip", marker=dict(color=color, size=13, line=dict(color="white", width=2)), )) fig.add_trace(go.Scatter( x=[lams[0], lams[-1]], y=[sigma**2, sigma**2], mode="lines", name="noise variance σ²", line=dict(color=INK, width=1, dash="dot"), hoverinfo="skip", )) style(fig, height=height) fig.update_xaxes(title="λ", type="log", exponentformat="power") fig.update_yaxes(title="mean squared error per pixel", type="log", exponentformat="power") return fig def slider_figure(lams, steps, estimate, residuals, cv, true_error): """Sharp photo, the estimate, and its cross-validated residuals, with a slider for lambda. `estimate(j)` is the estimate fit to all the data with lambda = lams[j]. """ cv_best, best = int(np.argmin(cv)), int(np.argmin(true_error)) def title(j): note = " (chosen by cross-validation)" if j == cv_best else ( " (best possible)" if j == best else "") return (f"CV error {cv[j]:.3e}, RMS error of the estimate " f"{np.sqrt(true_error[j]):.3f}{note}") # The sharp photograph is drawn once; each step shows an estimate and its residuals. fig = go.Figure(go.Image(source=png_uri(photo), hoverinfo="skip")) for j in steps: fig.add_trace(go.Image(source=png_uri(estimate(j)), x0=d + GAP, visible=j == cv_best, hoverinfo="skip")) fig.add_trace(go.Image(source=png_uri(residual_image(residuals[j])), x0=2 * (d + GAP), visible=j == cv_best, hoverinfo="skip")) fig.update_layout(sliders=[dict( active=steps.index(cv_best), currentvalue=dict(prefix="λ = ", font=dict(color=INK)), font=dict(color="rgba(0,0,0,0)"), # hide the crowded tick labels pad=dict(t=10), len=0.9, x=0.05, steps=[ dict(method="update", label=f"{lams[j]:.1e}", args=[{"visible": [True] + [i == j for i in steps for _ in range(2)]}, {"title.text": title(j)}]) for j in steps ], )]) style(fig, height=440, showlegend=False, title=title(cv_best)) label_images(fig, ["sharp", "estimate", "cross-validated residuals"]) fig.update_xaxes(range=[-0.5, 3 * d + 2 * GAP - 0.5]) fig.update_layout(margin=dict(l=20, r=20, t=50, b=90)) return fig ``` ### Holding out rows Our observations are the $d^2$ pixels of the blurred photograph. We hold out whole rows of it rather than scattered pixels, because then the training data still have the form "blur, then add noise", and ridge regression can be fit to them unchanged. @fig-fold shows one of the $K$ folds. ::: {.callout-note collapse="true"} ## Why whole rows? Deleting the rows in $\Fcal_k$ from $\bfY$ deletes the same rows from the column blur, so the remaining data are $$ \bfY_{-k}=\bfX_{-k}\bfTheta\bfX^\top+\bfE_{-k}, $$ where $\bfX_{-k}$ is $\bfX$ without the rows in $\Fcal_k$. This is the same kind of model with a shorter blur matrix, and it is fit the same way. The fit $\widehat{\bfTheta}^{(-k)}$ then predicts the held-out rows as $\bfX_k\widehat{\bfTheta}^{(-k)}\bfX^\top$, where $\bfX_k$ keeps only the rows in $\Fcal_k$, and the held-out residuals are $\bfY_k-\bfX_k\widehat{\bfTheta}^{(-k)}\bfX^\top$. Holding out scattered pixels instead would break the structure $\bfTheta\mapsto\bfX\bfTheta\bfX^\top$ that gives ridge regression its formula. ::: ```{python} #| label: fig-fold #| fig-cap: "One of the folds. Left: the training data, the blurred photograph with the held-out rows removed. Right: the held-out rows, which the fit must predict." #| fig-alt: "Two greyscale images side by side. The left one is the blurred photograph of the cat with a scattering of thin horizontal rows missing at random heights. The right one shows only those rows, as thin horizontal strips on an empty background." #| code-summary: "Split the rows into folds and plot one fold" K = 10 folds = np.array_split(rng.permutation(d), K) # K random sets of rows of the photograph held_out = np.isin(pixels, folds[0])[:, None] & np.ones(d, bool) image_row([np.where(held_out, np.nan, blurred), np.where(held_out, blurred, np.nan)], ["training data", "held-out rows"]).show() ``` For each $\lambda$, we fit to the training data, predict the held-out rows, and record the residuals. Doing this for all $K$ folds gives a residual at every pixel of the photograph, each one computed by a fit that never saw that pixel. We call this the *cross-validated residual image*; $\mathrm{CV}(\lambda)$ is the average over the folds of its mean square. @prp-prediction-error says what a good residual image looks like: when the fit predicts everything that can be predicted, what is left is the noise, so the residuals look like static with standard deviation $\sigma$. ### Ridge regression Flattening $\bfTheta$ into a vector in $\Rbb^{d^2}$ turns $\bfTheta\mapsto\bfX\bfTheta\bfX^\top$ into a $d^2\times d^2$ matrix with singular vectors $\bfv_i\bfv_j^\top$ and singular values $\sigma_i\sigma_j$. So by [ridge regression filters the singular values](../lectures/lecture-10.qmd#thm-ridge-svd), ridge regression for the photograph is $$ \widehat{\bfTheta}_\lambda=\sum_{i,j=1}^d\frac{(\sigma_i\sigma_j)^2}{(\sigma_i\sigma_j)^2+\lambda}\cdot\frac{\bfu_i^\top\bfY\bfu_j}{\sigma_i\sigma_j}\,\bfv_i\bfv_j^\top. $$ Since this is a formula, ridge regression is cheap to fit, and we can afford a fine, logarithmically spaced grid of $\lambda$ values. ```{python} #| code-summary: "Cross-validate ridge regression" lams = np.logspace(-9, -1, 65) # eight values of lambda per factor of ten S2 = np.outer(s, s) # entry (i, j) is the singular value sigma_i sigma_j def ridge(Y, lam): """Ridge regression for the whole photograph Y, by the SVD formula.""" return V @ ((U.T @ Y @ U) * S2 / (S2**2 + lam)) @ V.T def ridge_true_error(Y): """The true error of ridge regression for every lambda (needs the sharp photo). V is orthogonal, so the error can be measured in the coordinates v_i v_j^T. """ coeffs, target = U.T @ Y @ U, V.T @ photo @ V return np.array([np.mean((coeffs * S2 / (S2**2 + lam) - target)**2) for lam in lams]) def ridge_cv(Y, folds): """Cross-validate ridge regression on the blurred photograph Y. Returns the error on each fold for every lambda, a K x len(lams) array, and the cross-validated residual image for every lambda. """ fold_errors = np.zeros((K, len(lams))) residuals = np.zeros((len(lams), d, d)) for k, test in enumerate(folds): train = np.setdiff1d(pixels, test) # Ridge regression on the rows `train` is the same formula with the SVD of X[train]: # the estimate is V1 Z V^T, which predicts the held-out rows as (X[test] V1) Z (X V)^T. U1, s1, V1t = np.linalg.svd(X[train], full_matrices=False) S = np.outer(s1, s) coeffs = U1.T @ Y[train] @ U to_held_out = X[test] @ V1t.T for j, lam in enumerate(lams): Z = coeffs * S / (S**2 + lam) residuals[j, test] = Y[test] - to_held_out @ Z @ XV.T fold_errors[k, j] = np.mean(residuals[j, test]**2) return fold_errors, residuals fold_errors, ridge_residuals = ridge_cv(blurred, folds) cv = fold_errors.mean(axis=0) true_error = ridge_true_error(blurred) train_error = np.array([np.mean((blurred - X @ ridge(blurred, lam) @ X.T)**2) for lam in lams]) cv_best, best = int(np.argmin(cv)), int(np.argmin(true_error)) print(f"lambda chosen by cross-validation = {lams[cv_best]:.1e}, " f"RMS error {np.sqrt(true_error[cv_best]):.3f} per pixel") print(f"best lambda, using the sharp photo = {lams[best]:.1e}, " f"RMS error {np.sqrt(true_error[best]):.3f} per pixel") print(f"for comparison, the blurred photo itself has RMS error " f"{np.sqrt(np.mean((blurred - photo)**2)):.3f}") ``` @fig-ridge-cv shows every curve that goes into the choice. ```{python} #| label: fig-ridge-cv #| fig-cap: "Judging lambda for ridge regression, per pixel, on a logarithmic scale. Grey, dashed: the training error. Thin blue: the error on each held-out fold. Thick blue: their average, the cross-validation error. Orange: the true error, which needs the sharp photograph. Large dots mark the minima. The dotted line is the noise variance sigma^2." #| fig-alt: "A chart against lambda, with both axes logarithmic. A grey dashed curve, the training error, rises steadily from below the dotted noise line at the left, so its minimum is at the far left. An orange curve, the true error, falls steeply from enormous values at the left to a minimum near the right of the chart, then rises a little. Thin blue curves, one per fold, form a narrow band; each falls gently to a minimum and rises sharply after that, and most of their minima are close together. A thick blue curve runs through the middle of the band, with its minimum just above the dotted line, well to the left of the orange minimum." #| code-summary: "Plot the errors" error_figure(lams, fold_errors, true_error, train_error).show() ``` Three things stand out. - The training error does what @prp-ridge-training-error says: its minimum is at the left end of the grid, which is essentially least squares, whose true error is in the thousands. It even dips below $\sigma^2$, the sign of a fit that has absorbed some of the noise. - Each fold on its own gives a noisy verdict, but most folds have their minimum in the same neighbourhood, and averaging them makes a smooth curve with a clear minimum. - The minimum of the cross-validation error is in the right region, but to the left of the minimum of the true error. To see what cross-validation is looking at, drag the slider in @fig-ridge-slider. The middle image is the estimate, and the right image is its cross-validated residual image, with mid-grey for $0$, lighter for positive residuals, and darker for negative ones. ```{python} #| label: fig-ridge-slider #| fig-cap: "Left: the sharp photograph. Middle: the ridge estimate. Right: its cross-validated residual image, where mid-grey is 0 and lighter and darker pixels are positive and negative residuals. Drag the slider to change lambda; it starts at the value chosen by cross-validation." #| fig-alt: "An interactive figure with three greyscale images side by side and a slider for lambda. The left image is the sharp photograph of the cat. For small lambda, the middle image, the ridge estimate, is static, and the right image, the residuals, is strong static with horizontal streaks. At the lambda chosen by cross-validation, the estimate is a sharpened but mottled photograph and the residuals are even static with no visible structure. For large lambda, the estimate is blurry and the outlines of the cat's eyes and nose appear in the residuals." #| code-summary: "Build the ridge slider" steps = sorted(set(range(0, len(lams), 2)) | {cv_best, best}) # four values per factor of ten slider_figure(lams, steps, lambda j: ridge(blurred, lams[j]), ridge_residuals, cv, true_error).show() ``` Read the residuals from left to right. - **$\lambda$ too small.** The estimate is static. It predicts the held-out rows badly, because it fit the noise in the training rows, so the residuals are larger than the noise. - **$\lambda$ chosen by cross-validation.** The residuals look like pure static with no trace of the photograph: everything predictable has been predicted, and only the noise is left. - **$\lambda$ too large.** The estimate is too smooth and too dark, and the cat's eyes and nose appear in the residuals. That is the bias: there is predictable structure in the data that the fit failed to capture. Cross-validation picks the $\lambda$ that leaves the least behind. ### Why it stops a little short Move the slider from the cross-validation choice to the best $\lambda$, a few steps to the right. The estimate loses its mottled texture and gets visibly closer to the sharp photograph. But in the residuals, the only change is that a faint outline of the cat's eyes and nose starts to appear. Cross-validation sees the bias arriving, but not the variance leaving. @prp-prediction-error explains why. Cross-validation measures how well we predict *blurred* data, and an error in $\widehat{\bfTheta}$ along $\bfv_i\bfv_j^\top$ shows up in the predicted data shrunk by $\sigma_i\sigma_j$. The mottled texture is made of exactly those poorly measured directions, so it costs the estimate a lot and the predictions almost nothing. That is also why the blue curves in @fig-ridge-cv are so much flatter than the orange one on the left. Is this an accident of one noise draw and one split into folds? To find out, we repeat the whole experiment many times, each time with fresh noise and a fresh random split into folds. Each time we record the $\lambda$ that minimizes $\mathrm{CV}$, the $\lambda$ chosen by the one-standard-error rule, and the $\lambda$ that was actually best. ```{python} #| code-summary: "Repeat the experiment with fresh noise and fresh folds" def one_standard_error(fold_errors): """The largest lambda whose CV error is within one standard error of the minimum.""" cv = fold_errors.mean(axis=0) se = fold_errors.std(axis=0, ddof=1) / np.sqrt(len(fold_errors)) j = int(np.argmin(cv)) return int(np.flatnonzero(cv <= cv[j] + se[j]).max()) rules = ["minimum of CV", "one-standard-error rule", "best possible"] picks = {rule: [] for rule in rules} # the index of the chosen lambda, in each repetition rms_of = {rule: [] for rule in rules} # its RMS error repeat_rng = np.random.default_rng(2026) for _ in range(50): Y = X @ photo @ X.T + sigma * repeat_rng.standard_normal((d, d)) errors, _ = ridge_cv(Y, np.array_split(repeat_rng.permutation(d), K)) trial_true = ridge_true_error(Y) for rule, j in zip(rules, [np.argmin(errors.mean(axis=0)), one_standard_error(errors), np.argmin(trial_true)]): picks[rule].append(int(j)) rms_of[rule].append(np.sqrt(trial_true[j])) rms_of = {rule: np.array(values) for rule, values in rms_of.items()} for rule in rules: print(f"{rule:24s} average RMS error {rms_of[rule].mean():.4f}, worst {rms_of[rule].max():.4f}") ``` ```{python} #| label: fig-ridge-repeat #| fig-cap: "Repetitions of the ridge experiment, each with fresh noise and fresh folds. Left: how often each lambda was chosen by the minimum of the cross-validation error (blue), by the one-standard-error rule (grey), and how often it was the best (orange). Right: for each repetition, the RMS error with the best lambda against the RMS error with the lambda chosen by each rule; the dashed line is where the two would be equal." #| fig-alt: "Two charts. On the left, a bar chart against lambda on a logarithmic scale. Blue bars, the minima of the cross-validation error, are bunched together, with a thin tail reaching further left. Grey bars, the one-standard-error rule, sit a little to the right of the blue ones. Orange bars, the best values, sit further right still, in a narrow cluster. On the right, the RMS error with the best lambda, nearly the same in every repetition, against the RMS error with the chosen lambda. Grey dots cluster above the dashed diagonal. Blue dots lie higher, and a few lie far above." #| code-summary: "Plot the repetitions" colors = {"minimum of CV": BLUE, "one-standard-error rule": INK, "best possible": ORANGE} fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.14, subplot_titles=["which λ is chosen", "RMS error per pixel"]) for rule in rules: counts = np.bincount(picks[rule], minlength=len(lams)) shown = counts > 0 fig.add_trace(go.Bar( x=np.log10(lams[shown]), y=counts[shown], name=rule, marker_color=colors[rule], width=0.1, customdata=lams[shown], legendgroup=rule, hovertemplate="λ = %{customdata:.2e} %{y} of 50<extra></extra>", ), row=1, col=1) for rule in rules[:2]: fig.add_trace(go.Scatter( x=rms_of["best possible"], y=rms_of[rule], mode="markers", showlegend=False, legendgroup=rule, marker=dict(color=colors[rule], size=8, opacity=0.75, line=dict(color="white", width=1)), hovertemplate=f"{rule} best λ: %{{x:.4f}} this rule: %{{y:.4f}}<extra></extra>", ), row=1, col=2) low = 0.95 * rms_of["best possible"].min() high = 1.05 * max(rms_of[rule].max() for rule in rules) fig.add_trace(go.Scatter( x=[low, high], y=[low, high], mode="lines", showlegend=False, line=dict(color=INK, width=1, dash="dash"), hoverinfo="skip", ), row=1, col=2) style(fig, height=420) fig.update_layout(barmode="group") decades = np.arange(-4, 0) # the bars are drawn at log10(lambda), labelled as powers of ten fig.update_xaxes(title="λ", range=[-4, -1], tickvals=decades, ticktext=[f"10{k}" for k in decades], row=1, col=1) fig.update_yaxes(title="number of repetitions", row=1, col=1) fig.update_xaxes(title="with the best λ", range=[low, high], row=1, col=2) fig.update_yaxes(title="with the λ chosen from the data", range=[low, high], row=1, col=2) fig.update_annotations(font=dict(color=INK, size=14)) fig.show() ``` Most of the time the minimum of $\mathrm{CV}$ lands well to the left of the best $\lambda$, by the mechanism above. But the split into folds matters more than one might hope. In a few repetitions an unlucky split moves the minimum far to the left, and the estimate fills with static, sometimes ending up worse than the blurred photograph. The one-standard-error rule guards against exactly this. It never goes badly wrong, it moves the choice toward the best $\lambda$, and it gives a smaller RMS error on average and a much smaller one in the worst case. In every repetition, though, the best possible $\lambda$ stays out of reach: nothing computed from the blurred data alone can see the errors that the blur hides. ## Summary - A hyperparameter is set before fitting, and the training error cannot choose it, because the training error always prefers the most flexible model. - $K$-fold cross-validation (@alg-k-fold-cv) judges each candidate by how well fits on part of the data predict the rest, and picks the best. - Cross-validation estimates the prediction error plus the noise variance (@prp-prediction-error), not the error in $\bftheta$ itself. Errors along poorly measured directions are nearly invisible to it, so for an inverse problem like deblurring it tends to regularize too little. - It works for any method that can be fit to part of the data and used to predict the rest, and it can compare different methods as well as different values of $\lambda$. - The fold errors also say how uncertain $\mathrm{CV}(\lambda)$ is; the one-standard-error rule uses this to pick the most regularized $\lambda$ that is statistically tied with the minimum, which is more reliable from one data set to the next. - At a good choice of $\lambda$, the cross-validated residuals look like pure noise; structure left in them is bias, and residuals larger than the noise are variance. - Keep everything that depends on the data, including standardization, inside the folds, and keep a separate test set if you want an honest estimate of the final error.