lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

updated lecture notes

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2023-09-03 13:42:26 +02:00
parent e69df0de9c
commit 22c8178dfb
61 changed files with 5198 additions and 10085 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/src/week36/.ipynb_checkpoints/exercisesweek36-checkpoint.ipynb
+351
View File
@@ -0,0 +1,351 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "d4ee754f",
"metadata": {},
"source": [
"<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)\n",
"doconce format html exercisesweek36.do.txt -->\n",
"<!-- dom:TITLE: Exercises week 36 -->"
]
},
{
"cell_type": "markdown",
"id": "443836cf",
"metadata": {},
"source": [
"# Exercises week 36\n",
"**September 4-8, 2023**\n",
"\n",
"Date: **Deadline is Sunday September 10 at midnight**"
]
},
{
"cell_type": "markdown",
"id": "52627eae",
"metadata": {},
"source": [
"## Overarching aims of the exercises this week\n",
"\n",
"This set of exercises form an important part of the first project. The\n",
"analytical exercises deal with the material covered last week on the\n",
"mathematical interpretations of ordinary least squares and of Ridge\n",
"regression. The numerical exercises can be seen as a continuation of\n",
"exercise 3 from week 35, with the inclusion of Ridge regression. This\n",
"material enters also the discussions of the first project."
]
},
{
"cell_type": "markdown",
"id": "717c36eb",
"metadata": {},
"source": [
"## Exercise 1: Analytical exercises\n",
"\n",
"The aim here is to derive the expression for the optimal parameters\n",
"using Ridge regression. Furthermore, using the singular value\n",
"decomposition, we will analyze the difference between the ordinary\n",
"least squares approach and Ridge regression.\n",
"\n",
"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, that is \n",
"our optimization problem is"
]
},
{
"cell_type": "markdown",
"id": "5c734a1d",
"metadata": {},
"source": [
"$$\n",
"{\\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\\}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "6aaa8b49",
"metadata": {},
"source": [
"or we can state it as"
]
},
{
"cell_type": "markdown",
"id": "afbd0f2c",
"metadata": {},
"source": [
"$$\n",
"{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n",
"{\\mathbb{R}}^{p}}}\\frac{1}{n}\\sum_{i=0}^{n-1}\\left(y_i-\\tilde{y}_i\\right)^2=\\frac{1}{n}\\vert\\vert \\boldsymbol{y}-\\boldsymbol{X}\\boldsymbol{\\beta}\\vert\\vert_2^2,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "ebccf185",
"metadata": {},
"source": [
"where we have used the definition of a norm-2 vector, that is"
]
},
{
"cell_type": "markdown",
"id": "62b078ff",
"metadata": {},
"source": [
"$$\n",
"\\vert\\vert \\boldsymbol{x}\\vert\\vert_2 = \\sqrt{\\sum_i x_i^2}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "70bd79fa",
"metadata": {},
"source": [
"By minimizing the above equation with respect to the parameters\n",
"$\\boldsymbol{\\beta}$ we could then obtain an analytical expression for the\n",
"parameters $\\boldsymbol{\\beta}$.\n",
"\n",
"We can add a regularization parameter $\\lambda$ by\n",
"defining a new cost function to be optimized, that is"
]
},
{
"cell_type": "markdown",
"id": "114bf17c",
"metadata": {},
"source": [
"$$\n",
"{\\displaystyle \\min_{\\boldsymbol{\\beta}\\in\n",
"{\\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\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "2fcf0506",
"metadata": {},
"source": [
"which leads to the Ridge regression minimization problem where we\n",
"require that $\\vert\\vert \\boldsymbol{\\beta}\\vert\\vert_2^2\\le t$, where $t$ is\n",
"a finite number larger than zero."
]
},
{
"cell_type": "markdown",
"id": "ba91302c",
"metadata": {},
"source": [
"### a) Expression for Ridge regression\n",
"\n",
"Show that the optimal parameters"
]
},
{
"cell_type": "markdown",
"id": "4e86f509",
"metadata": {},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}}_{\\mathrm{Ridge}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}+\\lambda\\boldsymbol{I}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "94f90a0e",
"metadata": {},
"source": [
"with $\\boldsymbol{I}$ being a $p\\times p$ identity matrix with the constraint that"
]
},
{
"cell_type": "markdown",
"id": "8a814dcc",
"metadata": {},
"source": [
"$$\n",
"\\sum_{i=0}^{p-1} \\beta_i^2 \\leq t,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "3af42292",
"metadata": {},
"source": [
"with $t$ a finite positive number. \n",
"\n",
"The ordinary least squares result is"
]
},
{
"cell_type": "markdown",
"id": "c8625097",
"metadata": {},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}}_{\\mathrm{OLS}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "3a7f519d",
"metadata": {},
"source": [
"### b) The singular value decomposition\n",
"\n",
"Use the singular value decomposition of an $m\\times n$ matrix $\\boldsymbol{X}$ (our design matrix)"
]
},
{
"cell_type": "markdown",
"id": "a8ae7021",
"metadata": {},
"source": [
"$$\n",
"\\boldsymbol{X}=\\boldsymbol{U}\\boldsymbol{\\Sigma}\\boldsymbol{V}^T,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "e1d65ded",
"metadata": {},
"source": [
"where $\\boldsymbol{U}$ and $\\boldsymbol{V}$ are orthogonal matrices of dimensions\n",
"$m\\times m$ and $n\\times n$, respectively, and $\\boldsymbol{\\Sigma}$ is an\n",
"$m\\times n$ matrix which contains the ingular values only. This material was discussed during the lectures of week 35.\n",
"\n",
"Show that you can write the \n",
"OLS solutions in terms of the eigenvectors (the columns) of the orthogonal matrix $\\boldsymbol{U}$ as"
]
},
{
"cell_type": "markdown",
"id": "e147a177",
"metadata": {},
"source": [
"$$\n",
"\\tilde{\\boldsymbol{y}}_{\\mathrm{OLS}}=\\boldsymbol{X}\\boldsymbol{\\beta} = \\sum_{j=0}^{p-1}\\boldsymbol{u}_j\\boldsymbol{u}_j^T\\boldsymbol{y}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "749540e4",
"metadata": {},
"source": [
"For Ridge regression, show that the corresponding equation is"
]
},
{
"cell_type": "markdown",
"id": "da2f753f",
"metadata": {},
"source": [
"$$\n",
"\\tilde{\\boldsymbol{y}}_{\\mathrm{Ridge}}=\\boldsymbol{X}\\boldsymbol{\\beta}_{\\mathrm{Ridge}} = \\boldsymbol{U\\Sigma V^T}\\left(\\boldsymbol{V}\\boldsymbol{\\Sigma}^2\\boldsymbol{V}^T+\\lambda\\boldsymbol{I} \\right)^{-1}(\\boldsymbol{U\\Sigma V^T})^T\\boldsymbol{y}=\\sum_{j=0}^{p-1}\\boldsymbol{u}_j\\boldsymbol{u}_j^T\\frac{\\sigma_j^2}{\\sigma_j^2+\\lambda}\\boldsymbol{y},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "b322bc70",
"metadata": {},
"source": [
"with the vectors $\\boldsymbol{u}_j$ being the columns of $\\boldsymbol{U}$ from the SVD of the matrix $\\boldsymbol{X}$. \n",
"\n",
"Give an interpretation of the results. Section 3.4 of Hastie et al's textbook gives a good discussion of the results."
]
},
{
"cell_type": "markdown",
"id": "bfaa449c",
"metadata": {},
"source": [
"## Exercise 2: Adding Ridge Regression\n",
"\n",
"This exercise is a continuation of exercise 3 from last week. We will use the same function to\n",
"generate our data set, still staying with a simple function $y(x)$\n",
"which we want to fit using linear regression, but now extending the\n",
"analysis to include the Ridge regression method.\n",
"\n",
"In this exercise you need to include the same elements from last week\n",
"1. Scale your data by subtracting the mean value from each column in the design matrix.\n",
"\n",
"2. Perform a split of the data in training and test set.\n",
"\n",
"The addition to the analysis this time is the introduction of the hyperparameter $\\lambda$ when introducing Ridge regression.\n",
"\n",
"Extend the code from exercise 3 from week 35 to include Ridge regression with the hyperparameter $\\lambda$. The optimal parameters $\\hat{\\beta}$ for Ridge regression can be obtained by matrix inversion in a similar way as done for ordinary least squares. You need to add to your code the following equations"
]
},
{
"cell_type": "markdown",
"id": "ad23d47a",
"metadata": {},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}}_{\\mathrm{Ridge}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}+\\lambda\\boldsymbol{I}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "d49ef4fc",
"metadata": {},
"source": [
"The ordinary least squares result you encoded last week is given by"
]
},
{
"cell_type": "markdown",
"id": "45adbd8d",
"metadata": {},
"source": [
"$$\n",
"\\hat{\\boldsymbol{\\beta}}_{\\mathrm{OLS}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y},\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "ceac24b4",
"metadata": {},
"source": [
"Use these results to compute the mean squared error for ordinary least\n",
"squares and Ridge regression first for a polynomial of degree five\n",
"with $n=100$ data points and five selected values of\n",
"$\\lambda=[0.0001,0.001, 0.01,0.1,1.0]$. Compute thereafter the mean\n",
"squared error for the same values of $\\lambda$ for polynomials of degree ten\n",
"and $15$. Discuss your results for the training MSE and test MSE with\n",
"Ridge regression and ordinary least squares."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
doc/src/week36/exercisesweek36.do.txt
+15 -15
View File
@@ -19,14 +19,14 @@ using Ridge regression. Furthermore, using the singular value
decomposition, we will analyze the difference between the ordinary
least squares approach and 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, that is
our optimization problem is
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
!bt
\[
{\displaystyle \min_{\bm{\beta}\in {\mathbb{R}}^{p}}}\frac{1}{n}\left\{\left(\bm{y}-\bm{X}\bm{\beta}\right)^T\left(\bm{y}-\bm{X}\bm{\beta}\right)\right\}.
\]
!et
or we can state it as
which we can also write as
!bt
\[
{\displaystyle \min_{\bm{\beta}\in
@@ -54,9 +54,9 @@ defining a new cost function to be optimized, that is
\]
!et
which leads to the Ridge regression minimization problem where we
require that $\vert\vert \bm{\beta}\vert\vert_2^2\le t$, where $t$ is
a finite number larger than zero.
which leads to the Ridge regression minimization problem. One can require as part of the optimization problem
that $\vert\vert \bm{\beta}\vert\vert_2^2\le t$, where $t$ is
a finite number larger than zero. We will not implement that here.
=== a) Expression for Ridge regression ===
Show that the optimal parameters
@@ -85,7 +85,7 @@ The ordinary least squares result is
=== b) The singular value decomposition ===
Use the singular value decomposition of an $m\times n$ matrix $\bm{X}$ (our design matrix)
Use the singular value decomposition of an n\times p$ matrix $\bm{X}$ (our design matrix)
!bt
\[
\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^T,
@@ -93,8 +93,8 @@ Use the singular value decomposition of an $m\times n$ matrix $\bm{X}$ (our desi
!et
where $\bm{U}$ and $\bm{V}$ are orthogonal matrices of dimensions
$m\times m$ and $n\times n$, respectively, and $\bm{\Sigma}$ is an
$m\times n$ matrix which contains the ingular values only. This material was discussed during the lectures of week 35.
$n\times n$ and $p\times p$, respectively, and $\bm{\Sigma}$ is an
$n\times p$ matrix which contains the ingular values only. This material was discussed during the lectures of week 35.
Show that you can write the
OLS solutions in terms of the eigenvectors (the columns) of the orthogonal matrix $\bm{U}$ as
@@ -115,7 +115,7 @@ For Ridge regression, show that the corresponding equation is
with the vectors $\bm{u}_j$ being the columns of $\bm{U}$ from the SVD of the matrix $\bm{X}$.
Give an interpretation of the results. Section 3.4 of Hastie et al's textbook gives a good discussion of the results.
Give an interpretation of the results. "Section 3.4 of Hastie et al's textbook gives a good discussion of the above results":"https://link.springer.com/book/10.1007/978-0-387-84858-7".
@@ -123,18 +123,18 @@ Give an interpretation of the results. Section 3.4 of Hastie et al's textbook g
===== Exercise: Adding Ridge Regression =====
This exercise is a continuation of exercise 3 from last week. We will use the same function to
This exercise is a continuation of exercise 3 from week 35, see URL:"https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/exercisesweek35.html". We will use the same function to
generate our data set, still staying with a simple function $y(x)$
which we want to fit using linear regression, but now extending the
analysis to include the Ridge regression method.
In this exercise you need to include the same elements from last week
o Scale your data by subtracting the mean value from each column in the design matrix.
o Perform a split of the data in training and test set.
In this exercise you need to include the same elements from last week, that is
o scale your data by subtracting the mean value from each column in the design matrix.
o perform a split of the data in a training set and a test set.
The addition to the analysis this time is the introduction of the hyperparameter $\lambda$ when introducing Ridge regression.
Extend the code from exercise 3 from week 35 to include Ridge regression with the hyperparameter $\lambda$. The optimal parameters $\hat{\beta}$ for Ridge regression can be obtained by matrix inversion in a similar way as done for ordinary least squares. You need to add to your code the following equations
Extend the code from exercise 3 from "week 35":"https://compphysics.github.io/MachineLearning/doc/LectureNotes/_build/html/exercisesweek35.html" to include Ridge regression with the hyperparameter $\lambda$. The optimal parameters $\hat{\beta}$ for Ridge regression can be obtained by matrix inversion in a similar way as done for ordinary least squares. You need to add to your code the following equations
!bt
\[
\hat{\bm{\beta}}_{\mathrm{Ridge}} = \left(\bm{X}^T\bm{X}+\lambda\bm{I}\right)^{-1}\bm{X}^T\bm{y}.
doc/src/week36/week36.do.txt
+11 -607
View File
@@ -6,15 +6,16 @@ DATE: September 4-8, 2023
!split
===== Plans for week 36 =====
* Summary from last week on discussion of SVD, Ridge and Lasso linear regression.
* Linear Regression and links with Statistics, Resampling methods and presentation of first project
Recommended Reading:
o Lectures on Regression
o Bishop 1.1, 1.2, 2.1, 2.2, 2.3 and 3.1
o Hastie et al chapter 3
o Material for the active learning sessions on Tuesday and Wednesday
o Summary from last week on discussion of SVD, Ridge and Lasso linear regression.
o Recommended Reading: Hastie et al chapter 3, see URL:"https://link.springer.com/book/10.1007/978-0-387-84858-7"
o Presentation and discussion of first project
o Material for the lecture on Thursday September 7
o Linear Regression and links with Statistics, Resampling methods
o Recommended Reading: Hastie et al chapter 3, see URL:"https://link.springer.com/book/10.1007/978-0-387-84858-7"
!split
===== Material for the active learning sessions Tuesday and Wednesday =====
!split
===== Summary from last Week and discussion of SVD, Ridge and Lasso regression with examples =====
@@ -539,7 +540,7 @@ which leads to
\]
!et
Plotting these results ("figure in handwritten notes for week 36":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2021/NotesSeptember9.pdf") shows clearly that Lasso regression suppresses (sets to zero) values of $\beta_i$ for specific values of $\lambda$. Ridge regression reduces on the other hand the values of $\beta_i$ as function of $\lambda$.
Plotting these results shows clearly that Lasso regression suppresses (sets to zero) values of $\beta_i$ for specific values of $\lambda$. Ridge regression reduces on the other hand the values of $\beta_i$ as function of $\lambda$.
@@ -915,506 +916,7 @@ plt.show()
!split
===== To think about, first part =====
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.
!split
===== More thinking =====
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
"Standadscaler":"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.
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.
!split
===== Still thinking =====
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 follows
!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
!split
===== What does centering (subtracting the mean values) mean mathematically? =====
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 since this leads to an increase of the 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
!split
===== Further Manipulations =====
Let us special 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_1=\frac{1}{n}\sum_{i=0}^{n-1} (X_{i1},
\]
!et
and if we define 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_{1}.
\]
!et
In the general case, that is we have 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
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
!split
===== Wrapping it up =====
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.
!split
===== Linear Regression code, Intercept handling first =====
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 (wiothout 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 differeing 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.
!split
===== Code Examples =====
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 (thanks to Stian Bilek for inpouts here).
!split
===== Taking out the mean =====
!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. Note
also that the problem with the intercept occurs mainly in these type
of polynomial fitting problem.
!split
===== Friday September 9 =====
===== Material for lecture Thursday September 7 =====
!split
@@ -2061,101 +1563,3 @@ which is our Lasso cost function!
===== Exercise: mean values and variances in linear regression =====
This exercise deals with various mean values ad variances in linear regression method (here it may be useful to look up chapter 3, equation (3.8) of "Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning, Springer":"https://www.springer.com/gp/book/9780387848570").
The assumption we have made is
that there exists a function $f(\bm{x})$ and a normal distributed error $\bm{\varepsilon}\sim \mathcal{N}(0, \sigma^2)$
which describes our data
!bt
\[
\bm{y} = f(\bm{x})+\bm{\varepsilon}
\]
!et
We then approximate this function with our model from the solution of the linear regression equations (ordinary least squares OLS), that is our
function $f$ is approximated by $\bm{\tilde{y}}$ where we minimized $(\bm{y}-\bm{\tilde{y}})^2$, with
!bt
\[
\bm{\tilde{y}} = \bm{X}\bm{\beta}.
\]
!et
The matrix $\bm{X}$ is the so-called design matrix.
!bsubex
Show that the expectation value of $\bm{y}$ for a given element $i$
!bt
\begin{align*}
\mathbb{E}(y_i) & =\mathbf{X}_{i, \ast} \, \beta,
\end{align*}
!et
and that
its variance is
!bt
\begin{align*} \mbox{Var}(y_i) & = \sigma^2.
\end{align*}
!et
Hence, $y_i \sim \mathcal{N}( \mathbf{X}_{i, \ast} \, \bm{\beta}, \sigma^2)$, that is $\bm{y}$ follows a normal distribution with
mean value $\bm{X}\bm{\beta}$ and variance $\sigma^2$.
!esubex
!bsubex
With the OLS expressions for the parameters $\bm{\beta}$ show that
!bt
\[
\mathbb{E}(\bm{\beta}) = \bm{\beta}.
\]
!et
!esubex
!bsubex
Show finally that the variance of $\bm{\beta}$ is
!bt
\begin{eqnarray*}
\mbox{Var}(\bm{\beta}) & = & \sigma^2 \, (\mathbf{X}^{T} \mathbf{X})^{-1}.
\end{eqnarray*}
!et
!esubex
===== Exercise: Adding Ridge and Lasso Regression =====
This exercise is a continuation of the exercises from week 35.
We will
use the same function to generate our data set, still staying with a
simple function $y(x)$ which we want to fit using linear regression,
but now extending the analysis to include the Ridge and the Lasso
regression methods.
We will thus again generate our own dataset for a function $y(x)$ where
$x \in [0,1]$ and defined by random numbers computed with the uniform
distribution. The function $y$ is a quadratic polynomial in $x$ with
added stochastic noise according to the normal distribution $\cal{N}(0,1)$.
The following simple Python instructions define our $x$ and $y$ values (with 100 data points).
!bc pycod
x = np.random.rand(100)
y = 2.0+5*x*x+0.1*np.random.randn(100)
!ec
!bsubex
Write your own code for the Ridge method (see chapter 3.4 of Hastie *et al.*, equations (3.43) and (3.44)) and compute the parametrization for different values of $\lambda$. Study the dependence on $\lambda$ while also varying the strength of the noise in your expression for $y(x)$.
!esubex
!bsubex
Our next step is to study the variance of the parameters $\beta_1$ and $\beta_2$ (assuming that we are parameterizing our function with a second-order polynomial). We will use standard linear regression and the Ridge regression. You can now opt for either writing your own function or using _Scikit-Learn_ to find the parameters $\beta$. From your results calculate the variance of these parameters (recall that this is equal to the diagonal elements of the matrix $(\hat{X}^T\hat{X})+\lambda\hat{I})^{-1}$). Discuss the results of these variances as functions of $\lambda$. In particular, try to link your discussion with the discussion in Hastie *et al.* and their figures 3.10 and 3.11. _Scikit-Learn_ may not provide the variance of the parameters $\beta$. This needs to be checked. With your own code you can however do so.
!esubex