Skip to main content
Open In Colab Now that we have introduced somewhat more formally the learning problem and its notation lets us study a simple but instructive regression problem that is known in the statistics literature as shrinkage. Suppose that we are given the training set X={x(1),…,x(m)}\mathbb{X} = \{x^{(1)}, \ldots, x^{(m)}\} together with their labels {y(1),…,y(m)}\{y^{(1)}, \ldots, y^{(m)}\}. We need to construct a model such that a suitably chosen loss function is minimized for a different set of input data, the so-called test set. The ability to correctly predict when observing the test set, is called generalization. Since the output yy is a continuous variable, this supervised problem is a regression problem (otherwise it is a classification problem). To make it concrete, read the input xx as normalized time over a single tidal cycle (x=0x = 0 at one high tide, x=1x = 1 at the next) and the output yy as the water level relative to mean sea level, in meters. A tide gauge records the level at a sequence of times, and those readings are your training points. The water climbs to a high, falls to a low, and rises again: one smooth cycle. We draw the readings from sin⁡(2πx)+ϵ\sin(2 \pi x) + \epsilon, with measurement noise ϵ∼N(μ=0, σ=0.25)\epsilon \sim N(\mu=0,\, \sigma=0.25). That underlying tidal curve is completely unknown to you; you see only the noisy gauge readings, and your hypothesis has to recover the smooth tide beneath them. 
def create_toy_data(func, sample_size, std, domain=[0, 1]):
    x = np.linspace(domain[0], domain[1], sample_size)
    np.random.shuffle(x)
    y = func(x) + np.random.normal(scale=std, size=x.shape)
    return x, y


def sinusoidal(x):
    return np.sin(2 * np.pi * x)


x_train, y_train = create_toy_data(sinusoidal, 10, 0.25)
x_test = np.linspace(0, 1, 100)
y_test = sinusoidal(x_test)

plt.figure(figsize=[10, 8])
plt.scatter(
    x_train, y_train, facecolor="none", edgecolor="b", s=50, label="training data"
)
plt.plot(x_test, y_test, "-g", label="true tide (unknown)")
plt.xlabel("$x$ (normalized tidal time)")
plt.ylabel("$y$ (water level, m)")
plt.legend()
plt.show()
Output from cell 3 Let us now pick the hypothesis set that corresponds in general to the following sum, g(x; θ)=θ0+∑j=1M−1θj ϕj(x)=θ⊤ϕ(x)g(x;\, \boldsymbol{\theta}) = \theta_0 + \sum_{j=1}^{M-1} \theta_j\, \phi_j(x) = \boldsymbol{\theta}^\top \boldsymbol{\phi}(x) where ϕj(x)\phi_j(x) are known as basis functions. Evaluating the basis at every training input and stacking the results row by row gives the design (featurization) matrix Φ\boldsymbol{\Phi}, whose ii-th row is ϕ(x(i))⊤\boldsymbol{\phi}(x^{(i)})^\top. For the polynomial basis ϕj(x)=xj\phi_j(x) = x^j it is the Vandermonde matrix Φ=[ϕ(x(1))⊤ϕ(x(2))⊤⋮ϕ(x(m))⊤]=[1x(1)(x(1))2⋯(x(1))M1x(2)(x(2))2⋯(x(2))M⋮⋮⋮⋱⋮1x(m)(x(m))2⋯(x(m))M]∈Rm×(M+1),\boldsymbol{\Phi} = \begin{bmatrix} \boldsymbol{\phi}(x^{(1)})^\top \\ \boldsymbol{\phi}(x^{(2)})^\top \\ \vdots \\ \boldsymbol{\phi}(x^{(m)})^\top \end{bmatrix} = \begin{bmatrix} 1 & x^{(1)} & (x^{(1)})^2 & \cdots & (x^{(1)})^M \\ 1 & x^{(2)} & (x^{(2)})^2 & \cdots & (x^{(2)})^M \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x^{(m)} & (x^{(m)})^2 & \cdots & (x^{(m)})^M \end{bmatrix} \in \mathbb{R}^{m \times (M+1)}, so that all mm predictions collapse into a single matrix product y^=Φθ\hat{\boldsymbol{y}} = \boldsymbol{\Phi}\boldsymbol{\theta}. With this notation the ridge objective is ∥Φθ−y∥2+λ∥θ∥2\lVert \boldsymbol{\Phi}\boldsymbol{\theta} - \boldsymbol{y} \rVert^2 + \lambda \lVert \boldsymbol{\theta} \rVert^2, whose minimizer solves the normal equations (Φ⊤Φ+λI) θ=Φ⊤y(\boldsymbol{\Phi}^\top \boldsymbol{\Phi} + \lambda \boldsymbol{I})\,\boldsymbol{\theta} = \boldsymbol{\Phi}^\top \boldsymbol{y}. The fitting code below builds exactly this Φ\boldsymbol{\Phi} with np.vander, then standardizes its columns for numerical stability. A set of Polynomial, Gaussian and Sigmoidal basis functions are plotted below. 
x = np.linspace(-1, 1, 100)
X_polynomial = PolynomialFeature(11).transform(x[:, None])
X_gaussian = GaussianFeature(np.linspace(-1, 1, 11), 0.1).transform(x)
X_sigmoidal = SigmoidalFeature(np.linspace(-1, 1, 11), 10).transform(x)

