CS 109A/STAT 121A/AC 209A/CSCI E-109A¶

Lab 5: Regularization¶

Harvard University
Fall 2017
Instructors: Pavlos Protopapas, Kevin Rader, Rahul Dave, Margo Levine

Table of Contents¶

0. Learning Goals
1. Introduction to regularized regression
2. Ridge regression with one predictor on a grid
3. Ridge regression with polynomial features on a grid
4. Cross-validation
5. Refitting on full training set

6. Feature selection with LASSO regression - good for homework 4!

Part 0: Learning Goals¶

In this lab we continue where we left off in Lab 4, with regularized regression. By the end of this lab, you will be able to:

• Implement ridge and LASSO regression using sklearn.
• Interpret the results of ridge and LASSO regression, and compare to the results from simple and multiple linear regression.

This lab maps on to lectures 7 and 8 and homework 4.

Part 1: Introduction to regularized regression¶

Recall from lecture the main idea of regularization. In the ordinary least squares problem we minimize the loss function

\begin{equation} L(\mathbf{\beta}) = \frac{1}{n} \sum_{i = 1}^n |y_i - \mathbf{\beta}^T \mathbf{x}_i|^2, \end{equation}

to determine regression coefficients $\mathbf{\beta}$. Here $y_i$ is the response variable for observation $i$, and $\mathbf{x}_i$ is a vector from the predictor matrix corresponding to observation $i$.

The general idea behind regularization is to penalize the loss function to account for possibly very large values of the coefficients $\mathbf \beta$. The aforementioned optimization problem is then adjusted accordingly. Instead of minimizing $L(\mathbf{\beta})$, we minimize the regularized loss function

\begin{equation} L_{\mathrm{reg}}(\mathbf{\beta}) = L(\mathbf{\beta}) + \lambda R(\mathbf{\beta}), \end{equation}

where $R(\mathbf{\beta})$ is a penalty function and $\lambda$ is a scalar that weighs the relative importance of this penalty. In this lab we will explore two regularized regression models, ridge and LASSO. In ridge regression, the penalty function is the sum of the squares of the parameters, giving the regularized loss function

\begin{equation} L_{\mathrm{Ridge}}(\mathbf{\beta}) = \frac{1}{n} \sum_{i = 1}^n |y_i - \mathbf{\beta}^T \mathbf{x}_i|^2 + \lambda \sum_{j = 1}^d \beta_j^2. \end{equation}

In LASSO regression the penalty function is the sum of the magnitudes of the parameters, leading to

\begin{equation} L_{\mathrm{LASSO}}(\mathbf{\beta}) = \frac{1}{n} \sum_{i = 1}^n |y_i - \mathbf{\beta}^T \mathbf{x}_i|^2 + \lambda \sum_{j = 1}^d |\beta_j|. \end{equation}

We will show how these optimization problems can be solved with sklearn to determine the model parameters $\mathbf \beta$. We will also show how to choose $\lambda$ appropriately via cross-validation.

Let's continue working with our data from last time. We load and split the data as in Lab 4.

Here x and y are the predictor and measured response variables, and f is the true response.

We split the sample data into training and testing sets.

Part 2: Ridge regression with one predictor on a grid¶

To begin, we'll use sklearn to do simple linear regression on the sampled training data. We'll then do ridge regression with the same data, setting the penalty parameter $\lambda$ to zero. Setting $\lambda = 0$ reduces the ridge problem to the simple ordinary least squares problem, so we expect the results of these models to be identical.

We will use the above $\beta$ coefficients as a benchmark for comparision to ridge and LASSO methods. Let's see that we get the same coefficients with ridge regression.

For reference, here is the ridge regression documentation. Notice that the weight $\lambda$ is referred to as alpha in the documentation.

The snippet of code below implements the ridge regression with $\lambda = 0$.

The beta coefficients for linear and ridge regressions coincide for $\lambda = 0$, as expected. We plot the data and fits.

EXERCISE: Play around with the values of $\lambda$ in the ridge regression code. Increase $\lambda$ from 0 to .01, from 0.01 to 1, from 1 to 5. What do you observe? What happens as $\lambda$ goes to $\infty$?

YOUR DISCUSSION HERE:

Part 3: Ridge regression with polynomial features on a grid¶

Now we'll make a more complex model by adding polynomial features. Instead of building the linear model $y = \beta_0 + \beta_1 x$, we build a polynomial model $y = \beta_0 + \beta_1 x + \beta_2 x^2 + \ldots \beta_d x^d$ for some $d$ to be determined (see Lab 4 for details on choosing hyper-parameter $d$). This regression will be linear though, since we'll be treating $x^2, \ldots, x^d$ themselves as predictors in the linear model. We did this in Lab 4 but it's worth a review.

