10. Exercise: Linear Regression

10.1. Import modules

# Common imports
import numpy as np
import os

# To plot pretty figures
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt

10.2. Generate data

# Let us generate noisy data with a linear feature

# to make this notebook's output stable across runs
np.random.seed(42)

# X are picked uniform random [0,2]
X = 2 * np.random.rand(100, 1)
# Linear relation to the predicted value, but with Gaussian noise (mean=0, variance=1)
y = 2 + 4 * X + np.random.randn(100, 1)

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(X, y, "b.")
ax.set_xlabel(r'$x_1$')
ax.set_ylabel(r'$y$');

10.3. Linear regression

10.3.1. Linear regression using scikit-learn

# Fit the linear regression model with sklearn
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)

# Print the linear regression parameters (with a specified precision)
with np.printoptions(precision=4):
    print(lin_reg.intercept_, lin_reg.coef_)

[2.2151] [[3.7701]]

X_predict = np.array([[0], [2]])
y_predict = lin_reg.predict(X_predict)

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(X_predict, y_predict, "r-",label="Prediction")
ax.plot(X, y, "b.",label="Data")
ax.set_xlabel('X')
ax.set_ylabel('y')
ax.legend(loc='best');

10.3.2. Linear regression using the Normal Equation

Let us create the design matrix X_d for the linear model. It is a linear model with two features corresponding to the two terms in a polynom of order 1:

• an \(X^0\) (bias)

• an \(X^1\) (linear feature)

m = len(X) # number of instances
# The design matrix will have shape (m x 2)
# It is created by adding X^0 = 1 to each instance
X_d = np.c_[np.ones((m, 1)), X]  

X_d.shape

(100, 2)

### EXERCISE
#
# (a) Solve the normal equation to compute the fit parameters
# (b) Print the best fit parameters
# (c) Plot the fitted linear model together with the data
#

10.3.3. Comment on the base for the sklearn LinearRegression method

The LinearRegression class is based on the scipy.linalg.lstsq() function (the name stands for “least squares”), which you could call directly:

theta_best_svd, residuals, rank, s = np.linalg.lstsq(X_d, y, rcond=1e-6)
theta_best_svd

array([[2.21509616],
       [3.77011339]])

This function computes \(\mathbf{X}^+\mathbf{y}\), where \(\mathbf{X}^{+}\) is the pseudoinverse of \(\mathbf{X}\) (specifically the Moore-Penrose inverse). You can use np.linalg.pinv() to compute the pseudoinverse directly:

np.linalg.pinv(X_d).dot(y)

array([[2.21509616],
       [3.77011339]])

10.3.4. Linear regression using batch gradient descent

Gradient descent optimization is not really needed for linear regression since we can solve the normal equation exactly. However, it is instructive to implement it yourself and to be able to compare results with a known solution.

• Define a function to perform gradient descent optimization on the batch of data \((\boldsymbol{X},\boldsymbol{y})\) that was produced above.

• For linear regression, the gradient of the cost function can be expressed as a matrix equation and is therefore easy to evaluate. See the chapter Linear Regression in the lecture notes.

• Perform up to 1000 iterations where each iteration is one update of the parameter vector $\( \theta = \theta - \eta \nabla_\theta\chi^2 \)$

• Use a learning hyperparameter eta = 0.1

• Start from an initial guess \(\theta_0 = (0,0)\), or some other choice.

### EXERCISE
#
# (a) Implement batch gradient descent as described above and find the best fit parameters
# (b) Print the best fit parameters
# (c) Plot the fitted linear model together with the data
#

# EXERCISE
#
# Compute a learning curve. How fast does the MSE decrease as a function of iteration number?

# Explore the learning rates with different learning hyperparameter 
# For example: eta = 0.001, 0.01, 0.1, 0.5
#
# You can also try to use a learning schedule such that the rate decreases with iteration number 
# (taking smaller and smaller steps)
# A popular choice is $\eta(t) = eta_0  t_0 / (t+t_0)$, where $\eta_0$ and $t_0$ are learning hyperparameters.

10.3.5. Advanced topic: Linear regression using Stochastic Gradient Descent

Implement an algorithm that performs stochastic gradient descent.

By convention we proceed in epochs. All data will be used once every epoch. The gradient update is performed iteratively for each training sample and the data is usually shuffled for each new epoch to prevent cycles. We will use much fewer epochs that we did for the batch gradient descent approach, thereby losing accuracy but gaining speed.

