ERM for polynomial regression: model mismatch and noise

This notebook supports the exercises on linear and affine empirical risk minimization. It keeps the same observation model,

$$ y_i=f_0(x_i)+\sigma\varepsilon_i+\Delta\delta_i, \qquad \varepsilon_i\sim\mathcal N(0,1), \qquad \delta_i\sim\mathrm{Bernoulli}(p), $$

but allows the true signal $f_0$ and the fitted model to be polynomials of different degrees.

The main parameters are collected in short cells. The fitting and plotting code is supplied, so the work consists mainly of changing values, interpreting plots, and connecting them to the formulas in the handout.

import numpy as np
import matplotlib.pyplot as plt

plt.style.use("seaborn-v0_8-whitegrid")
np.set_printoptions(precision=3, suppress=True)

SEED = 12
rng = np.random.default_rng(SEED)

Tools used throughout

Coefficients are stored in increasing order: [b, w, q] represents

$$b+wx+qx^2.$$

The empirical risk is the quadratic risk from the handout,

$$\widehat R_n=\frac{1}{2n}\sum_{i=1}^n (y_i-f(x_i))^2.$$

You do not need to modify the helper functions.

def polynomial_matrix(x, degree, include_intercept=True):
    """Columns 1, x, ..., x**degree; optionally remove the constant column."""
    matrix = np.vander(np.asarray(x), N=degree + 1, increasing=True)
    return matrix if include_intercept else matrix[:, 1:]


def polynomial_values(x, coefficients):
    return polynomial_matrix(x, len(coefficients) - 1) @ np.asarray(coefficients)


def generate_data(x, true_coefficients, sigma, Delta, p, rng):
    gaussian_noise = rng.normal(size=len(x))
    contamination = rng.binomial(1, p, size=len(x))
    clean_signal = polynomial_values(x, true_coefficients)
    y = clean_signal + sigma * gaussian_noise + Delta * contamination
    return y, clean_signal, gaussian_noise, contamination


def fit_polynomial(x, y, degree, include_intercept=True):
    X = polynomial_matrix(x, degree, include_intercept)
    coefficients = np.linalg.lstsq(X, y, rcond=None)[0]
    return coefficients


def predict_polynomial(x, coefficients, include_intercept=True):
    if include_intercept:
        return polynomial_matrix(x, len(coefficients) - 1) @ coefficients
    return polynomial_matrix(x, len(coefficients), False) @ coefficients


def empirical_risk(y, prediction):
    return 0.5 * np.mean((y - prediction) ** 2)


def describe_fit(name, y, prediction, coefficients):
    print(f"{name:24s} risk={empirical_risk(y, prediction):.4f}  coefficients={coefficients}")

1. The two noise sources

The four panels below use the same inputs and random draws. This isolates the effects of:

• Gaussian noise, controlled by $\sigma$;
• Bernoulli contamination, controlled by its size $\Delta$ and probability $p$.

The default true model is quadratic, while later sections fit models of other degrees.

# Parameters to change
n = 60
true_coefficients = np.array([1.0, -1.5, 0.8])  # 1 - 1.5x + 0.8x²
sigma = 0.45
Delta = 5.0
p = 0.10

def plot_noise_source_effect(n, true_coefficients, sigma, Delta, p):
    x = np.sort(rng.uniform(-2.0, 2.0, n))
    x_plot = np.linspace(-2.2, 2.2, 500)

    # Draw the random variables once so all panels are directly comparable.
    eps = rng.normal(size=n)
    delta = rng.binomial(1, p, size=n)
    clean = polynomial_values(x, true_coefficients)

    datasets = {
        "no noise": clean,
        "Gaussian only": clean + sigma * eps,
        "contamination only": clean + Delta * delta,
        "both sources": clean + sigma * eps + Delta * delta,
    }

    fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True, sharey=True)
    for ax, (name, values) in zip(axes.ravel(), datasets.items()):
        ax.scatter(x[delta == 0], values[delta == 0], s=24, label="ordinary point")
        if np.any(delta == 1):
            ax.scatter(x[delta == 1], values[delta == 1], s=45, marker="x",
                    color="C3", label=r"$\delta_i=1$")
        ax.plot(x_plot, polynomial_values(x_plot, true_coefficients),
                color="black", lw=2, label="true signal")
        ax.set_title(name)
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        ax.legend(fontsize=8)
    plt.tight_layout()
    plt.show()

    print(f"Observed contamination proportion: {delta.mean():.3f} (target p={p})")

plot_noise_source_effect(n, true_coefficients, sigma, Delta, p)

Questions

1. Set sigma = 0, then set Delta = 0. Which source moves every observation, and which moves only a random subset?

# Change variables here and use plot_noise_source_effect to visualize the effect of different noise sources
plot_noise_source_effect(n, true_coefficients, sigma, Delta, p)

2. Increase $p$ while decreasing $\Delta$. Can many moderate contaminations resemble a change in intercept?

# Change variables here and use plot_noise_source_effect to visualize the effect of different noise sources
plot_noise_source_effect(n, true_coefficients, sigma, Delta, p)

2. Linear through the origin versus affine regression

This section directly supports Exercises 1 and 2. When the true intercept is nonzero, the model $f_w(x)=wx$ is misspecified. An affine model can represent the signal exactly when the true signal is affine.

# Change b0 and rerun. Noise is initially disabled to isolate model mismatch.
w0 = 1.5
b0 = 2.0
xmin = -2.0
xmax = 2.0