We map $x$ to $1, x, x^2, \ldots, x^d$, and then build a linear regression model on this linear function of polynomial features. To do this, we use sklearn to build what is known as the Vandermonde matrix, the generalizaiton of the predictor matrix $X$ discussed in Lab 3. For example, if we have three observations

\begin{equation*}\begin{pmatrix} x_1 \\ x_2 \\ x_3\\ \end{pmatrix}, \end{equation*}

and we want polynomial features up to and including degree 4, we build the predictor matrix

\begin{equation*}\begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ \end{pmatrix} \rightarrow X = \begin{bmatrix} x_1^0 & x_1^1 & x_1^2 & x_1^3 & x_1^4\\ x_2^0 & x_2^1 & x_2^2 & x_2^3 & x_2^4\\ x_3^0 & x_3^1 & x_3^2 & x_3^3 & x_3^4\\ \end{bmatrix} = \begin{bmatrix} 1& x_1^1 & x_1^2 & x_1^3 & x_1^4\\ 1 & x_2^1 & x_2^2 & x_2^3 & x_2^4\\ 1 & x_3^1 & x_3^2 & x_3^3 & x_3^4\\ \end{bmatrix}. \end{equation*}

EXERCISE: Before we continue working with the data, make a toy vector called toy, where

\begin{equation} \mathrm{toy} = \begin{pmatrix} 0 \\ 2 \\ 5 \\ \end{pmatrix} . \end{equation} Build the feature matrix up to (and including) degree 4. Confirm that the entries in the matrix are what you'd expect based on the above discussion.

We now continue working with our data. We write a function to make polynomial features of given degrees as we did in Lab 4, and we store the features in a dictionary.

EXERCISE: Fill in the code below to perform ridge regression on the training data for the given set of $\lambda$. Then predict on the grid and store the results in ypredict_ridge.

As you can see, as we increase $\lambda$ from 0 to 1, we start out overfitting, then doing well, and then our fits develop a mind of their own irrespective of data, as the penalty term dominates.

EXERCISE: What would you expect if you compared a performance metric between these models on a grid?

YOUR DISCUSSION HERE:

Part 4: Cross-validation¶

Let's use cross-validation to determine the critical value of $\lambda$, which we'll refer to as $\lambda^*$. To do this we use the concept of a meta-estimator from scikit-learn. As the API paper from Lab 4 explains:

In scikit-learn, model selection is supported in two distinct meta-estimators, GridSearchCV and RandomizedSearchCV. They take as input an estimator (basic or composite), whose hyper-parameters must be optimized, and a set of hyperparameter settings to search through.

The concept of a meta-estimator allows us to wrap, for example, cross-validation, or methods that build and combine simpler models or schemes. For example:

est = Ridge()
parameters = {"alpha": [1e-8, 1e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 1e-2, 1e-1, 1.0]}
gridclassifier=GridSearchCV(est, param_grid=parameters, cv=4, scoring="neg_mean_squared_error")

The GridSearchCV replaces the manual iteration over the folds using KFolds and the averaging we did in Lab 4, doing it all for us. It takes a hyper-parameter grid in the shape of a dictionary as input, and sets $\lambda$ to the values you want to try, one by one. It then trains the model using cross-validation, and gets the error for each value of the hyper-parameter $\lambda$. Finally it compares the errors for the different $\lambda$'s, and picks the best choice model.

EXERCISE: Use the function above to fit the model on the training set with 4-fold cross validation. Save the fit as the variable fitmodel.

We also output the mean cross-validation error at different $\lambda$ (with a negative sign, as scikit-learn likes to maximize negative error which is equivalent to minimizing error).

EXERCISE: Plot log10-log10 plot of -fit_scores versus fit_lambdas.

Part 5: Refitting on full training set¶

We now refit the estimator on the training set, and calculate and plot the test set error and the polynomial coefficients. Notice how many of these coefficients have been pushed to lower values or 0.

EXERCISE: Assign to variable est the classifier obtained by fitting the entire training set using the best $\lambda$ found above. Assign the predictions to the variable ypredict_ridge_best.

END OF LAB¶

Part 6: Feature selection with LASSO regression¶

Below is a completely worked example of feature selection with LASSO, which will be helpful for homework 4. For reference here is the documentation for LASSO.