Demo: Double Descent

We fit a nonlinear function using increasingly many nonlinear features. The fitted weights still come from least squares. We will compare training and test error as the number of weights passes the number of training examples.

import numpy as np
import matplotlib.pyplot as plt

plt.style.use("../figure_style.mplstyle")
n, d_max, n_test = 100, 400, 2000
input_dim, frequency_std, noise_std = 5, 1.0, 0.05
seeds = list(range(100, 116))
n_trials = len(seeds)
feature_counts = np.array([1, 5, 10, 15, 19, 20, 21, 25, 30, 50, 70, 85,
                           90, 95, 99, 100, 101, 105, 110, 130, 160, 200, 300, 400])

Fit a Nonlinear Function

Draw five inputs independently and uniformly from \([-1,1]\). The true function is \[ f(\mathbf x)=\sin(3x_1)+\sin(3x_2)+x_3^2. \] The observed label is \(y=f(\mathbf x)+\varepsilon\), with independent normal noise of standard deviation \(0.05\).

We use features \(\phi_j(\mathbf x)=\sqrt{2}\cos(\boldsymbol\omega_j^\top\mathbf x+b_j)\), where \(\boldsymbol\omega_j\sim N(\mathbf0,\mathbf I)\) and \(b_j\) is uniform on \([0,2\pi]\). Each model uses the first \(d\) columns of the same feature matrix, so increasing \(d\) adds features while preserving the existing ones. The model is nonlinear in the inputs and linear in the weights.

def true_function(x):
    return np.sin(3*x[:, 0]) + np.sin(3*x[:, 1]) + x[:, 2]**2


def run_curve(seed, label_noise_std=noise_std):
    rng = np.random.default_rng(seed)
    x = rng.uniform(-1, 1, (n, input_dim))
    x_test = rng.uniform(-1, 1, (n_test, input_dim))
    omega = rng.normal(0, frequency_std, (d_max, input_dim))
    phase = rng.uniform(0, 2*np.pi, d_max)
    X_full = np.sqrt(2)*np.cos(x @ omega.T + phase)
    X_test_full = np.sqrt(2)*np.cos(x_test @ omega.T + phase)
    y = true_function(x) + rng.normal(0, label_noise_std, n)
    y_test = true_function(x_test) + rng.normal(0, label_noise_std, n_test)
    results = {"train": [], "test": [], "smallest_singular_value": [],
               "noise_weight_change": []}
    for d in feature_counts:
        X, X_test = X_full[:, :d], X_test_full[:, :d]
        U, singular_values, Vt = np.linalg.svd(X, full_matrices=False)
        w = Vt.T @ ((U.T @ y) / singular_values)
        results["train"].append(float(np.mean((X @ w - y)**2)))
        results["test"].append(float(np.mean((X_test @ w - y_test)**2)))
        results["smallest_singular_value"].append(float(singular_values[-1]))
        # Expected squared difference between fits to noisy and clean labels.
        results["noise_weight_change"].append(float(label_noise_std**2 * np.sum(1 / singular_values**2)))
    return results

Where Does Test Error Peak?

Predict what happens when \(d\) crosses \(n=100\). Can a model with more than 100 weights predict better than every model with fewer than 100 weights?

We repeat the experiment 16 times with independent inputs, features, and noise. Lines show median errors, and shading spans the 25th to 75th percentiles across repetitions. These bands show variation across repetitions. The vertical axis is logarithmic above \(10^{-4}\) and linear below it to include zero.

trials = [run_curve(seed) for seed in seeds]
arrays = {key: np.array([trial[key] for trial in trials]) for key in trials[0]}
medians = {key: np.median(value, axis=0) for key, value in arrays.items()}
bands = {key: np.quantile(value, [.25, .75], axis=0) for key, value in arrays.items()}

def plot_errors(widths, medians, n, d_max, mobile=False, bands=None):
    """Median errors with interquartile bands and exact zeros on one axis."""
    fig, ax = plt.subplots(figsize=(4.1, 3.6) if mobile else (7.2, 3.8), layout="constrained")
    maximum = 0.0
    for key, color, style, label in [("train", "#44AA99", "-", "Training error"),
                                      ("test", "#332288", ":", "Test error")]:
        values = np.asarray(medians[key]).copy()
        if key == "train":
            values[widths >= n] = 0.0
        maximum = max(maximum, float(values.max()))
        if bands is not None:
            lower, upper = np.asarray(bands[key]).copy()
            if key == "train":
                lower[widths >= n] = upper[widths >= n] = 0.0
            ax.fill_between(widths, lower, upper, color=color, alpha=.18, linewidth=0)
            maximum = max(maximum, float(upper.max()))
        ax.plot(widths, values, color=color, linewidth=2.2, linestyle=style,
                dash_capstyle="round", label=label)
    ax.set_yscale("symlog", linthresh=1e-4, linscale=.5)
    ax.set(xlim=(0, d_max), ylim=(-1e-5, maximum * 1.8),
           xticks=[0, n, 200, 300, 400], xlabel="$d$", ylabel="MSE")
    powers = np.arange(-4, int(np.floor(np.log10(maximum * 1.8))) + 1)
    ax.set_yticks([0, *10.0**powers], labels=["0", *[f"$10^{{{k}}}$" for k in powers]])
    ax.legend(loc="upper right")
    return fig


fig = plot_errors(feature_counts, medians, n, d_max, bands=bands)
plt.show()

What happened? Training error reaches zero once there are enough independent features to fit every label. Test error rises near this point, then falls as more nonlinear features become available. Compare the test error at the largest width with its best value below 100 weights.

Check the Small Singular Values

At which width is the smallest singular value smallest? What does that imply for the change in the weights caused by noise?

fig, ax = plt.subplots(figsize=(7, 3), layout="constrained")
ax.plot(feature_counts, medians["smallest_singular_value"])
ax.set(xlabel="$d$", ylabel="Smallest singular value", yscale="log")
plt.show()

What happened? Small singular values near the interpolation threshold make the fitted weights sensitive to label noise. Beyond the threshold, adding columns increases the nonzero singular values and gives more choices of weights with the same training predictions.

Change the Noise Level

Keep the inputs and features fixed within each repetition, and increase the noise standard deviation from 0.05 to 0.15. Predict how this changes the peak and the errors away from it.

noise_std_new = 0.15
trials_new = [run_curve(seed, noise_std_new) for seed in seeds]
test_new = np.median([trial["test"] for trial in trials_new], axis=0)

fig, ax = plt.subplots(figsize=(7, 3.2), layout="constrained")
ax.plot(feature_counts, medians["test"], label="Noise SD 0.05")
ax.plot(feature_counts, test_new, linestyle=":", label="Noise SD 0.15")
ax.set(xlabel="$d$", ylabel="Median test MSE", yscale="log")
ax.legend()
plt.show()

What happened? More label noise raises the errors. The interpolation threshold stays at \(d=n\), because the number of examples and the ranks of the feature matrices stay the same. The size of the second descent depends on both the target function and the noise level.