It is common to use a learning schedule in which we first take large steps, and then shorter and shorter ones. This approach is similar in spirit to the optimization algorithm simulated annealing.

m = len(X)

n_epochs = 50

# We will store the iterations in a list to monitor the convergence
theta_path_sgd = []

def learning_schedule(t):
    eta0, t0 = 1, 50  # learning schedule hyperparameters
    # implement proper function here
    pass

# random initialization
np.random.seed(42)
theta = np.random.randn(2,1)  

for epoch in range(n_epochs):
    random_indices=np.arange(m)
    # shuffle the order of indices randomly
    np.random.shuffle(random_indices)
    # Loop through data and perform the stochastic gradient descent.

# Plot the convergence of optimal parameters.
# Possibly compare with BGD and the solution to the normal equation.

# The Stochastic Gradient Descent optimizer is built-in scikit-learn
from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(max_iter=50, tol=1e-5, penalty=None, eta0=0.1, random_state=42)
sgd_reg.fit(X, y.ravel())

SGDRegressor(eta0=0.1, max_iter=50, penalty=None, random_state=42, tol=1e-05)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

sgd_reg.intercept_, sgd_reg.coef_

(array([2.26548882]), array([3.80567955]))

10.3.6. Advanced topic: Mini-batch gradient descent

Instead of computing the gradients based on the full training set (as in batch gradient descent, BGD), or based on a single instance (as in stochastic gradient descent, SGD) there is a third alternative: mini-batch gradient descent, MGD:

• For each iteration, use a random subset (a mini-batch) of the training data (say 10-25%) and compute the gradient based on that.

• An epoch is defined as several such steps after which you have used all instances of training data. E.g. if the mini-batch size is 10% of all instances, then 10 non-overlapping steps would exhaust all data and would correspond to one epoch.

• It is common to use a learning schedule with smaller and smaller steps.

The advantages of the mini-batch gradient descent are the following:

• The convergence is less erratic than with SGD, while still being much faster than full BGD.

• The gradient computation for a mini-batch is a matrix operation, which means that the algorithm can get a performance boost from hardware optimization, especially when using GPUs.

10.4. Important note on gradient descent methods

Gradient descent is a general optimization algorithm. However, there are several important issues that should be known before using it:

1. It requires the computation of partial derivatives of the cost function. This is straight-forward for the linear regression method, but can be difficult for other models. The use of automatic differentiation is very popular in the ML community,and is well worth exploring.

2. In principle, gradient descent works well for convex cost functions, i.e. where the gradient will eventually direct you to the position of the global minimum. Again, the linear regression problem is favorable because you can show that the cost function has that property. However, most cost functions—in particular in many dimensions—correspond to very complicated surfaces with many local minima. In those cases, gradient descent is often not a good method.

10.5. Linear regression with higher-degree polynomials

Any basis-expansion model that is linear in the parameters is a linear model.

Here we will consider a case in which we use at least quadratic features in our model. This is sometimes known as polynomial regression, but this label is really a misnomer since it is still a linear regression problem.

import numpy as np

# Let us generate some cubic data
m = 100
minX = -3
maxX = 3
np.random.seed(1)
X = (maxX-minX) * np.random.rand(m, 1) + minX
# up to cubic features, plus random noise
theta_true = np.array([2, 1, 0.5, -0.25])
eps_noise = 1.
y = eps_noise * np.random.randn(m, 1)
for order in range(len(theta_true)):
    y += theta_true[order] * X**order

# Plot the data. It is pretty tricky to see the features
fig,ax = plt.subplots(1,1)

ax.plot(X, y, "b.")
ax.set_xlabel("$x_1$")
ax.set_ylabel("$y$");

10.5.1. using scikit-learn

In the following we use the LinearRegression class from scikit-learn. You can equally well construct the design matrix and solve the normal equation explicitly with linear algebra.

The PolynomialFeatures class offers a fast way to construct the design matrix with all features up to a specified degree.

Note that problems with more than one input dimension would then contain mixtures up to that degree. E.g., degree-two polynomial features for a problem with \(x_1\) and \(x_2\) would contain the following features: 1 ; \(x_1\) ; \(x_2\) ; \(x_1^2\) ; \(x_2^2\) ; \(x_1 x_2\)