plt.figure(figsize=(20, 5))
for i, X in enumerate([X_polynomial, X_gaussian, X_sigmoidal]):
    plt.subplot(1, 3, i + 1)
    for j in range(12):
        plt.plot(x, X[:, j])
txt = "Polynomial, Gaussian and Sigmoidal basis functions"
plt.figtext(0.5, 0.01, txt, wrap=True, horizontalalignment="center", fontsize=12)

Text(0.5, 0.01, 'Polynomial, Gaussian and Sigmoidal basis functions')
Output from cell 4 Using the polynomial basis, we fit candidate hypotheses g(x; θ)g(x;\, \boldsymbol{\theta}) of increasing order MM to the training data. The figure below overlays the fits for M=0,1,2,3,9M = 0, 1, 2, 3, 9 on a single plot. Low-order models underfit, while M=9M=9 wiggles through every training point. 
plt.figure(figsize=[10, 8])
plt.scatter(x_train, y_train, facecolor="none", edgecolor="b", s=50, label="training data")
for M in [0, 1, 2, 3, 9]:
    coeffs = np.polyfit(x_train, y_train, M)
    plt.plot(x_test, np.polyval(coeffs, x_test), label=f"$M={M}$")
plt.ylim(-2.5, 2.5)
plt.xlabel("$x$ (normalized tidal time)")
plt.ylabel("$y$ (water level, m)")
plt.legend()
plt.show()
Output from cell 5 Our job is to find θ\boldsymbol{\theta} such that the polynomial expansion above fits the data we are given - as we will see there are multiple hypothesis that can satisfy this requirement. Consistent with the block diagram we need to define a metric, an figure of merit or loss function in fact, that is also a common metric in regression problems of this nature. This is the Mean Squared Error (MSE) function, which here plays the role of the empirical risk on the training set: Remp(θ)=1m∑i=1m(g(x(i); θ)−y(i))2R_{\text{emp}}(\boldsymbol{\theta}) = \frac{1}{m} \sum_{i=1}^m \big(g(x^{(i)};\, \boldsymbol{\theta}) - y^{(i)}\big)^2 We can visualize this loss for a specific hypothesis. The green segments below are the displacements (residuals) between each training point and the order-M=2M=2 hypothesis gg; the MSE is the average of their squares. 
M = 2
coeffs = np.polyfit(x_train, y_train, M)
y_hat = np.polyval(coeffs, x_train)

plt.figure(figsize=[10, 8])
plt.plot(x_test, np.polyval(coeffs, x_test), "-r", label=f"order-{M} hypothesis $g$")
plt.scatter(x_train, y_train, facecolor="none", edgecolor="b", s=50, label="training data")
for xi, yi, yh in zip(x_train, y_train, y_hat):
    plt.plot([xi, xi], [yi, yh], "-g")
