Overfitting, the Bias-Variance tradeoff, Regularization, and Double Descent

We have seen a number of supervised learning algorithms, which learn deterministic mappings from input to output. In this chapter, we will analyze overfitting and underfitting in supervised learning, in order to understand the general conditions under which we can expect a supervised learning algorithm to perform well.

Preamble: Run the cells below to import the necessary Python packages

Open this notebook in Google Colab:

Polynomial Regression¶

Recall that, in supervised learning, we fit a linear model to a dataset $X \in \mathbb{R}^{N_\text{data} \times N_\text{features}}$ and a vector of labels $y \in \mathbb{R}^{N_\text{data}}$. Here, we will consider the simple one-dimensional case $X \in \mathbb{R}^{N_\text{data} \times 1}$ and $y \in \mathbb{R}^{N_\text{data}}$

One option to find the relationship between the features and the labels is to fit a linear model of the form

where $\boldsymbol{\theta} \in \mathbb{R}^{1}$ is a scalar, because each data point has only one feature, and we want to map it to one regression target value.

Bias as a feature in linear regression¶

You might recall that univariate linear regression models can be written as

where $\mathbf{X}, \mathbf{y} \in \mathbb{R}^{N_\text{data}}$ and all model weights are scalars. We can “lift” the dimensionality of the feature space by adding a constant feature to the data, so that the model becomes

where $\boldsymbol{\theta} \in \mathbb{R}^{2}$ and $\mathbf{X} \in \mathbb{R}^{N_\text{data} \times 2}$, with the first column of $\mathbf{X}$ being all ones.

Nonlinear features are more expressive¶

What if we wanted to fit a polynomial instead? We could instead solve a polynomial regression problem of the form

where $\Phi(\mathbf{X})$ is a matrix of features that are polynomial transformations of the original features,

Thus, in this lifted feature space the labels remain the same $\mathbf{y} \in \mathbb{R}^{N_\text{data}}$, but the transformed features $\Phi(\mathbf{X}) \in \mathbb{R}^{N_\text{data} \times N_\text{features}}$ where $N_\text{features} > 1$. As a result, the model weights $\boldsymbol{\theta} \in \mathbb{R}^{N_\text{features}}$.

Geature-wise transformations are sometimes called basis functions or kernels. For a single datapoint, the resulting regression problem has the form

that is linear in $\theta$ but non-linear in $x$ because the features

are non-linear. Using these features, we can fit any polynomial of degree $p$ by first transforming the features and then fitting a linear model using standard methods

Polynomial Regression¶

Instead of directly performing nonlinear transformations of the input features, we can use the PolynomialFeatures class from sklearn.preprocessing to compute the lifted feature matrix for us. Because it’s built into scikit-learn, it automatically has the right API to be used within a standard machine learning model. After transforming the input data, we will pass the transformed data to scikit-learn’s LinearRegression or Ridge classes to fit a linear model to the transformed data

scitkit-learn allows us to combine these operations into a single combined model with the Pipeline class. We will score the quality of the model using the $R^2$ score on both the training and test data

Why not always lift the features? Diminishing returns of the Vandermonde matrix¶

The Vandermonde matrix $\Phi \in \mathbb{R}^{N \times p}$ is a matrix whose rows are the features $\phi(x)$ for a set of points $x_1,\ldots,x_n$. For example, if we have $n=3$ points $x_1=1$, $x_2=2$, $x_3=3$, and $p=2$, then the Vandermonde matrix is $$

$$

Below, we plot the Vandermonde matrix for $N_\text{data}=3$ points and $p=0,2,4,9$

Solving yet another least-squares problem¶

In principle, we can simply fit a polynomial of degree $p$ to a set of points $(x_1,y_1),\ldots,(x_n,y_n)$ by solving the linear system $\Phi\theta = y$.

$$

=

$$

Notice how we “linearized” our problem by applying nonlinear functions to our features. This is the basis of many “kernel methods” in machine learning.

Properties of the Vandermonde matrix¶

• The Vandermonde matrix $\Phi$ is a square matrix of size $n\times p$, where $n$ is the number of datapoints and $p$ is one plus the degree of the polynomial that we want to fit. In this sense, the “features” of the model are just higher powers of existing features, making the features redundant.

