Python Code can be found at https://github.uio.no/larsbog/FYS-STK4155
Learning goals¶
After completing these exercises, you will know how to
- Take more derivatives of simple products between vectors and matrices
- Implement Ridge regression using the analytical expressions
- Scale data appropriately for linear regression
- Evaluate a model across two different hyperparameters
Exercise 1 - Choice of model and degrees of freedom¶
a) How many degrees of freedom does an OLS model fit to the features $x, x^2, x^3$ and the intercept have?
A ordinary least squares model fit with $n$ features has $n$ free parameters. In this case $n=4$, as the intercept can be interpreted as a separate feature.
b) Why is it bad for a model to have too many degrees of freedom?
c) Why is it bad for a model to have too few degrees of freedom?
d) Read chapter 3.4.1 of Hastie et al.'s book. What is the expression for the effective degrees of freedom of the ridge regression fit?
$$ n_\text{effective} = \sum_{i=1}^{n_\text{features}} \frac{d_i^2}{d_i^2 + \lambda} $$
where $d_i$ is the singular value of variable $i$.
e) Why might we want to use Ridge regression instead of OLS?
f) Why migth we want to use OLS instead of Ridge regression?
Exercise 2 - Deriving the expression for Ridge Regression¶
The aim here is to derive the expression for the optimal parameters using Ridge regression.
The expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, was given by the optimization problem
$$ {\displaystyle \min_{\boldsymbol{\beta}\in {\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)^T\left(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\right)\right\}. $$
By minimizing the above equation with respect to the parameters $\boldsymbol{\beta}$ we could then obtain an analytical expression for the parameters $\boldsymbol{\hat\beta_{OLS}}$.
We can add a regularization parameter $\lambda$ by defining a new cost function to be optimized, that is
$$ {\displaystyle \min_{\boldsymbol{\beta}\in {\mathbb{R}}^{p}}}\frac{1}{n}\vert\vert \boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\vert\vert_2^2+\lambda\vert\vert \boldsymbol{\beta}\vert\vert_2^2 $$
which leads to the Ridge regression minimization problem. (One can require as part of the optimization problem that $\vert\vert \boldsymbol{\beta}\vert\vert_2^2\le t$, where $t$ is a finite number larger than zero. We will not implement that in this course.)
a) Expression for Ridge regression¶
Show that the optimal parameters $$ \hat{\boldsymbol{\beta}}_{\mathrm{Ridge}} = \left(\boldsymbol{X}^T\boldsymbol{X}+\lambda\boldsymbol{I}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}, $$ with $\boldsymbol{I}$ being a $p\times p$ identity matrix.
The ordinary least squares result is $$ \hat{\boldsymbol{\beta}}_{\mathrm{OLS}} = \left(\boldsymbol{X}^T\boldsymbol{X}\right)^{-1}\boldsymbol{X}^T\boldsymbol{y}, $$
Differentiate the objective and set the gradient to zero: $$ \nabla_{\beta}\big(\|y-X\beta\|_2^2+\lambda\|\beta\|_2^2\big) = -2X^{\!T}(y-X\beta)+2\lambda\beta = -2X^{\!T}y+2X^{\!T}X\beta+2\lambda\beta = 0 . $$ This yields the equation $$ \big(X^{\!T}X+\lambda I_p\big)\beta = X^{\!T}y . $$ For $\lambda>0$, $X^{\!T}X+\lambda I_p$ is positive definite (hence invertible), so $$ \hat{\beta}_{\mathrm{Ridge}} =\big(X^{\!T}X+\lambda I_p\big)^{-1}X^{\!T}y $$ which reduces to OLS, $\hat\beta_{\mathrm{OLS}}=(X^{\!T}X)^{-1}X^{\!T}y$, when $\lambda=0$ (assuming $X^{\!T}X$ invertible).
Exercise 3 - Scaling data¶
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
n = 100
x = np.linspace(-3, 3, n)
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2) + np.random.normal(0, 0.1)
a) Adapt your function from last week to only include the intercept column if the boolean argument intercept is set to true.
def polynomial_features(x, p, intercept=False):
n = len(x)
if intercept:
P = np.arange(p+1)
else:
P = np.arange(1, p+1)
X = np.power(x[:, np.newaxis], P)
return X
b) Split your data into training and test data(80/20 split)
X = polynomial_features(x, 3)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
x_train = X_train[:, 0] # These are used for plotting
x_test = X_test[:, 0] # These are used for plotting
c) Scale your design matrix with the sklearn standard scaler, though based on the mean and standard deviation of the training data only.
scaler = StandardScaler()
scaler.fit(X_train)
X_train_s = scaler.transform(X_train)
X_test_s = scaler.transform(X_test)
y_offset = np.mean(y_train)
x_s = scaler.transform(X)[:, 0] # These are used for plotting
Exercise 4 - Implementing Ridge Regression¶
a) Implement a function for computing the optimal Ridge parameters using the expression from 2a).
def Ridge_parameters(X, y, lam = 0.01):
# Assumes X is scaled and has no intercept column
return np.linalg.inv(X.T @ X + lam * np.eye(X.shape[1])) @ X.T @ y
beta = Ridge_parameters(X_train_s, y_train)
b) Fit a model to the data, and plot the prediction using both the training and test x-values extracted before scaling, and the y_offset.
plt.plot(x, y)
plt.scatter(x_train, X_train_s @ beta + y_offset)
plt.scatter(x_test, X_test_s @ beta + y_offset)
<matplotlib.collections.PathCollection at 0x7f4f71186490>
Exercise 4 - Testing multiple hyperparameters¶
a) Compute the MSE of your ridge model for polynomials of degrees 1 to 5 with lambda set to 0.01. Plot the MSE as a function of polynomial degree.
from sklearn.metrics import mean_squared_error
def evaluate_model(p, lam, intercept=False):
X_train = polynomial_features(x_train, p, intercept)
X_test = polynomial_features(x_test, p, intercept)
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)
beta = Ridge_parameters(X_train, y_train, lam)
y_pred_train = X_train @ beta
y_pred = X_test @ beta
return mean_squared_error(y_train, y_pred_train), mean_squared_error(y_test, y_pred)
mse_values = []
for p in range(1, 6):
mse_train, mse_test = evaluate_model(p, 0.01, True)
mse_values.append((mse_train, mse_test))
mse_values = np.array(mse_values)
plt.plot(range(1, 6), mse_values[:, 0], marker='o', label='Train')
plt.plot(range(1, 6), mse_values[:, 1], marker='o', label='Test')
plt.xlabel('Polynomial Degree')
plt.ylabel('Mean Squared Error')
plt.legend()
plt.show()
b) Compute the MSE of your ridge model for a polynomial with degree 3, and with lambdas from $10^{-1}$ to $10^{-5}$ on a logarithmic scale. Plot the MSE as a function of lambda.
mse_values = []
lam_values = np.logspace(-5, -1, 10)
for lam in lam_values:
mse_train, mse_test = evaluate_model(p, lam, True)
mse_values.append((mse_train, mse_test))
mse_values = np.array(mse_values)
plt.plot(lam_values, mse_values[:, 0], marker='o', label='Train')
plt.plot(lam_values, mse_values[:, 1], marker='o', label='Test')
plt.xscale('log')
plt.xlabel('Regularization Strength (lambda)')
plt.ylabel('Mean Squared Error')
plt.legend()
plt.show()
c) Compute the MSE of your ridge model for polynomials of degrees 1 to 5, and with lambdas from $10^{-1}$ to $10^{-5}$ on a logarithmic scale. Plot the MSE as a function of polynomial degree and lambda using a heatmap.
p_values = np.arange(1, 6)
lam_values = np.logspace(-5, -1, 8)
mse_values = np.zeros((len(p_values), len(lam_values), 2))
for i, p in enumerate(p_values):
for j, lam in enumerate(lam_values):
mse_train, mse_test = evaluate_model(p, lam=lam, intercept=True)
mse_values[i, j, 0] = mse_train
mse_values[i, j, 1] = mse_test
mse_values = np.array(mse_values)
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 5))
heatmap = ax1.imshow(mse_values[:, :, 0], aspect='auto', interpolation='nearest', cmap='viridis')
ax1.set_title('Train MSE')
fig.colorbar(heatmap, ax=ax1)
heatmap = ax2.imshow(mse_values[:, :, 1], aspect='auto', interpolation='nearest', cmap='viridis')
ax2.set_title('Test MSE')
fig.colorbar(heatmap, ax=ax2)
for ax in (ax1, ax2):
ax.set_xlabel('Regularization Strength (lambda)')
ax.set_ylabel('Polynomial Degree')
ax.set_xticks(np.arange(len(lam_values)))
ax.set_xticklabels([f"{lam:.2e}" for lam in lam_values], rotation=45)
ax.set_yticks(np.arange(len(p_values)))
ax.set_yticklabels([f"Degree {p}" for p in p_values])
plt.tight_layout()
plt.show()
# Optional bigger range
p_values = np.arange(1, 20)
lam_values = np.logspace(-5, 0, 19)
mse_values = np.zeros((len(p_values), len(lam_values), 2))
for i, p in enumerate(p_values):
for j, lam in enumerate(lam_values):
mse_train, mse_test = evaluate_model(p, lam=lam, intercept=True)
mse_values[i, j, 0] = mse_train
mse_values[i, j, 1] = mse_test
mse_values = np.array(mse_values)
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 5))
heatmap = ax1.imshow(mse_values[:, :, 0], aspect='auto', interpolation='nearest', cmap='viridis')
ax1.set_title('Train MSE')
fig.colorbar(heatmap, ax=ax1)
heatmap = ax2.imshow(mse_values[:, :, 1], aspect='auto', interpolation='nearest', cmap='viridis')
ax2.set_title('Test MSE')
fig.colorbar(heatmap, ax=ax2)
for ax in (ax1, ax2):
ax.set_xlabel('Regularization Strength (lambda)')
ax.set_ylabel('Polynomial Degree')
ax.set_xticks(np.arange(len(lam_values)))
ax.set_xticklabels([f"{lam:.2e}" for lam in lam_values], rotation=90)
ax.set_yticks(np.arange(len(p_values)))
ax.set_yticklabels([f"Degree {p}" for p in p_values])
plt.tight_layout()
plt.show()
