addition
This commit is contained in:
Morten Hjorth-Jensen
@@ -2320,3 +2320,505 @@ We can redefine $\lambda$ to absorb the constant $n/2$ and we rewrite the last e
|
||||
This equation does not lead to a nice analytical equation as in either Ridge regression or ordinary least squares. This equation can however be solved by using standard convex optimization algorithms using for example the Python package "CVXOPT":"https://cvxopt.org/". We will discuss how to code LASSO regression next week, when we have introduced gradient methods.
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== Material for exercises week 35 =====
|
||||
|
||||
|
||||
===== Important technicalities: More on Rescaling data =====
|
||||
|
||||
|
||||
When you are comparing your own code with for example _Scikit-Learn_'s
|
||||
library, there are some technicalities to keep in mind. The examples
|
||||
here demonstrate some of these aspects with potential pitfalls.
|
||||
|
||||
The discussion here focuses on the role of the intercept, how we can
|
||||
set up the design matrix, what scaling we should use and other topics
|
||||
which tend confuse us.
|
||||
|
||||
The intercept can be interpreted as the expected value of our
|
||||
target/output variables when all other predictors are set to zero.
|
||||
Thus, if we cannot assume that the expected outputs/targets are zero
|
||||
when all predictors are zero (the columns in the design matrix), it
|
||||
may be a bad idea to implement a model which penalizes the intercept.
|
||||
Furthermore, in for example Ridge and Lasso regression, the default solutions
|
||||
from the library _Scikit-Learn_ (when not shrinking $\beta_0$) for the unknown parameters
|
||||
$\bm{\beta}$, are derived under the assumption that both $\bm{y}$ and
|
||||
$\bm{X}$ are zero centered, that is we subtract the mean values.
|
||||
|
||||
|
||||
If our predictors represent different scales, then it is important to
|
||||
standardize the design matrix $\bm{X}$ by subtracting the mean of each
|
||||
column from the corresponding column and dividing the column with its
|
||||
standard deviation. Most machine learning libraries do this as a default. This means that if you compare your code with the results from a given library,
|
||||
the results may differ.
|
||||
|
||||
The
|
||||
"Standardscaler":"https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html"
|
||||
function in _Scikit-Learn_ does this for us. For the data sets we
|
||||
have been studying in our various examples, the data are in many cases
|
||||
already scaled and there is no need to scale them. You as a user of different machine learning algorithms, should always perform a
|
||||
survey of your data, with a critical assessment of them in case you need to scale the data.
|
||||
|
||||
If you need to scale the data, not doing so will give an *unfair*
|
||||
penalization of the parameters since their magnitude depends on the
|
||||
scale of their corresponding predictor.
|
||||
|
||||
The _Scikit-Learn_ site URL:"https://scikit-learn.org/stable/auto_examples/preprocessing/plot_all_scaling.html#plot-all-scaling-standard-scaler-section" has a good discussion of different ways of preprocessing data.
|
||||
|
||||
Suppose as an example that you
|
||||
you have an input variable given by the heights of different persons.
|
||||
Human height might be measured in inches or meters or
|
||||
kilometers. If measured in kilometers, a standard linear regression
|
||||
model with this predictor would probably give a much bigger
|
||||
coefficient term, than if measured in millimeters.
|
||||
This can clearly lead to problems in evaluating the cost/loss functions.
|
||||
|
||||
|
||||
|
||||
Keep in mind that when you transform your data set before training a model, the same transformation needs to be done
|
||||
on your eventual new data set before making a prediction. If we translate this into a Python code, it would could be implemented as
|
||||
|
||||
!bc pycod
|
||||
"""
|
||||
#Model training, we compute the mean value of y and X
|
||||
y_train_mean = np.mean(y_train)
|
||||
X_train_mean = np.mean(X_train,axis=0)
|
||||
X_train = X_train - X_train_mean
|
||||
y_train = y_train - y_train_mean
|
||||
|
||||
# The we fit our model with the training data
|
||||
trained_model = some_model.fit(X_train,y_train)
|
||||
|
||||
|
||||
#Model prediction, we need also to transform our data set used for the prediction.
|
||||
X_test = X_test - X_train_mean #Use mean from training data
|
||||
y_pred = trained_model(X_test)
|
||||
y_pred = y_pred + y_train_mean
|
||||
"""
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
Let us try to understand what this may imply mathematically when we
|
||||
subtract the mean values, also known as *zero centering*. For
|
||||
simplicity, we will focus on ordinary regression, as done in the above example.
|
||||
|
||||
The cost/loss function for regression is
|
||||
!bt
|
||||
\[
|
||||
C(\beta_0, \beta_1, ... , \beta_{p-1}) = \frac{1}{n}\sum_{i=0}^{n} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij}\beta_j\right)^2,.
|
||||
\]
|
||||
!et
|
||||
|
||||
Recall also that we use the squared value. This expression can lead to an
|
||||
increased penalty for higher differences between predicted and
|
||||
output/target values.
|
||||
|
||||
What we have done is to single out the $\beta_0$ term in the
|
||||
definition of the mean squared error (MSE). The design matrix $X$
|
||||
does in this case not contain any intercept column. When we take the
|
||||
derivative with respect to $\beta_0$, we want the derivative to obey
|
||||
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial C}{\partial \beta_j} = 0,
|
||||
\]
|
||||
!et
|
||||
|
||||
for all $j$. For $\beta_0$ we have
|
||||
|
||||
!bt
|
||||
\[
|
||||
\frac{\partial C}{\partial \beta_0} = -\frac{2}{n}\sum_{i=0}^{n-1} \left(y_i - \beta_0 - \sum_{j=1}^{p-1} X_{ij} \beta_j\right).
|
||||
\]
|
||||
!et
|
||||
Multiplying away the constant $2/n$, we obtain
|
||||
!bt
|
||||
\[
|
||||
\sum_{i=0}^{n-1} \beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} \sum_{j=1}^{p-1} X_{ij} \beta_j.
|
||||
\]
|
||||
!et
|
||||
|
||||
Let us specialize first to the case where we have only two parameters $\beta_0$ and $\beta_1$.
|
||||
Our result for $\beta_0$ simplifies then to
|
||||
!bt
|
||||
\[
|
||||
n\beta_0 = \sum_{i=0}^{n-1}y_i - \sum_{i=0}^{n-1} X_{i1} \beta_1.
|
||||
\]
|
||||
!et
|
||||
We obtain then
|
||||
!bt
|
||||
\[
|
||||
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \beta_1\frac{1}{n}\sum_{i=0}^{n-1} X_{i1}.
|
||||
\]
|
||||
!et
|
||||
If we define
|
||||
!bt
|
||||
\[
|
||||
\mu_{\bm{x}_1}=\frac{1}{n}\sum_{i=0}^{n-1} X_{i1},
|
||||
\]
|
||||
!et
|
||||
and the mean value of the outputs as
|
||||
!bt
|
||||
\[
|
||||
\mu_y=\frac{1}{n}\sum_{i=0}^{n-1}y_i,
|
||||
\]
|
||||
!et
|
||||
we have
|
||||
!bt
|
||||
\[
|
||||
\beta_0 = \mu_y - \beta_1\mu_{\bm{x}_1}.
|
||||
\]
|
||||
!et
|
||||
In the general case with more parameters than $\beta_0$ and $\beta_1$, we have
|
||||
!bt
|
||||
\[
|
||||
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=1}^{p-1} X_{ij}\beta_j.
|
||||
\]
|
||||
!et
|
||||
|
||||
We can rewrite the latter equation as
|
||||
!bt
|
||||
\[
|
||||
\beta_0 = \frac{1}{n}\sum_{i=0}^{n-1}y_i - \sum_{j=1}^{p-1} \mu_{\bm{x}_j}\beta_j,
|
||||
\]
|
||||
!et
|
||||
where we have defined
|
||||
!bt
|
||||
\[
|
||||
\mu_{\bm{x}_j}=\frac{1}{n}\sum_{i=0}^{n-1} X_{ij},
|
||||
\]
|
||||
!et
|
||||
the mean value for all elements of the column vector $\bm{x}_j$.
|
||||
|
||||
|
||||
|
||||
Replacing $y_i$ with $y_i - y_i - \overline{\bm{y}}$ and centering also our design matrix results in a cost function (in vector-matrix disguise)
|
||||
!bt
|
||||
\[
|
||||
C(\boldsymbol{\beta}) = (\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta})^T(\boldsymbol{\tilde{y}} - \tilde{X}\boldsymbol{\beta}).
|
||||
\]
|
||||
!et
|
||||
|
||||
|
||||
|
||||
If we minimize with respect to $\bm{\beta}$ we have then
|
||||
|
||||
!bt
|
||||
\[
|
||||
\hat{\bm{\beta}} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\boldsymbol{\tilde{y}},
|
||||
\]
|
||||
!et
|
||||
|
||||
where $\boldsymbol{\tilde{y}} = \boldsymbol{y} - \overline{\bm{y}}$
|
||||
and $\tilde{X}_{ij} = X_{ij} - \frac{1}{n}\sum_{k=0}^{n-1}X_{kj}$.
|
||||
|
||||
For Ridge regression we need to add $\lambda \boldsymbol{\beta}^T\boldsymbol{\beta}$ to the cost function and get then
|
||||
!bt
|
||||
\[
|
||||
\hat{\bm{\beta}} = (\tilde{X}^T\tilde{X} + \lambda I)^{-1}\tilde{X}^T\boldsymbol{\tilde{y}}.
|
||||
\]
|
||||
!et
|
||||
|
||||
What does this mean? And why do we insist on all this? Let us look at some examples.
|
||||
|
||||
|
||||
This code shows a simple first-order fit to a data set using the above transformed data, where we consider the role of the intercept first, by either excluding it or including it (*code example thanks to Øyvind Sigmundson Schøyen*). Here our scaling of the data is done by subtracting the mean values only.
|
||||
Note also that we do not split the data into training and test.
|
||||
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from sklearn.linear_model import LinearRegression
|
||||
|
||||
|
||||
np.random.seed(2021)
|
||||
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
|
||||
|
||||
def fit_beta(X, y):
|
||||
return np.linalg.pinv(X.T @ X) @ X.T @ y
|
||||
|
||||
|
||||
true_beta = [2, 0.5, 3.7]
|
||||
|
||||
x = np.linspace(0, 1, 11)
|
||||
y = np.sum(
|
||||
np.asarray([x ** p * b for p, b in enumerate(true_beta)]), axis=0
|
||||
) + 0.1 * np.random.normal(size=len(x))
|
||||
|
||||
degree = 3
|
||||
X = np.zeros((len(x), degree))
|
||||
|
||||
# Include the intercept in the design matrix
|
||||
for p in range(degree):
|
||||
X[:, p] = x ** p
|
||||
|
||||
beta = fit_beta(X, y)
|
||||
|
||||
# Intercept is included in the design matrix
|
||||
skl = LinearRegression(fit_intercept=False).fit(X, y)
|
||||
|
||||
print(f"True beta: {true_beta}")
|
||||
print(f"Fitted beta: {beta}")
|
||||
print(f"Sklearn fitted beta: {skl.coef_}")
|
||||
ypredictOwn = X @ beta
|
||||
ypredictSKL = skl.predict(X)
|
||||
print(f"MSE with intercept column")
|
||||
print(MSE(y,ypredictOwn))
|
||||
print(f"MSE with intercept column from SKL")
|
||||
print(MSE(y,ypredictSKL))
|
||||
|
||||
|
||||
plt.figure()
|
||||
plt.scatter(x, y, label="Data")
|
||||
plt.plot(x, X @ beta, label="Fit")
|
||||
plt.plot(x, skl.predict(X), label="Sklearn (fit_intercept=False)")
|
||||
|
||||
|
||||
# Do not include the intercept in the design matrix
|
||||
X = np.zeros((len(x), degree - 1))
|
||||
|
||||
for p in range(degree - 1):
|
||||
X[:, p] = x ** (p + 1)
|
||||
|
||||
# Intercept is not included in the design matrix
|
||||
skl = LinearRegression(fit_intercept=True).fit(X, y)
|
||||
|
||||
# Use centered values for X and y when computing coefficients
|
||||
y_offset = np.average(y, axis=0)
|
||||
X_offset = np.average(X, axis=0)
|
||||
|
||||
beta = fit_beta(X - X_offset, y - y_offset)
|
||||
intercept = np.mean(y_offset - X_offset @ beta)
|
||||
|
||||
print(f"Manual intercept: {intercept}")
|
||||
print(f"Fitted beta (without intercept): {beta}")
|
||||
print(f"Sklearn intercept: {skl.intercept_}")
|
||||
print(f"Sklearn fitted beta (without intercept): {skl.coef_}")
|
||||
ypredictOwn = X @ beta
|
||||
ypredictSKL = skl.predict(X)
|
||||
print(f"MSE with Manual intercept")
|
||||
print(MSE(y,ypredictOwn+intercept))
|
||||
print(f"MSE with Sklearn intercept")
|
||||
print(MSE(y,ypredictSKL))
|
||||
|
||||
plt.plot(x, X @ beta + intercept, "--", label="Fit (manual intercept)")
|
||||
plt.plot(x, skl.predict(X), "--", label="Sklearn (fit_intercept=True)")
|
||||
plt.grid()
|
||||
plt.legend()
|
||||
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
The intercept is the value of our output/target variable
|
||||
when all our features are zero and our function crosses the $y$-axis (for a one-dimensional case).
|
||||
|
||||
Printing the MSE, we see first that both methods give the same MSE, as
|
||||
they should. However, when we move to for example Ridge regression,
|
||||
the way we treat the intercept may give a larger or smaller MSE,
|
||||
meaning that the MSE can be penalized by the value of the
|
||||
intercept. Not including the intercept in the fit, means that the
|
||||
regularization term does not include $\beta_0$. For different values
|
||||
of $\lambda$, this may lead to different MSE values.
|
||||
|
||||
To remind the reader, the regularization term, with the intercept in Ridge regression, is given by
|
||||
!bt
|
||||
\[
|
||||
\lambda \vert\vert \bm{\beta} \vert\vert_2^2 = \lambda \sum_{j=0}^{p-1}\beta_j^2,
|
||||
\]
|
||||
!et
|
||||
but when we take out the intercept, this equation becomes
|
||||
!bt
|
||||
\[
|
||||
\lambda \vert\vert \bm{\beta} \vert\vert_2^2 = \lambda \sum_{j=1}^{p-1}\beta_j^2.
|
||||
\]
|
||||
!et
|
||||
|
||||
For Lasso regression we have
|
||||
!bt
|
||||
\[
|
||||
\lambda \vert\vert \bm{\beta} \vert\vert_1 = \lambda \sum_{j=1}^{p-1}\vert\beta_j\vert.
|
||||
\]
|
||||
!et
|
||||
|
||||
It means that, when scaling the design matrix and the outputs/targets,
|
||||
by subtracting the mean values, we have an optimization problem which
|
||||
is not penalized by the intercept. The MSE value can then be smaller
|
||||
since it focuses only on the remaining quantities. If we however bring
|
||||
back the intercept, we will get a MSE which then contains the
|
||||
intercept.
|
||||
|
||||
|
||||
Armed with this wisdom, we attempt first to simply set the intercept equal to _False_ in our implementation of Ridge regression for our well-known vanilla data set.
|
||||
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.model_selection import train_test_split
|
||||
from sklearn import linear_model
|
||||
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
|
||||
|
||||
# A seed just to ensure that the random numbers are the same for every run.
|
||||
# Useful for eventual debugging.
|
||||
np.random.seed(3155)
|
||||
|
||||
n = 100
|
||||
x = np.random.rand(n)
|
||||
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
|
||||
|
||||
Maxpolydegree = 20
|
||||
X = np.zeros((n,Maxpolydegree))
|
||||
#We include explicitely the intercept column
|
||||
for degree in range(Maxpolydegree):
|
||||
X[:,degree] = x**degree
|
||||
# We split the data in test and training data
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
p = Maxpolydegree
|
||||
I = np.eye(p,p)
|
||||
# Decide which values of lambda to use
|
||||
nlambdas = 6
|
||||
MSEOwnRidgePredict = np.zeros(nlambdas)
|
||||
MSERidgePredict = np.zeros(nlambdas)
|
||||
lambdas = np.logspace(-4, 2, nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lambdas[i]
|
||||
OwnRidgeBeta = np.linalg.pinv(X_train.T @ X_train+lmb*I) @ X_train.T @ y_train
|
||||
# Note: we include the intercept column and no scaling
|
||||
RegRidge = linear_model.Ridge(lmb,fit_intercept=False)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
# and then make the prediction
|
||||
ytildeOwnRidge = X_train @ OwnRidgeBeta
|
||||
ypredictOwnRidge = X_test @ OwnRidgeBeta
|
||||
ytildeRidge = RegRidge.predict(X_train)
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
|
||||
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
|
||||
print("Beta values for own Ridge implementation")
|
||||
print(OwnRidgeBeta)
|
||||
print("Beta values for Scikit-Learn Ridge implementation")
|
||||
print(RegRidge.coef_)
|
||||
print("MSE values for own Ridge implementation")
|
||||
print(MSEOwnRidgePredict[i])
|
||||
print("MSE values for Scikit-Learn Ridge implementation")
|
||||
print(MSERidgePredict[i])
|
||||
|
||||
# Now plot the results
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'r', label = 'MSE own Ridge Test')
|
||||
plt.plot(np.log10(lambdas), MSERidgePredict, 'g', label = 'MSE Ridge Test')
|
||||
|
||||
plt.xlabel('log10(lambda)')
|
||||
plt.ylabel('MSE')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
The results here agree when we force _Scikit-Learn_'s Ridge function to include the first column in our design matrix.
|
||||
We see that the results agree very well. Here we have thus explicitely included the intercept column in the design matrix.
|
||||
What happens if we do not include the intercept in our fit?
|
||||
Let us see how we can change this code by zero centering.
|
||||
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
from sklearn.model_selection import train_test_split
|
||||
from sklearn import linear_model
|
||||
from sklearn.preprocessing import StandardScaler
|
||||
|
||||
def MSE(y_data,y_model):
|
||||
n = np.size(y_model)
|
||||
return np.sum((y_data-y_model)**2)/n
|
||||
# A seed just to ensure that the random numbers are the same for every run.
|
||||
# Useful for eventual debugging.
|
||||
np.random.seed(315)
|
||||
|
||||
n = 100
|
||||
x = np.random.rand(n)
|
||||
y = np.exp(-x**2) + 1.5 * np.exp(-(x-2)**2)
|
||||
|
||||
Maxpolydegree = 20
|
||||
X = np.zeros((n,Maxpolydegree-1))
|
||||
|
||||
for degree in range(1,Maxpolydegree): #No intercept column
|
||||
X[:,degree-1] = x**(degree)
|
||||
|
||||
# We split the data in test and training data
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
|
||||
#For our own implementation, we will need to deal with the intercept by centering the design matrix and the target variable
|
||||
X_train_mean = np.mean(X_train,axis=0)
|
||||
#Center by removing mean from each feature
|
||||
X_train_scaled = X_train - X_train_mean
|
||||
X_test_scaled = X_test - X_train_mean
|
||||
#The model intercept (called y_scaler) is given by the mean of the target variable (IF X is centered)
|
||||
#Remove the intercept from the training data.
|
||||
y_scaler = np.mean(y_train)
|
||||
y_train_scaled = y_train - y_scaler
|
||||
|
||||
p = Maxpolydegree-1
|
||||
I = np.eye(p,p)
|
||||
# Decide which values of lambda to use
|
||||
nlambdas = 6
|
||||
MSEOwnRidgePredict = np.zeros(nlambdas)
|
||||
MSERidgePredict = np.zeros(nlambdas)
|
||||
|
||||
lambdas = np.logspace(-4, 2, nlambdas)
|
||||
for i in range(nlambdas):
|
||||
lmb = lambdas[i]
|
||||
OwnRidgeBeta = np.linalg.pinv(X_train_scaled.T @ X_train_scaled+lmb*I) @ X_train_scaled.T @ (y_train_scaled)
|
||||
intercept_ = y_scaler - X_train_mean@OwnRidgeBeta #The intercept can be shifted so the model can predict on uncentered data
|
||||
#Add intercept to prediction
|
||||
ypredictOwnRidge = X_test_scaled @ OwnRidgeBeta + y_scaler
|
||||
RegRidge = linear_model.Ridge(lmb)
|
||||
RegRidge.fit(X_train,y_train)
|
||||
ypredictRidge = RegRidge.predict(X_test)
|
||||
MSEOwnRidgePredict[i] = MSE(y_test,ypredictOwnRidge)
|
||||
MSERidgePredict[i] = MSE(y_test,ypredictRidge)
|
||||
print("Beta values for own Ridge implementation")
|
||||
print(OwnRidgeBeta) #Intercept is given by mean of target variable
|
||||
print("Beta values for Scikit-Learn Ridge implementation")
|
||||
print(RegRidge.coef_)
|
||||
print('Intercept from own implementation:')
|
||||
print(intercept_)
|
||||
print('Intercept from Scikit-Learn Ridge implementation')
|
||||
print(RegRidge.intercept_)
|
||||
print("MSE values for own Ridge implementation")
|
||||
print(MSEOwnRidgePredict[i])
|
||||
print("MSE values for Scikit-Learn Ridge implementation")
|
||||
print(MSERidgePredict[i])
|
||||
|
||||
|
||||
# Now plot the results
|
||||
plt.figure()
|
||||
plt.plot(np.log10(lambdas), MSEOwnRidgePredict, 'b--', label = 'MSE own Ridge Test')
|
||||
plt.plot(np.log10(lambdas), MSERidgePredict, 'g--', label = 'MSE SL Ridge Test')
|
||||
plt.xlabel('log10(lambda)')
|
||||
plt.ylabel('MSE')
|
||||
plt.legend()
|
||||
plt.show()
|
||||
!ec
|
||||
We see here, when compared to the code which includes explicitely the
|
||||
intercept column, that our MSE value is actually smaller. This is
|
||||
because the regularization term does not include the intercept value
|
||||
$\beta_0$ in the fitting. This applies to Lasso regularization as
|
||||
well. It means that our optimization is now done only with the
|
||||
centered matrix and/or vector that enter the fitting procedure.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user