• If $p = n$, then the Vandermonde matrix is square and invertible, and the linear system has a unique solution via $\Phi^{-1}$. Notice how the Vandermond matrix always has full rank in this case, if the points $x_1,\ldots,x_n$ are all distinct. However, the matrix becomes singular if any of the datapoints are repeated.

• However, if $p > n$, then the Vandermonde matrix is rectangular and the linear system is not guaranteed to have a unique solution. Conversely, if $p < n$, then the Vandermonde matrix is rank-deficient and the linear system is not guaranteed to have a solution. In this case, we can find the least-squares solution using the Moore-Penrose pseudoinverse $(\Phi^\top\Phi)^{-1}\Phi^\top$.

• Important, as the number of points increases, the Vandermonde matrix becomes ill-conditioned, which degrades the accuracy of the solution.

How can we interpret this effect?¶

Let’s try defining a data matrix with $N_\text{data}=2$. What do the column vectors of a high-degree Vandermonde matrix look like?

The largest element of our data vector $x$ begins to dominate at higher powers of $x$, causing the column vectors to become more and more similar to each other, and resulting in the least-squares problem becoming ill-conditioned.

Nothing comes for free! If we don’t have sufficient features or measurement channels, we can’t always just make up new ones as functions of our existing features. Eventually, redundancy catches up with us. We can’t convert an underdetermined problem to an overdetermined problem by adding redundant features.

Residual Scatter: Determining model complexity without looking at the test dataset¶

• Although fitting a linear model does not work well, quadratic or cubic polynomials seem to improve the fit to both the training and test data. Increasing the model expressivity therefore seems to improve the model’s ability to fit the training data, and generalize to unseen test data.

• We can also assess whether a given model is underfitting by plotting the residuals, which are the differences between the predicted values and the true values. If the residuals are not random, then the model is underfitting.

  If the true values are $y$ and the predicted values are $\hat{y}$, then the residuals are $y - \hat{y}$. The reason we expect uniform scatter is that our linear regression model can be seen as defining a data-generating process of the form

  $$y = \theta^\top \phi(x) + \epsilon$$

  (12)

  where $\epsilon$ is a random variable with zero mean and constant variance. The residuals are therefore given by

  $$y - \theta^\top \phi(x) = \epsilon$$

  (13)

We can see one heuristic for determining whether we have chosen a sufficiently-expressive model: Uniform scatter in the residuals

Hypothetically, we want to view our data generating process as comprising a deterministic function $f$ and a random noise term $\epsilon$. We want to fit a model that is close to $f$ and that is able to capture everything except the variability in $\epsilon$.

where $f(x_i) = \theta_0 + \theta_1 x_i + ...$ and $\epsilon_i$ is a random variable with mean 0 and variance $\sigma^2$. When we compute our residuals, we are computing the difference between the true values $y_i$ and our model’s predictions $\hat{y}_i = f(x_i)$.

Thus we want to $r_i$ to look like a zero-centered random distribution, consistent with our implicit model of the generating process for this dataset.

What happens when our model complexity is too high?¶

• Increase the degree of our polynomial improved our fit accuracy by producing a model that explained more of the variance in the data.
• Our residuals appeared more uniform as well, suggesting that our model stopped underfitting

What happens if we further increase the degree of the polynomial?

Runge’s phenomenon: can fit an $N$-point dataset perfectly with a polynomial of degree $N$.

Let’s quantify this effect¶

Model complexity¶

• This is a nebulous concept in machine learning. In principle, it’s a function of the number of trainable parameters in a model

• However, we’ve seen before that regularizers and constraints can reduce the effective number of degrees of freedom, and that not all parameter combinations are “reachable” during training.

• Nonetheless, we’ll treat the number of parameters as an upper bound on the complexity of a model.

• Other measures of complexity include the VC dimension, the Rademacher complexity, and the margin of a classifier.

Plotting train/test error vs. model complexity¶

• Usually, model complexity is controlled by a hyperparameter, such as the degree of a polynomial, the number of layers in a neural network, the number of trees in a random forest, or the strength of a regularizer.

• We can see how the training versus validation error changes as we vary the hyperparameter, to see how overfitting emerges.

• By convention, we normally plot error instead of accuracy, so that the error is always decreasing as the model improves. In this case, we will record $1 - R^2$ as our error metric.

Looking at this graph, we can see several patterns