def compare_linear_vs_affine_fit(w0, b0, xmin, xmax):
    x_linear = np.sort(rng.uniform(xmin, xmax, 35))
    y_linear = w0 * x_linear + b0

    w_origin = fit_polynomial(x_linear, y_linear, degree=1, include_intercept=False)
    ab_affine = fit_polynomial(x_linear, y_linear, degree=1, include_intercept=True)

    grid = np.linspace(xmin - 0.2, xmax + 0.2, 400)
    pred_origin = predict_polynomial(grid, w_origin, include_intercept=False)
    pred_affine = predict_polynomial(grid, ab_affine, include_intercept=True)

    fig, ax = plt.subplots(figsize=(8, 5))
    ax.scatter(x_linear, y_linear, color="black", s=28, label="data")
    ax.plot(grid, w0 * grid + b0, "k--", lw=2, label="true affine signal")
    ax.plot(grid, pred_origin, lw=2, label="fit through origin")
    ax.plot(grid, pred_affine, lw=2, label="affine fit")
    ax.set(xlabel="x", ylabel="y", title="Effect of omitting the intercept")
    ax.legend()
    plt.show()

    describe_fit("through origin", y_linear,
                predict_polynomial(x_linear, w_origin, False), w_origin)
    describe_fit("affine", y_linear,
                predict_polynomial(x_linear, ab_affine, True), ab_affine)

compare_linear_vs_affine_fit(w0, b0, xmin, xmax)

Questions

1. Try $b_0\in\{0,1,3\}$. Which estimator changes even though the true slope stays fixed?

# Change variables here and use compare_linear_vs_affine_fit
compare_linear_vs_affine_fit(w0, b0, xmin, xmax)

2. Shift every input to the right by increasing xmin and xmax. How does it affect the linear estimator?

# Change variables here and use compare_linear_vs_affine_fit
compare_linear_vs_affine_fit(w0, b0, xmin, xmax)

3. True polynomial degree versus fitted degree

Now choose the degree of the true polynomial and the degree used by ERM independently. The same noisy observations are used for every fitted degree, so differences come only from the model class.

• fit_degree < true_degree: underfitting or approximation error;
• fit_degree = true_degree: correctly specified model;
• fit_degree > true_degree: extra flexibility and possible sensitivity to noise.

from turtle import color


def experiment_polynomial_fit(true_coefficients, fit_degrees, n, sigma, Delta, p):
    experiment_rng = np.random.default_rng(40)
    x_train = np.sort(experiment_rng.uniform(-1, 1, n))
    y_train, clean_train, eps_train, delta_train = generate_data(
        x_train, true_coefficients, sigma, Delta, p, experiment_rng
    )

    x_test = np.linspace(-1.1, 1.1, 600)
    y_test_clean = polynomial_values(x_test, true_coefficients)

    fig, axes = plt.subplots(1, len(fit_degrees), figsize=(15, 4.3), sharex=True, sharey=True)
    for ax, degree in zip(axes, fit_degrees):
        coefficients = fit_polynomial(x_train, y_train, degree, True)
        train_prediction = predict_polynomial(x_train, coefficients, True)
        test_prediction = predict_polynomial(x_test, coefficients, True)
        train_risk = empirical_risk(y_train, train_prediction)
        clean_test_mse = np.mean((y_test_clean - test_prediction) ** 2)

        ax.scatter(x_train[delta_train == 0], y_train[delta_train == 0], s=20)
        ax.scatter(x_train[delta_train == 1], y_train[delta_train == 1],
                marker="x", s=45, color="C3")
        ax.plot(x_test, y_test_clean, "k--", lw=2, label="true signal")
        ax.plot(x_test, test_prediction, color="C1", lw=2, label="ERM fit")
        ax.set_title(f"fit degree {degree}\ntrain risk={train_risk:.3f}, clean MSE={clean_test_mse:.3f}")
        ax.set_xlabel("x")
        ax.legend(fontsize=8)
    axes[0].set_ylabel("y")
    axes[0].set_ybound(y_train.min() - 1, y_train.max() + 1)
    plt.tight_layout()
    plt.show()

Questions

1. Use the previous function for displaying the best fit polynomial for this experiment, using degrees up to 30. What do you observe as the degree gets higher?

# Put the parameters for the experiment here and run it
true_coefficients = np.array([0.8, -1.2, 0.0, 5.2])  # true degree 3
n = 55
sigma = 0.40
Delta = 4.0
p = 0.08

fit_degrees = [1, 3, 8, 30]

experiment_polynomial_fit(true_coefficients, fit_degrees, n, sigma, Delta, p)

2. Set sigma = Delta = 0. Can a fitted degree smaller than the true degree attain zero empirical risk? Can a larger degree attain zero empirical risk?

# Put the parameters for the experiment here and run it

experiment_polynomial_fit(true_coefficients, fit_degrees, n, sigma, Delta, p)

3. Restore the noise and try the same experiments with randomly generated true coefficients for higher degrees.

# Put the parameters for the experiment here and run it
true_degree = 15
true_coefficients = np.random.uniform(-2, 2, true_degree + 1)

experiment_polynomial_fit(true_coefficients, fit_degrees, n, sigma, Delta, p)

4. Increase the sample size while keeping the degrees fixed. Which fitted model benefits most?

# Put the parameters for the experiment here and run it
true_degree = 15
true_coefficients = np.random.uniform(-2, 2, true_degree + 1)
n = 500

fit_degrees = [1, 3, 8]
experiment_polynomial_fit(true_coefficients, fit_degrees, n, sigma, Delta, p)