from sklearn.preprocessing import PolynomialFeatures
poly_features = PolynomialFeatures(degree=3, include_bias=False)
X_poly = poly_features.fit_transform(X)

# Note that the first element of the instance vector is
X[0]

array([-0.49786797])

# while the first element of the design matrix is
# NOTE that we have excluded the bias term (x^0). 
# This term will later be added in the linear regression
X_poly[0]

array([-0.49786797,  0.24787252, -0.12340779])

# perform the fit
lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)

# note the bias term, which is the attribute `intercept_` from the fit
lin_reg.intercept_, lin_reg.coef_

(array([1.96984865]), array([[ 1.02594731,  0.53415279, -0.26334408]]))

# EXERCISE
#
# Compare to the "known" amplitudes of the different features. 
# Why doesn't the fit parameters agree better?
#
#

# Make predictions with the linear regression model on a linear grid of new points
# The PolynomialFeatures.transfom method is useful for preparing the new data for a prediction,
# but it is picky with the shape of the input vector.
X_new=np.linspace(minX, maxX, 100).reshape(100, 1)
X_new_poly = poly_features.transform(X_new)
y_new = lin_reg.predict(X_new_poly)

# EXERCISE
#
# Plot the data and the prediction
#
#

10.6. Over- and underfitting

We will explore fitting to models that have both too many and too few features.

# For these fits we will employ scaling of the data
# We use the built-in StandardScaler to rescale the data to zero mean and unit variance.
# This will make the fit more stable
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

fig,ax = plt.subplots(1,1)

for style, degree in (("g-", 100), ("b--", 3), ("r-.", 1)):
    polybig_features = PolynomialFeatures(degree=degree, include_bias=False)
    std_scaler = StandardScaler()
    lin_reg = LinearRegression()
    # Here we use a Pipeline that assembles several steps that we
    # previously applied sequentially:
    # 1. The data is transformed to the chosen polynomial features.
    # 2. The data is transformed to mean=0 and variance=1 (usually makes it numerically more stable)
    # 3. Perform the linear regression fit
    polynomial_regression = Pipeline([
            ("poly_features", polybig_features),
            ("std_scaler", std_scaler),
            ("lin_reg", lin_reg),
        ])
    polynomial_regression.fit(X, y)
    y_newbig = polynomial_regression.predict(X_new)
    ax.plot(X_new, y_newbig, style, label=f'{degree:>3}')


ax.plot(X, y, "b.")
ax.legend(loc="best")
ax.set_xlabel("$x_1$")
ax.set_ylim([-10,30])
ax.set_ylabel("$y$");

• Note how the high-degree polynomial produces a very wiggly curve that tries very hard to go through the training data. The model explodes near the edges where there is no more training data.

• The first degree polynomial, on the other hand, fails to pick up some trends in the data that is clearly there.

# EXERCISE
#
# Which of these models would you label as **overfitting** 
# and which one as **underfitting** the data?

10.6.1. Learning curves

In order to gauge a model’s generalization performance (predictive power) it is common to split the data into a training set and a validation set. We will also see examples of a third set called the test set.

Learning curves are plots of the model’s performance on both the training and the validation sets, measured by some performance metric such as the mean squared error. This measure is plotted as a function of the size of the training set, or alternatively as a function of the training iterations.

10.6.1.1. Performance metric

# EXERCISE
#
# Before moving on, let us first make sure to understand the mean_squared_error metric.
# Perform the training of a first-order polynomial model to the cubic data that we just generated.

lin_reg = LinearRegression()
lin_reg.fit(X, y)
y_predict = lin_reg.predict(X)

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(X, y_predict, "r-",label="Prediction")
ax.plot(X, y, "b.",label="Data");

# EXERCISE
#
# Write your own function that evaluates the mean-squared-error 
# metric on the training set. 

\(\mathrm{MSE} = \frac{1}{m}\sum_{i=1}^m \left(y_i - y_{\mathrm{predict},i} \right)^2\)

# EXERCISE
#
# Then import the built-in convenience function from sckikit-learn 
# for computing the MSE metric
from sklearn.metrics import mean_squared_error
# and use it to compute the same metric. The numbers should agree.

10.6.1.2. Learning curves

# built-in convenience function for computing the MSE metric
from sklearn.metrics import mean_squared_error
# built-in convenience function for splitting data
from sklearn.model_selection import train_test_split