• The training error decreases as the model complexity increases. This is because the model is able to fit the training data more closely as it becomes more complex.

• The test error initially decreases as the model complexity increases, as the model reaches a model class that is able to fit the data more closely.

• The test error then increases as the model complexity increases further. This is because the model is overfitting to the training data, and is not able to generalize to the test data.

We’d say that the best model is the one that minimizes the test error.¶

The Bias-Variance Tradeoff¶

We will return to our convention that each dataset comes from an underlying generator, which is unknown. We can think of the data as being generated by a function $f(x)$, plus some random noise $\epsilon$.

Suppose we fit a model $\hat{f}_\theta(x)$ to the data, where $\theta$ are the parameters of the model. We can write the overall mea squared error of the model as

where we are taking the expectation over the data points (indexed by $i$) in the dataset.

We can expand this term by adding and subtracting the unknown true generating function $f(x)$,

to get

Simplifying,

The last term is the cross-term between the noise and the model error. We can simplify this term by noting that the noise is independent of the model error, so that the expectation of the product is the product of the expectations. This term is therefore zero, and we have

We can now expand the second term in the sum by adding and subtracting the expected value of the model error $\mathbb{E}_i[\hat{f}_\theta(x_i)]$, which we denote as $\bar{f}$,

Substituting this back into the equation for the MSE, we get

The last cross term is again zero, because the model error is independent of the expected model error (there is no covariance between the two). We can simplify our expression to

The first term is the irreducible error, which is the error that we get no matter how good our model is. The second term is the bias of the model, which the error associated with making a naive prediction that all data points are equal to the expected value of the model. This term is minimized if the data has a very simple underlying structure (like a constant). The third term is the variance of the model across all datapoints.

[1] Hastie, T., Tibshirani, R., & Friedman, J. (2009). The elements of statistical learning: data mining, inference, and prediction. Springer Science & Business Media.

[2] Breiman, L. (2001). Statistical modeling: The two cultures. Statistical Science, 16(3), 199-231.

Overfitting¶

• A very expressive model (e.g., a high degree polynomial) fits the training dataset perfectly.

• But the model makes highly incorrect predictions outside this dataset, and doesn’t generalize.

• We would say that the overfit model exhibits high variance in its score across different hypothetical testing sets

• We could also say that the model has low bias because it will fit any training data very well.

Underfitting¶

• A small model (e.g. a straight line), will not fit the training data well.

• Therefore, it will also not be accurate on new data.

• The model will, however, exhibit low variance in its error on random testing sets.

• We say that it has high bias due to its strong tendency to do a limited number of things.

Determining Overfitting vs. Underfitting¶

We can diagnose overfitting and underfitting by measuring performance on a the held out test dataset (data that was not used for training).

• If training perforance is high but test performance is low, we are overfitting.

• If training and testing performance are both low, we are underfitting.

The gap between train and test scores is therefore one measure of overfitting. Larger gaps indicate more overfitting. How do we choose the correct model complexity? For our polynomial example, $p$ is a hyperparameter. We need to tune it on a validation set, or run cross-validation on our training partition

How to Fix Underfitting¶

• Feature engineering: find better features that will make the dataset easier to fit.

• Pick a more expressive model class (random forests or neural networks instead of linear models).

• Pick an optimization algorithm that can find better parameters (SGD instead of gradient descent).

How to Fix Overfitting¶

• Use a simpler model class (a linear model instead of a neural net)

• Use fewer trainable parameters (a smaller neural net)

• Keep the same model, but collect more training data

• Regularization: modify the training process to penalize overly complex models.

Regularization¶

• Regularization is a technique to prevent overfitting (or increase bias) by adding a penalty term to the loss function.

• We will try applying L2 or Ridge regularization to our polynomial regression model.

• We will also calculate the accuracy on the training set and the test set. We will again use the coefficient of determination $R^2$ as our metric.

• Why are we using L2 regularization instead of L1 (Lasso)?

L2 Regularization for Polynomial Regression¶

Let’s consider an application to the polynomial model we have seen so far. Given polynomial features $\phi(x)$, we optimize the following objective:

Regularization reduces the variance of the model¶

We implement regularized and polynomial regression of degree 15 on three random training sets sampled from the same distribution.

We can get a better idea of what’s going on by printing the learned weights $\theta$ for each model.