plt.xlabel("$x$ (normalized tidal time)")
plt.ylabel("$y$ (water level, m)")
plt.legend()
plt.show()
Output from cell 6 The loss function chosen for this regression problem corresponds to the sum of the squares of the displacements of each data point and our hypothesis. The sum of squares in the case of Gaussian errors gives rise to an (unbiased) Maximum Likelihood estimate of the model parameters. Contrast this to the sum of absolute differences. Now our job has become to choose two things: the parameter vector θ∗\boldsymbol{\theta}^\ast and MM the order of the polynomial. Both define our hypothesis. If you think about it, the order MM defines the model complexity in the sense that the larger MM becomes the more the number of parameters we need to estimate and store. Obviously this is a trivial example and storage is not a concern here but treat this example as instructive for that it applies in far for complicated settings. Sweeping the model order MM and recording the MSE on the training set and on a held-out test set exposes the tension between fitting and generalization: the training error keeps decreasing with MM, but the test error eventually rises, the signature of overfitting. 
# held-out noisy test set drawn from the same distribution
x_eval = np.linspace(0, 1, 100)
y_eval = np.sin(2 * np.pi * x_eval) + np.random.RandomState(0).normal(scale=0.25, size=x_eval.size)

orders = list(range(0, 10))
train_mse = [np.mean((np.polyval(np.polyfit(x_train, y_train, M), x_train) - y_train) ** 2) for M in orders]
test_mse = [np.mean((np.polyval(np.polyfit(x_train, y_train, M), x_eval) - y_eval) ** 2) for M in orders]

plt.figure(figsize=[10, 8])
plt.plot(orders, train_mse, "-o", label="training")
plt.plot(orders, test_mse, "-o", label="test")
plt.xlabel("$M$ (polynomial order)")
plt.ylabel("MSE")
plt.legend()
plt.show()
Output from cell 7 Obviously you can reduce the training error to almost zero by selecting a model that is complicated enough (M=9) to perfectly fit the training data (if m is small). But this is not what you want to do. Because when met with test data, the model will perform far worse than a less complicated model that is closer to the true model (e.g. M=3). This is a central observation in statistical learning called overfitting. In addition, you may not have the time to iterate over M (very important in online learning settings). 
M = 9
coeffs = np.polyfit(x_train, y_train, M)[::-1]  # theta_0 .. theta_M

plt.figure(figsize=[10, 8])
plt.plot(range(M + 1), coeffs, "-o")
plt.xlabel("coefficient index $j$")
plt.ylabel(r"$\theta_j$")
plt.title("Unregularized $M=9$ coefficients")
plt.show()
Output from cell 8 To avoid overfitting we have multiple strategies. One straightforward one is evident by observing the wild oscillations of the θ\boldsymbol{\theta} elements as the model complexity increases. We can penalize such oscillations by introducing the ℓ2\ell_2 norm of θ\boldsymbol{\theta} in our loss function. Remp(θ)=1m∑i=1m(g(x(i); θ)−y(i))2+λ∥θ∥2R_{\text{emp}}(\boldsymbol{\theta}) = \frac{1}{m} \sum_{i=1}^m \big(g(x^{(i)};\, \boldsymbol{\theta}) - y^{(i)}\big)^2 + \lambda \|\boldsymbol{\theta}\|^2 This type of solution is called regularization and because we effectively shrink the parameter dynamic range it is also called in statistics shrinkage or ridge regression. We have introduced a new parameter λ\lambda that regulates the relative importance of the penalty term as compared to the MSE. This parameter together with the polynomial order is what we call hyperparameters and we need to optimize them as both are needed for the determination of our final hypothesis gg. The graph below show the results of each search iteration on the λ\lambda hyperparameter. Sweeping the regularization strength λ\lambda for the M=9M=9 model and recording the MSE on the training set and on a held-out test set reproduces Bishop’s λ\lambda-optimization curve. At very small λ\lambda the model overfits (low training error, high test error); too large a λ\lambda underfits; the test error is minimized at an intermediate λ∗\lambda^\ast. 
# fit M=9 with ridge on standardized polynomial features (numerically stable),
# sweeping the penalty lambda and recording MSE on the training and a held-out test set.
def std_ridge_predict(M, lam):
    P = np.vander(x_train, M + 1, increasing=True)[:, 1:]
    mu, sd = P.mean(0), P.std(0) + 1e-12
    Z = (P - mu) / sd
    theta = np.linalg.solve(Z.T @ Z + lam * np.eye(M), Z.T @ (y_train - y_train.mean()))
    return lambda xq: ((np.vander(xq, M + 1, increasing=True)[:, 1:] - mu) / sd) @ theta + y_train.mean()