def plot_learning_curves(model, X, y, ax=None):
    # split the data into training and validation sets
    X_train, X_val, y_train, y_val = train_test_split(X, y, train_size=0.7, random_state=42)
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])
        y_train_predict = model.predict(X_train[:m])
        y_val_predict = model.predict(X_val)
        train_errors.append(mean_squared_error(y_train[:m], y_train_predict))
        val_errors.append(mean_squared_error(y_val, y_val_predict))

    if not ax:
        fig,ax = plt.subplots(1,1)
    ax.plot(np.sqrt(train_errors), "r-+", label="train")
    ax.plot(np.sqrt(val_errors), "b-", label="validation")
    ax.legend(loc="best")
    ax.set_xlabel("Training set size")
    ax.set_ylabel("MSE")

Let us use a first-order polynomial to model the training data and plot the learning curve. Recall that a low mean-square error implies that the model predicts the data very well.

lin_reg = LinearRegression()
fig,ax = plt.subplots(1,1)
plot_learning_curves(lin_reg, X, y, ax=ax)
ax.set_ylim([0,4]);

Several features deserves to be mentioned:

1. The performance on the training set starts at zero when only 1-2 data are in the training set.

2. The error on the training set then increases steadily as more data is added.

3. It finally reaches a plateau.

4. The validation error is initially very high, but reaches a plateau that is very close to the training error.

# EXERCISE
#
# This features in a learning curve are typical for a model that underfits. 
# Can you explain / understand why that is?

Now let us try a very high degree polynomial, which should be overfitting the data.

from sklearn.pipeline import Pipeline

polynomial_regression = Pipeline([
        ("poly_features", PolynomialFeatures(degree=15, include_bias=False)),
        ("lin_reg", LinearRegression()),
    ])

fig,ax = plt.subplots(1,1)
plot_learning_curves(polynomial_regression, X, y, ax=ax)
ax.set_ylim([0,4]);

These learning curves are similar to the underfitting model; but there are some important differences:

1. The training error is much smaller than with the linear model.

2. There is no clear plateau.

3. There is a gap between the curves, which implies that the model performs significantly better on the training data than on the validation set.

Both these examples that we have just studied demonstrate the so called bias-variance tradeoff.

• A high bias model has a relatively large error, most probably due to wrong assumptions about the data features.

• A high variance model is excessively sensitive to small variations in the training data.

• The irreducible error is due to the noisiness of the data itself. It can only be reduced by obtaining better data.

We seek a more systematic way of distinguishing between under- and overfitting models, and for quantification of the different kinds of errors.

We will find that Bayesian statistics has the promise to deliver on that ultimate goal.

First, however, we study a common approach to avoid overfitting–namely regularization. We will later provide a Bayesian interpretation of this approach.

10.7. Regularized models

10.7.1. Ridge regression

from sklearn.linear_model import Ridge

# Let's generate some noisy data with a linear feature.
np.random.seed(1)
m = 20
X = 3 * np.random.rand(m, 1)
y = 1 + 0.5 * X + 0.5 * np.random.randn(m, 1) 
X_new = np.linspace(0, 3, 100).reshape(100, 1)

def train_ridge_model(X_train, y_train, alpha, X_predict=None, degree=1, **model_kargs):
    model = Ridge(alpha, **model_kargs) if alpha > 0 else LinearRegression()
    model = Pipeline([
        ("poly_features", PolynomialFeatures(degree=degree, include_bias=False)),
        ("std_scaler", StandardScaler()),
        ("regul_reg", model),
        ])
    model.fit(X_train, y_train)
    if not len(X_predict):
        X_predict=X_train
    return model.predict(X_predict)

fig,axs = plt.subplots(1,2,figsize=(8,5))

alphas=(0, 1e-5,10, 100)
for i,degree in enumerate((1,10)):
    ax = axs[i]
    for alpha, style in zip(alphas, ("b-", "k-.", "g--", "r:")):
        y_new_regul = train_ridge_model(X, y, alpha, X_predict=X_new, degree=degree, random_state=42)
        ax.plot(X_new, y_new_regul, style, label=r'$\alpha={}$'.format(alpha))
    ax.plot(X, y, "b.")
    ax.legend(loc="upper left")
    ax.set_xlabel("$x_1$")
    ax.axis([0, 3, 0, 4])

axs[0].set_ylabel("$y$");