The bias-variance tradeoff¶

• I can either have a jack-of-all-trades model that is good at everything, but which requires a lot of data to train, or I can have a model that is good at a specific task, and which requires less data to train because it has high bias.

• Underfitting: high bias, low variance. The model has a strong “bias” towards a certain shape, and is inflexible to capture all of the peculiarities of the training data.

• Overfitting: low bias, high variance. The model is flexible, but overfits because it attempts to capture all of the peculiarities of the training data, some of which may result from noise or other factors that are not present in the test data.

• In physics, we often have “inductive bias” towards certain models. For example, we might know that energy is conserved, or that certain symmetries exist. This is a form of bias that can help us learn from less data, by restricting model space.

Hyperparameter Search¶

• The Ridge penalty parameter $\lambda$ as a hyperparameter, because it’s a high-level parameter that controls the properties of the model, rather than a parameter that is learned from the data.

• We can set $\lambda$ by assessing performance on a held-out validation dataset, or we can use cross-validation to select $\lambda$ from subsets of the training data.

L1 vs. L2 Regularization¶

Recall that L2 regularization (Ridge or Tikhonov) is defined as

L1 regularization (Lasso or Sparse) is defined as

Where $||\theta||_1 = \sum_{j=1}^d |\theta_j|$ is the L1 norm of $\theta$.

These two regularizers look pretty similar: both deter the values of $\theta$ from becoming too large. Mathematically, both are convex, both are non-negative, and both are differentiable almost everywhere.

However, in practice, they turn out to have very different effects on the learned model. To see this, we will re-run the polynomial regression example with both L1 and L2 regularization, and then see how the learned weights change as $\lambda$ is varied.

What happened here?¶

• We can see that L1 regularization has a tendency to set many of the weights to zero, while L2 regularization tends to smoothly shrink the weights towards zero.

• In practice, L2 regularization leads to a smoother model. in which all of the features are used, but none of them are used too much. These models often have smoother decision boundaries than L1-regularized models, and they are more robust to noise or sampling bias in the training data.

• However, L1 regularizion forces the model to avoid using many of the features, which can make the model more interpretable. For example, if two features contain redundant information, L1 regularization will help the model to choose one of them. This can be useful in feature selection or in reducing the complexity of a model.

• In modern deep learning, L2 regularization is used more often. However, in smaller-scale machine learning, when feature selection matters it is commong to first fit an L1-regularized model, and then drop all features from the training dataset that have zero weights (this is equivalent to dropping the corresponding columns of the feature matrix). Then an L2-regularized model is fit on this reduced dataset. This procedure is called *sequential feature selection.

An emerging topic: double descent¶

• We’ve seen how regularization helps us to control the complexity of a model, which is equivalent to navigating the bias-variance tradeoff.

• However, what happens if we increase the model complexity well beyond the point where we start to overfit?

• This example is based on Daniela Witten’s description of double descent

We will now fit a polynomial regression model, and visually inspect how the model behaves as we increase the degree of the polynomial.

We can see that the model slightly overfits the data, but at some point the overfitting stops getting any worse. Our intuition is that, once the model has enough capacity to fit all $N$ training datapoints, the extra degrees of freedom in the model are underspecified by the data.

What’s going on?¶

• Double descent was first described in 2018!

• It seemingly violates classical statistics, where too many parameters lead to overfitting and thus rising test error.

• When $p > n$ (more features than data), the model is overparameterized and the training error can be made arbitrarily small.

• Unexpectedly, the test error starts decreasing again as $p = n$ (the interpolation threshold, where the number of features is equal to the number of data points)

Header from Belkin et al 2019

Fig. 1 from Belkin et al 2019

Fig. 2 from Belkin et al 2019

How does this happen? The overparameterized regime¶

• The least-squares problem is underdetermined, and so there are multiple possible solutions once we pass the interpolation threshold.

• Given many equivalent solutions, scikit-learn’s Ridge chooses the solution that minimizes the L2 norm of the coefficients

• When $p$ gets arbitrarily large, then there are more possible solutions to choose from. The one that minimize the L2 norm of the coefficients ends up being closer to the correct solution, and the test error decreases again.

• Double descent is a loose argument for why neural networks work so well, despite often having many more parameters than data points. Overparametrization gives them more flexibility to find a good solution.

Implicit Regularization¶