x_eval = np.linspace(0, 1, 100)
y_eval = np.sin(2 * np.pi * x_eval) + np.random.RandomState(0).normal(scale=0.25, size=x_eval.size)

lambdas = np.logspace(-7, 3, 30)
train_mse, test_mse = [], []
for lam in lambdas:
    g = std_ridge_predict(9, lam)
    train_mse.append(np.mean((g(x_train) - y_train) ** 2))
    test_mse.append(np.mean((g(x_eval) - y_eval) ** 2))

plt.figure(figsize=[10, 8])
plt.plot(np.log(lambdas), train_mse, "-o", label="training")
plt.plot(np.log(lambdas), test_mse, "-o", label="test")
plt.xlabel(r"$\ln \lambda$")
plt.ylabel("MSE")
plt.legend()
plt.show()
Output from cell 9 
lambda_star = lambdas[int(np.argmin(test_mse))]
g = std_ridge_predict(9, lambda_star)

plt.figure(figsize=[10, 8])
plt.scatter(x_train, y_train, facecolor="none", edgecolor="b", s=50, label="training data")
plt.plot(x_test, g(x_test), "-r", label=r"regularized hypothesis ($M=9$, $\lambda=\lambda^\ast$)")
plt.plot(x_test, np.sin(2 * np.pi * x_test), "-g", label="true tide (unknown)")
plt.xlabel("$x$ (normalized tidal time)")
plt.ylabel("$y$ (water level, m)")
plt.legend()
plt.show()
Output from cell 10 The table below lists the fitted coefficient coordinates for M=3,6,9M = 3, 6, 9, without regularization and with ridge regularization (λ=0.1\lambda = 0.1). Without regularization the polynomial coefficients grow explosively with the model order, reaching the order of 10610^{6} at M=9M = 9 (the wild oscillation plotted above). Ridge penalizes the whole coefficient vector (typically excluding the intercept), so it reduces the overall coefficient norm; in this polynomial example the poorly constrained high-order behavior is what is most visibly suppressed. Individual coordinates, especially with correlated or standardized polynomial features, should not be read as each moving monotonically toward zero. Ridge also does not set individual high-order terms exactly to zero; that kind of sparsity is associated with L1 / Lasso.
coordinateM=3M=3M=3M=3 (ridge)M=6M=6M=6M=6 (ridge)M=9M=9M=9M=9 (ridge)
θ0\theta_0-0.100-0.0430.013-0.0430.020-0.043
θ1\theta_19.282-0.0264.5660.523340.170.477
θ2\theta_2-27.567-1.224-0.845-1.150-7.71e+03-0.890
θ3\theta_318.4160.925-34.211-0.8876.86e+04-0.713
θ4\theta_442.758-0.213-3.18e+05-0.320
θ5\theta_5-14.2940.4698.50e+050.004
θ6\theta_62.1431.056-1.37e+060.217
θ7\theta_71.30e+060.332
θ8\theta_8-6.73e+050.372
θ9\theta_91.47e+050.358
From shrinkage to pruning in real models. Ridge only shrinks the coefficients toward zero; it keeps all of them. In a large neural network the analogous, more aggressive step is pruning: the weights that regularization has driven close to zero are set to exactly zero, so they can be dropped to compress and speed up the model (this is closer in spirit to the L1 / Lasso sparsity mentioned above than to ridge). This is a standard step in production: NVIDIA’s TensorRT Model Optimizer prunes, imposes 2:4 structured sparsity (half the weights in every group zeroed), and quantizes (INT8 / FP8), and the TensorRT runtime then executes the resulting sparse, low-precision network on the GPU’s Sparse Tensor Cores. So “drop the coefficients close to zero” is not only a textbook idea, it is a core step in deploying real models.
Key references: (, n.d.; , n.d.; , n.d.)

​
References

  • (n.d.). Pattern Recognition and Machine Learning.
  • (n.d.). PRML: PRML algorithms implemented in Python.
  • (n.d.). Elements of Statistical Learning.
Edit this page on GitHub or file an issue.