From 30d664a82bf65bcbf3ef04a165485474ec6fbfdb Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Mon, 13 Sep 2021 14:00:16 +0200 Subject: [PATCH] update --- doc/BookChapters/chapter1.dlog | 27 + doc/BookChapters/chapter1.do.txt | 299 +++++---- doc/BookChapters/chapter2.do.txt | 90 ++- doc/BookChapters/chapter3.do.txt | 311 ++++----- doc/BookChapters/chapter4.do.txt | 712 +++++++++++++++++---- doc/LectureNotes/chapter1.ipynb | 366 ++++++----- doc/LectureNotes/chapter2.ipynb | 110 +++- doc/LectureNotes/chapter3.ipynb | 387 +++++------ doc/LectureNotes/chapter4.ipynb | 1032 +++++++++++++++++++++++++++--- 9 files changed, 2392 insertions(+), 942 deletions(-) diff --git a/doc/BookChapters/chapter1.dlog b/doc/BookChapters/chapter1.dlog index aaeebdc52..6162c8df1 100644 --- a/doc/BookChapters/chapter1.dlog +++ b/doc/BookChapters/chapter1.dlog @@ -7,3 +7,30 @@ Translating doconce text in chapter1.do.txt to ipynb Failed to remove ans_at_end environment Failed to remove sol_at_end environment output in chapter1.ipynb +Translating doconce text in chapter1.do.txt to ipynb +*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax) + +*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments. + +*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments. +Failed to remove ans_at_end environment +Failed to remove sol_at_end environment +output in chapter1.ipynb +Translating doconce text in chapter1.do.txt to ipynb +*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax) + +*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments. + +*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments. +Failed to remove ans_at_end environment +Failed to remove sol_at_end environment +output in chapter1.ipynb +Translating doconce text in chapter1.do.txt to ipynb +*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax) + +*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments. + +*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments. +Failed to remove ans_at_end environment +Failed to remove sol_at_end environment +output in chapter1.ipynb diff --git a/doc/BookChapters/chapter1.do.txt b/doc/BookChapters/chapter1.do.txt index 43367c5de..e822f162e 100644 --- a/doc/BookChapters/chapter1.do.txt +++ b/doc/BookChapters/chapter1.do.txt @@ -7,38 +7,37 @@ -Our emphasis throughout this series of lectures -is on understanding the mathematical aspects of -different algorithms used in the fields of data analysis and machine learning. +Our emphasis throughout this series of lectures is on understanding +the mathematical aspects of different algorithms used in the fields of +data analysis and machine learning. -However, where possible we will emphasize the -importance of using available software. We start thus with a hands-on -and top-down approach to machine learning. The aim is thus to start with -relevant data or data we have produced -and use these to introduce statistical data analysis -concepts and machine learning algorithms before we delve into the -algorithms themselves. The examples we will use in the beginning, start with simple -polynomials with random noise added. We will use the Python -software package "Scikit-Learn":"http://scikit-learn.org/stable/" and -introduce various machine learning algorithms to make fits of -the data and predictions. We move thereafter to more interesting -cases such as data from say experiments (below we will look at experimental nuclear binding energies as an example). -These are examples where we can easily set up the data and -then use machine learning algorithms included in for example -_Scikit-Learn_. +However, where possible we will emphasize the importance of using +available software. We start thus with a hands-on and top-down +approach to machine learning. The aim is thus to start with relevant +data or data we have produced and use these to introduce statistical +data analysis concepts and machine learning algorithms before we delve +into the algorithms themselves. The examples we will use in the +beginning, start with simple polynomials with random noise added. We +will use the Python software package +"Scikit-Learn":"http://scikit-learn.org/stable/" and introduce various +machine learning algorithms to make fits of the data and +predictions. We move thereafter to more interesting cases such as data +from say experiments (below we will look at experimental nuclear +binding energies as an example). These are examples where we can +easily set up the data and then use machine learning algorithms +included in for example _Scikit-Learn_. These examples will serve us the purpose of getting started. Furthermore, they allow us to catch more than two birds with a stone. They will allow us to bring in some programming specific -topics and tools as well as showing the power of various Python -libraries for machine learning and statistical data analysis. +topics and tools as well as showing the power of various Python +libraries for machine learning and statistical data analysis. -Here, we will mainly focus on two -specific Python packages for Machine Learning, Scikit-Learn and -Tensorflow (see below for links etc). Moreover, the examples we -introduce will serve as inputs to many of our discussions later, as -well as allowing you to set up models and produce your own data and -get started with programming. +Here, we will mainly focus on two specific Python packages for Machine +Learning, Scikit-Learn and Tensorflow (see below for links etc). +Moreover, the examples we introduce will serve as inputs to many of +our discussions later, as well as allowing you to set up models and +produce your own data and get started with programming. @@ -275,7 +274,17 @@ quality figures. Feel free to explore the extensive "gallery":"https://matplotlib.org/gallery/index.html" of examples. In this example we plot our original values of $x$ and $y$ as well as the prediction _ypredict_ ($\tilde{y}$), which attempts at fitting our -data with a straight line. +data with a straight line. Note also that _Scikit-Learn_ requires a +matrix as input for the input values $x$ and $y$. In the above code we +have solved this by declaring $x$ and $y$ as arrays of dimension +$n\times 1$. + +In the code here we have also made a new array for $x\in [0,1]$. Our +prediction is computed for these values, meaning that they were not +included in the data set used to *train* (or fit) the model. +This is a recurrring theme in machine learning and data analysis. We would like to train a model on a specific given data set. +Thereafter we wish to apply it to data which were not included in the training. Below we will encounter this again in the so-called *train-validate-test* spliting. We will typically split our data into different sets, oen for training, one for validation and finally, our data from the untouched test vault! + The Python code follows here. !bc pycod @@ -288,6 +297,7 @@ x = np.random.rand(100,1) y = 2*x+np.random.randn(100,1) linreg = LinearRegression() linreg.fit(x,y) +# This is our new x-array to which we test our model xnew = np.array([[0],[1]]) ypredict = linreg.predict(xnew) @@ -394,8 +404,7 @@ plt.show() Depending on the parameter in front of the normal distribution, we may have a small or larger relative error. Try to play around with different training data sets and study (graphically) the value of the -relative error. Note also that _Scikit-Learn_ requires a matrix as input for the input values $x$ and $y$. In the above code we have -solved this by declaring $x$ and $y$ as arrays of dimension $n\times 1$. +relative error. As mentioned above, _Scikit-Learn_ has an impressive functionality. We can for example extract the values of $\alpha$ and $\beta$ and @@ -756,7 +765,7 @@ Energies = Masses['Ebinding'] print(Masses) !ec The next step, and we will define this mathematically later, is to set up the so-called _design matrix_. We will throughout call this matrix $\bm{X}$. -It has dimensionality $p\times n$, where $n$ is the number of data points and $p$ are the so-called predictors. In our case here they are given by the number of polynomials in $A$ we wish to include in the fit. +It has dimensionality $n\times p$, where $n$ is the number of data points and $p$ are the so-called predictors. In our case here they are given by the number of polynomials in $A$ we wish to include in the fit. !bc pycod # Now we set up the design matrix X X = np.zeros((len(A),5)) @@ -767,9 +776,12 @@ X[:,3] = A**(-1.0/3.0) X[:,4] = A**(-1.0) !ec -Note well that we have made life simple here. We perform a fit in terms of the number of nucleons only. A more sophisticated fit can be done by including an explicit dependence on the number of protons and neutrons in the asymmetry and Coulomb terms. +Note well that we have made life simple here. We perform a fit in +terms of the number of nucleons only. A more sophisticated fit can be +done by including an explicit dependence on the number of protons and +neutrons in the asymmetry and Coulomb terms. We leave this as an exercise to you the reader. -With _scikitlearn_ we are now ready to use linear regression and fit our data. +With _Scikit-Learn_ we are now ready to use linear regression and fit our data. !bc pycod clf = skl.LinearRegression().fit(X, Energies) fity = clf.predict(X) @@ -799,7 +811,7 @@ plt.show() !ec -As a teaser, let us now see how we can do this with decision trees using _scikit-learn_. Later we will switch to so-called _random forests_! +As a teaser, let us now see how we can do this with decision trees using _Scikit-Learn_. Later we will switch to so-called _random forests_! !bc pycod @@ -835,6 +847,16 @@ print(Masses) print(np.mean( (Energies-y_1)**2)) !ec +With a deeper and deeper tree level, we can almost reproduce every +single data point by increasing the max depth of the tree. +We can actually decide to make a decision tree which fits every single point. +As we will +see later, this has the benefit that we can really train a model which +traverses every single data point. However, the price we pay is that +we will easily overfit. That is, if we apply our model to unseen data, +we will most likely fail miserably in our attempt at making +predictions. As an exercise, try to make the tree level larger by adjusting the maximum depth variable. When printing out the predicition, you will note that the binding energy of every nucleus is accurately reproduced. + The _seaborn_ package allows us to visualize data in an efficient way. Note that we use _scikit-learn_'s multi-layer perceptron (or feed forward neural network) functionality. @@ -1270,12 +1292,12 @@ matrices as upper case boldfaced letters. !bt \[ -\frac{\partial (\bm{b}^T\bm{a})}{\partial\bm{a}}=\bm{b}, +\frac{\partial\bm{b}^T\bm{a}}{\partial\bm{a}}=\bm{b}, \] !et !bt \[ -\frac{\partial (\bm{a}^T\bm{A}\bm{a})}{\partial\bm{a}}=(\bm{A}+\bm{A}^T)\bm{a}, +\frac{\partial\bm{a}^T\bm{A}\bm{a}}{\partial\bm{a}}=(\bm{A}+\bm{A}^T)\bm{a}, \] !et !bt @@ -1289,6 +1311,30 @@ matrices as upper case boldfaced letters. \] !et +We can then compute the second derivative of the cost function, which in our case is the second derivative +of the means squared error. This leads to + +!bt +\[ +\frac{\partial^2 C(\bm{\beta})}{\partial \bm{\beta}^T\partial \bm{\beta}} =\frac{2}{n}\bm{X}^T\bm{X}. +\] +!et +This quantity defines was what is called the Hessian matrix (the second derivative of a function we want to optimize). + +The Hessian matrix plays an important role and is defined for the mean squared error as + +!bt +\[ +\bm{H}=\bm{X}^T\bm{X}. +\] +!et + + +The Hessian matrix for ordinary least squares is also proportional to +the covariance matrix. As we will see in the chapter on Ridge and Lasso regression, This means that we can use the Singular Value Decomposition of a matrix to find +the eigenvalues of the covariance matrix and the Hessian matrix in +terms of the singular values. + The residuals $\bm{\epsilon}$ are in turn given by !bt @@ -1546,9 +1592,8 @@ before, with the same initializations and declarations. We use also _pandas_ again, rather extensively in order to organize our data. The difference now is that we use _Scikit-Learn's_ regression tools -instead of our own matrix inversion implementation. Furthermore, we -sneak in _Ridge_ regression (to be discussed below) which includes a -hyperparameter $\lambda$, also to be explained below. +instead of our own matrix inversion implementation. + !bc pycod @@ -1612,18 +1657,6 @@ print('Variance score: %.2f' % r2_score(Energies, ytilde)) print('Mean absolute error: %.2f' % mean_absolute_error(Energies, ytilde)) print(clf.coef_, clf.intercept_) -# The Ridge regression with a hyperparameter lambda = 0.1 -_lambda = 0.1 -clf_ridge = skl.Ridge(alpha=_lambda).fit(X, Energies) -yridge = clf_ridge.predict(X) -EoS['Eridge'] = yridge -# The mean squared error -print("Mean squared error: %.2f" % mean_squared_error(Energies, yridge)) -# Explained variance score: 1 is perfect prediction -print('Variance score: %.2f' % r2_score(Energies, yridge)) -# Mean absolute error -print('Mean absolute error: %.2f' % mean_absolute_error(Energies, yridge)) -print(clf_ridge.coef_, clf_ridge.intercept_) fig, ax = plt.subplots() ax.set_xlabel(r'$\rho[\mathrm{fm}^{-3}]$') @@ -1632,8 +1665,6 @@ ax.plot(EoS['Density'], EoS['Energy'], alpha=0.7, lw=2, label='Theoretical data') ax.plot(EoS['Density'], EoS['Eols'], alpha=0.7, lw=2, c='m', label='OLS') -ax.plot(EoS['Density'], EoS['Eridge'], alpha=0.7, lw=2, c='g', - label='Ridge $\lambda = 0.1$') ax.legend() save_fig("EoSfitting") plt.show() @@ -1642,13 +1673,11 @@ plt.show() The above simple polynomial in density $\rho$ gives an excellent fit to the data. -We note also that there is a small deviation between the -standard OLS and the Ridge regression at higher densities. We discuss this in more detail -below. ===== Splitting our Data in Training and Test data ===== + It is normal in essentially all Machine Learning studies to split the data in a training set and a test set (sometimes also an additional validation set). _Scikit-Learn_ has an own function for this. There @@ -1660,6 +1689,81 @@ our discussion of the so-called _bias-variance_ tradeoff. Here we limit ourselves to repeat the above equation of state fitting example but now splitting the data into a training set and a test set. +Let us study some examples. The first code here takes a simple +one-dimensional second-order polynomial and we fit it to a +second-order polynomial. Depending on the strength of the added noise, +the various measures like the $R2$ score or the mean-squared error, +the fit becomes better or worse. + +!bc pycod +import os +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from sklearn.model_selection import train_test_split + + +def R2(y_data, y_model): + return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) +def MSE(y_data,y_model): + n = np.size(y_model) + return np.sum((y_data-y_model)**2)/n + +x = np.random.rand(100) +y = 2.0+5*x*x+0.1*np.random.randn(100) + + +# The design matrix now as function of a given polynomial +X = np.zeros((len(x),3)) +X[:,0] = 1.0 +X[:,1] = x +X[:,2] = x**2 +# 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) +# matrix inversion to find beta +beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train +print(beta) +# and then make the prediction +ytilde = X_train @ beta +print("Training R2") +print(R2(y_train,ytilde)) +print("Training MSE") +print(MSE(y_train,ytilde)) +ypredict = X_test @ beta +print("Test R2") +print(R2(y_test,ypredict)) +print("Test MSE") +print(MSE(y_test,ypredict)) +!ec + +Alternatively, you could write your own test-train splitting function as shown here. + +!bc pycod +# equivalently in numpy +def train_test_split_numpy(inputs, labels, train_size, test_size): + n_inputs = len(inputs) + inputs_shuffled = inputs.copy() + labels_shuffled = labels.copy() + + np.random.shuffle(inputs_shuffled) + np.random.shuffle(labels_shuffled) + + train_end = int(n_inputs*train_size) + X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:] + Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:] + + return X_train, X_test, Y_train, Y_test +!ec + +But since _scikit-learn_ has its own function for doing this and since +it interfaces easily with _tensorflow_ and other libraries, we +normally recommend using the latter functionality. + + +As another example, we apply the training and testing split to +to the above equation of state fitting example +but now splitting the data into a training set and a test set. + !bc pycod import os import numpy as np @@ -1881,91 +1985,6 @@ plt.show() !ec -===== Splitting our Data in Training and Test data ===== - - -It is normal in essentially all Machine Learning studies to split the -data in a training set and a test set (sometimes also an additional -validation set). _Scikit-Learn_ has an own function for this. There -is no explicit recipe for how much data should be included as training -data and say test data. An accepted rule of thumb is to use -approximately $2/3$ to $4/5$ of the data as training data. We will -postpone a discussion of this splitting to the end of these notes and -our discussion of the so-called _bias-variance_ tradeoff. Here we -limit ourselves to repeat the above equation of state fitting example -but now splitting the data into a training set and a test set. - -Let us study some examples. The first code here takes a simple -one-dimensional second-order polynomial and we fit it to a -second-order polynomial. Depending on the strength of the added noise, -the various measures like the $R2$ score or the mean-squared error, -the fit becomes better or worse. - -!bc pycod -import os -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -from sklearn.model_selection import train_test_split - - -def R2(y_data, y_model): - return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2) -def MSE(y_data,y_model): - n = np.size(y_model) - return np.sum((y_data-y_model)**2)/n - -x = np.random.rand(100) -y = 2.0+5*x*x+0.1*np.random.randn(100) - - -# The design matrix now as function of a given polynomial -X = np.zeros((len(x),3)) -X[:,0] = 1.0 -X[:,1] = x -X[:,2] = x**2 -# 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) -# matrix inversion to find beta -beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train -print(beta) -# and then make the prediction -ytilde = X_train @ beta -print("Training R2") -print(R2(y_train,ytilde)) -print("Training MSE") -print(MSE(y_train,ytilde)) -ypredict = X_test @ beta -print("Test R2") -print(R2(y_test,ypredict)) -print("Test MSE") -print(MSE(y_test,ypredict)) -!ec - -Alternatively, you could write your own test-train splitting function as shown here. - -!bc pycod -# equivalently in numpy -def train_test_split_numpy(inputs, labels, train_size, test_size): - n_inputs = len(inputs) - inputs_shuffled = inputs.copy() - labels_shuffled = labels.copy() - - np.random.shuffle(inputs_shuffled) - np.random.shuffle(labels_shuffled) - - train_end = int(n_inputs*train_size) - X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:] - Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:] - - return X_train, X_test, Y_train, Y_test -!ec - -But since _scikit-learn_ has its own function for doing this and since -it interfaces easily with _tensorflow_ and other libraries, we -normally recommend using the latter functionality. - - ===== Reducing the number of degrees of freedom, overarching view ===== diff --git a/doc/BookChapters/chapter2.do.txt b/doc/BookChapters/chapter2.do.txt index dcb02629f..4d02b5560 100644 --- a/doc/BookChapters/chapter2.do.txt +++ b/doc/BookChapters/chapter2.do.txt @@ -7,11 +7,11 @@ What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD). -We have shown that in ordinary least squares the optimal parameters $\beta$ are given by +We have shown that in ordinary least squares (OLS) the optimal parameters $\beta$ are given by !bt \[ -\hat{\bm{\beta}} = \left(\bm{X}^T\bm{X}\right)^{-1}\bm{X}^T\bm{y}. +\hat{\bm{\beta}}_{\mathrm{OLS}} = \left(\bm{X}^T\bm{X}\right)^{-1}\bm{X}^T\bm{y}. \] !et @@ -39,7 +39,7 @@ We can rewrite \] !et -The matrix $\bm{A}$ has the important property that $\bm{A}^2=\bm{A}$. This is the definition of a projection matrix. +The matrix $\bm{A}$ has the important property that $\bm{A}^2=\bm{A}$. This is the definition of a "projection matrix":"https://en.wikipedia.org/wiki/Projection_matrix". We can then interpret our optimal model $\tilde{\bm{y}}$ as being represented by an orthogonal projection of $\bm{y}$ onto a space defined by the column vectors of $\bm{X}$. In our case here the matrix $\bm{A}$ is a square matrix. If it is a general rectangular matrix we have an oblique projection matrix. @@ -352,7 +352,7 @@ example ===== Code for SVD and Inversion of Matrices ===== -How do we use the SVD to invert a matrix $\bm{X}^\bm{X}$ which is singular or near singular? +How do we use the SVD to invert a matrix $\bm{X}^T\bm{X}$ which is singular or near singular? The simple answer is to use the linear algebra function for the pseudoinverse, that is !bc pycod #Ainv = np.linlag.pinv(A) @@ -2057,7 +2057,7 @@ y = matrix( [4, 2, 3]) x = l1regls(X,y) !ec - +_More text will be added to this example._ ===== Linking the regression analysis with a statistical interpretation ===== @@ -2279,7 +2279,9 @@ p(\bm{D}\vert\bm{\beta})=\prod_{i=0}^{n-1}\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\lef \] !et -It is a conditional probability (see below) and reads as the likelihood of a domain of events $\bm{D}$ given a set of parameters $\bm{\beta}$. +It is a conditional probability (see below) and reads as the +likelihood of a domain of events $\bm{D}$ given a set of parameters +$\bm{\beta}$. In statistics, maximum likelihood estimation (MLE) is a method of @@ -2348,47 +2350,40 @@ Assume we have two domains of events $X=[x_0,x_1,\dots,x_{n-1}]$ and $Y=[y_0,y_1 We define also the likelihood for $X$ and $Y$ as $p(X)$ and $p(Y)$ respectively. The likelihood of a specific event $x_i$ (or $y_i$) is then written as $p(X=x_i)$ or just $p(x_i)=p_i$. -!bblock Union of events is given by +The union of events is given by !bt \[ p(X \cup Y)= p(X)+p(Y)-p(X \cap Y). \] !et -!eblock -!bblock The product rule (aka joint probability) is given by +The product rule (aka joint probability) is given by !bt \[ p(X \cup Y)= p(X,Y)= p(X\vert Y)p(Y)=p(Y\vert X)p(X), \] !et where we read $p(X\vert Y)$ as the likelihood of obtaining $X$ given $Y$. -!eblock If we have independent events then $p(X,Y)=p(X)p(Y)$. The marginal probability is defined in terms of only one of the set of variables $X,Y$. For a discrete probability we have -!bblock !bt \[ p(X)=\sum_{i=0}^{n-1}p(X,Y=y_i)=\sum_{i=0}^{n-1}p(X\vert Y=y_i)p(Y=y_i)=\sum_{i=0}^{n-1}p(X\vert y_i)p(y_i). \] !et -!eblock - The conditional probability, if $p(Y) > 0$, is -!bblock !bt \[ p(X\vert Y)= \frac{p(X,Y)}{p(Y)}=\frac{p(X,Y)}{\sum_{i=0}^{n-1}p(Y\vert X=x_i)p(x_i)}. \] !et -!eblock @@ -2413,21 +2408,27 @@ evaluated for the observed data $Y$ and can be viewed as a function of the parameter space represented by $X$. This function is not necesseraly normalized and is normally called the likelihood function. -The function $p(X)$ on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution. +The function $p(X)$ on the right hand side is called the prior while +the function on the left hand side is the called the posterior +probability. The denominator on the right hand side serves as a +normalization factor for the posterior distribution. Let us try to illustrate Bayes' theorem through an example. -Let us suppose that you are undergoing a series of mammography scans in -order to rule out possible breast cancer cases. We define the +Let us suppose that you are undergoing a series of mammography scans +in order to rule out possible breast cancer cases. We define the sensitivity for a positive event by the variable $X$. It takes binary values with $X=1$ representing a positive event and $X=0$ being a negative event. We reserve $Y$ as a classification parameter for -either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing). +either a negative or a positive breast cancer confirmation. (Short +note on wordings: positive here means having breast cancer, although +none of us would consider this being a positive thing). We let $Y=1$ represent the the case of having breast cancer and $Y=0$ as not. -Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$, that is we have +Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$ (the numbers here are all made up), +that is we have !bt \[ @@ -2435,12 +2436,22 @@ p(X=1\vert Y=1) =0.8. \] !et -This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of $80\%$ for having cancer. -It is however not correct, as the following Bayesian analysis shows. +This obviously sounds scary since many would conclude that if the test +is positive, there is a likelihood of $80\%$ for having cancer. It is +however not correct, as the following Bayesian analysis shows. The correct question to pose is *what is the probability of having breast cancer in case of a positive test?* +We are thus interested in +!bt +\[ +p(Y=1\vert X=1), +\] +!et +instead of $p(X=1\vert Y=1)$. -If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number. -Let us assume that the prior probability in the population as a whole is + +If we look at various national surveys on breast cancer, the general +likelihood of developing breast cancer is a very small number. Let us +assume that the prior probability in the population as a whole is !bt \[ @@ -2448,19 +2459,24 @@ p(Y=1) =0.004. \] !et -We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have +We need also to account for the fact that the test may produce a false +positive result (false alarm). Let us here assume that we have + !bt \[ p(X=1\vert Y=0) =0.1. \] !et -Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute +Using Bayes' theorem we can then find the posterior probability that +the person has breast cancer in case of a positive test, that is we +can compute !bt -\[ -p(Y=1\vert X=1)=\frac{p(X=1\vert Y=1)p(Y=1)}{p(X=1\vert Y=1)p(Y=1)+p(X=1\vert Y=0)p(Y=0)}=\frac{0.8\times 0.004}{0.8\times 0.004+0.1\times 0.996}=0.031. -\] +\begin{align} +p(Y=1\vert X=1)=\frac{p(X=1\vert Y=1)p(Y=1)}{p(X=1\vert Y=1)p(Y=1)+p(X=1\vert Y=0)p(Y=0)}&= \\ +& \frac{0.8\times 0.004}{0.8\times 0.004+0.1\times 0.996}=0.031. +\end{align} !et That is, in case of a positive test, there is only a $3\%$ chance of having breast cancer! @@ -2480,9 +2496,12 @@ We will play around with a study of the values for the optimal parameters $\bm{\beta}$ using OLS, Ridge and Lasso regression. For OLS, you will notice as function of the noise and polynomial degree, that the parameters $\beta$ will fluctuate from order to order in the -polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS. +polynomial fit and that for larger and larger polynomial degrees of +freedom, the parameters will tend to increase in value for OLS. -For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one. +For Ridge and Lasso regression, the higher order parameters will +typically be reduced, providing thereby less fluctuations from one +order to another one. !bc pycod import numpy as np @@ -2538,8 +2557,8 @@ for i in range(nlambdas): # Compute the MSE and print it MSERidgePredict[i] = MSE(y_test,ypredictRidge) MSELassoPredict[i] = MSE(y_test,ypredictLasso) - print(lmb,RegRidge.coef_) - print(lmb,RegLasso.coef_) +# print(lmb,RegRidge.coef_) +# print(lmb,RegLasso.coef_) # Now plot the results plt.figure() plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test') @@ -2552,10 +2571,13 @@ plt.show() !ec -How can we understand this? +How can we understand this? _More text to be added_. +===== Linking Bayes' Theorem with Ridge and Lasso Regression ===== + + Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression. For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\bm{D}$ (one-dimensional case) diff --git a/doc/BookChapters/chapter3.do.txt b/doc/BookChapters/chapter3.do.txt index 16031ac61..45ea40e7d 100644 --- a/doc/BookChapters/chapter3.do.txt +++ b/doc/BookChapters/chapter3.do.txt @@ -427,160 +427,6 @@ plt.show() !ec - -===== Various steps in cross-validation ===== - -When the repetitive splitting of the data set is done randomly, -samples may accidently end up in a fast majority of the splits in -either training or test set. Such samples may have an unbalanced -influence on either model building or prediction evaluation. To avoid -this $k$-fold cross-validation structures the data splitting. The -samples are divided into $k$ more or less equally sized exhaustive and -mutually exclusive subsets. In turn (at each split) one of these -subsets plays the role of the test set while the union of the -remaining subsets constitutes the training set. Such a splitting -warrants a balanced representation of each sample in both training and -test set over the splits. Still the division into the $k$ subsets -involves a degree of randomness. This may be fully excluded when -choosing $k=n$. This particular case is referred to as leave-one-out -cross-validation (LOOCV). - - -* Define a range of interest for the penalty parameter. - -* Divide the data set into training and test set comprising samples $\{1, \ldots, n\} \setminus i$ and $\{ i \}$, respectively. - -* Fit the linear regression model by means of ridge estimation for each $\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\bm{\sigma}_{-i}^2(\lambda)$, as -!bt -\begin{align*} -\bm{\beta}_{-i}(\lambda) & = ( \bm{X}_{-i, \ast}^{T} -\bm{X}_{-i, \ast} + \lambda \bm{I}_{pp})^{-1} -\bm{X}_{-i, \ast}^{T} \bm{y}_{-i} -\end{align*} -!et - -* Evaluate the prediction performance of these models on the test set by $\log\{L[y_i, \bm{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}$. Or, by the prediction error $|y_i - \bm{X}_{i, \ast} \bm{\beta}_{-i}(\lambda)|$, the relative error, the error squared or the R2 score function. - -* Repeat the first three steps such that each sample plays the role of the test set once. - -* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as -!bt -\begin{align*} -\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}. -\end{align*} -!et - - -For the various values of $k$ - -o shuffle the dataset randomly. -o Split the dataset into $k$ groups. -o For each unique group: - o Decide which group to use as set for test data - o Take the remaining groups as a training data set - o Fit a model on the training set and evaluate it on the test set - o Retain the evaluation score and discard the model -o Summarize the model using the sample of model evaluation scores - - - - -The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial. -!bc pycod -import numpy as np -import matplotlib.pyplot as plt -from sklearn.model_selection import KFold -from sklearn.linear_model import Ridge -from sklearn.model_selection import cross_val_score -from sklearn.preprocessing import PolynomialFeatures - -# A seed just to ensure that the random numbers are the same for every run. -# Useful for eventual debugging. -np.random.seed(3155) - -# Generate the data. -nsamples = 100 -x = np.random.randn(nsamples) -y = 3*x**2 + np.random.randn(nsamples) - -## Cross-validation on Ridge regression using KFold only - -# Decide degree on polynomial to fit -poly = PolynomialFeatures(degree = 6) - -# Decide which values of lambda to use -nlambdas = 500 -lambdas = np.logspace(-3, 5, nlambdas) - -# Initialize a KFold instance -k = 5 -kfold = KFold(n_splits = k) - -# Perform the cross-validation to estimate MSE -scores_KFold = np.zeros((nlambdas, k)) - -i = 0 -for lmb in lambdas: - ridge = Ridge(alpha = lmb) - j = 0 - for train_inds, test_inds in kfold.split(x): - xtrain = x[train_inds] - ytrain = y[train_inds] - - xtest = x[test_inds] - ytest = y[test_inds] - - Xtrain = poly.fit_transform(xtrain[:, np.newaxis]) - ridge.fit(Xtrain, ytrain[:, np.newaxis]) - - Xtest = poly.fit_transform(xtest[:, np.newaxis]) - ypred = ridge.predict(Xtest) - - scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred) - - j += 1 - i += 1 - - -estimated_mse_KFold = np.mean(scores_KFold, axis = 1) - -## Cross-validation using cross_val_score from sklearn along with KFold - -# kfold is an instance initialized above as: -# kfold = KFold(n_splits = k) - -estimated_mse_sklearn = np.zeros(nlambdas) -i = 0 -for lmb in lambdas: - ridge = Ridge(alpha = lmb) - - X = poly.fit_transform(x[:, np.newaxis]) - estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold) - - # cross_val_score return an array containing the estimated negative mse for every fold. - # we have to the the mean of every array in order to get an estimate of the mse of the model - estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds) - - i += 1 - -## Plot and compare the slightly different ways to perform cross-validation - -plt.figure() - -plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score') -plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold') - -plt.xlabel('log10(lambda)') -plt.ylabel('mse') - -plt.legend() - -plt.show() - -!ec - - - ===== The bias-variance tradeoff ===== @@ -957,6 +803,163 @@ plt.show() + + + +===== Cross-validation ===== + +When the repetitive splitting of the data set is done randomly, +samples may accidently end up in a fast majority of the splits in +either training or test set. Such samples may have an unbalanced +influence on either model building or prediction evaluation. To avoid +this $k$-fold cross-validation structures the data splitting. The +samples are divided into $k$ more or less equally sized exhaustive and +mutually exclusive subsets. In turn (at each split) one of these +subsets plays the role of the test set while the union of the +remaining subsets constitutes the training set. Such a splitting +warrants a balanced representation of each sample in both training and +test set over the splits. Still the division into the $k$ subsets +involves a degree of randomness. This may be fully excluded when +choosing $k=n$. This particular case is referred to as leave-one-out +cross-validation (LOOCV). + + +* Define a range of interest for the penalty parameter. + +* Divide the data set into training and test set comprising samples $\{1, \ldots, n\} \setminus i$ and $\{ i \}$, respectively. + +* Fit the linear regression model by means of ridge estimation for each $\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\bm{\sigma}_{-i}^2(\lambda)$, as +!bt +\begin{align*} +\bm{\beta}_{-i}(\lambda) & = ( \bm{X}_{-i, \ast}^{T} +\bm{X}_{-i, \ast} + \lambda \bm{I}_{pp})^{-1} +\bm{X}_{-i, \ast}^{T} \bm{y}_{-i} +\end{align*} +!et + +* Evaluate the prediction performance of these models on the test set by $\log\{L[y_i, \bm{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}$. Or, by the prediction error $|y_i - \bm{X}_{i, \ast} \bm{\beta}_{-i}(\lambda)|$, the relative error, the error squared or the R2 score function. + +* Repeat the first three steps such that each sample plays the role of the test set once. + +* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as +!bt +\begin{align*} +\frac{1}{n} \sum_{i = 1}^n \log\{L[y_i, \mathbf{X}_{i, \ast}; \bm{\beta}_{-i}(\lambda), \bm{\sigma}_{-i}^2(\lambda)]\}. +\end{align*} +!et + + +For the various values of $k$ + +o shuffle the dataset randomly. +o Split the dataset into $k$ groups. +o For each unique group: + o Decide which group to use as set for test data + o Take the remaining groups as a training data set + o Fit a model on the training set and evaluate it on the test set + o Retain the evaluation score and discard the model +o Summarize the model using the sample of model evaluation scores + + + + +The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial. +!bc pycod +import numpy as np +import matplotlib.pyplot as plt +from sklearn.model_selection import KFold +from sklearn.linear_model import Ridge +from sklearn.model_selection import cross_val_score +from sklearn.preprocessing import PolynomialFeatures + +# A seed just to ensure that the random numbers are the same for every run. +# Useful for eventual debugging. +np.random.seed(3155) + +# Generate the data. +nsamples = 100 +x = np.random.randn(nsamples) +y = 3*x**2 + np.random.randn(nsamples) + +## Cross-validation on Ridge regression using KFold only + +# Decide degree on polynomial to fit +poly = PolynomialFeatures(degree = 6) + +# Decide which values of lambda to use +nlambdas = 500 +lambdas = np.logspace(-3, 5, nlambdas) + +# Initialize a KFold instance +k = 5 +kfold = KFold(n_splits = k) + +# Perform the cross-validation to estimate MSE +scores_KFold = np.zeros((nlambdas, k)) + +i = 0 +for lmb in lambdas: + ridge = Ridge(alpha = lmb) + j = 0 + for train_inds, test_inds in kfold.split(x): + xtrain = x[train_inds] + ytrain = y[train_inds] + + xtest = x[test_inds] + ytest = y[test_inds] + + Xtrain = poly.fit_transform(xtrain[:, np.newaxis]) + ridge.fit(Xtrain, ytrain[:, np.newaxis]) + + Xtest = poly.fit_transform(xtest[:, np.newaxis]) + ypred = ridge.predict(Xtest) + + scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred) + + j += 1 + i += 1 + + +estimated_mse_KFold = np.mean(scores_KFold, axis = 1) + +## Cross-validation using cross_val_score from sklearn along with KFold + +# kfold is an instance initialized above as: +# kfold = KFold(n_splits = k) + +estimated_mse_sklearn = np.zeros(nlambdas) +i = 0 +for lmb in lambdas: + ridge = Ridge(alpha = lmb) + + X = poly.fit_transform(x[:, np.newaxis]) + estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold) + + # cross_val_score return an array containing the estimated negative mse for every fold. + # we have to the the mean of every array in order to get an estimate of the mse of the model + estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds) + + i += 1 + +## Plot and compare the slightly different ways to perform cross-validation + +plt.figure() + +plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score') +plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold') + +plt.xlabel('log10(lambda)') +plt.ylabel('mse') + +plt.legend() + +plt.show() + +!ec + + +More examples of the application of cross-validation follow here. + !bc pycod # Common imports import os diff --git a/doc/BookChapters/chapter4.do.txt b/doc/BookChapters/chapter4.do.txt index c7a2771c6..94dbf7b95 100644 --- a/doc/BookChapters/chapter4.do.txt +++ b/doc/BookChapters/chapter4.do.txt @@ -1,8 +1,6 @@ ======= Logistic Regression ======= -"Video of Lecture":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK3155/h20/forelesningsvideoer/LectureSeptember18.mp4?vrtx=view-as-webpage" - ===== Logistic Regression ===== @@ -622,7 +620,7 @@ plt.show() -===== Optimization, the central part of any Machine Learning algortithm ===== +======= Optimization, the central part of any Machine Learning algortithm ======= Almost every problem in machine learning and data science starts with a dataset $X$, a model $g(\beta)$, which is a function of the @@ -635,7 +633,6 @@ some approximative/numerical method to compute the minimum. -===== Revisiting our Logistic Regression case ===== In our discussion on Logistic Regression we studied the case of @@ -653,7 +650,6 @@ p(y_i=0|x_i,\bm{\beta}) &= 1 - p(y_i=1|x_i,\bm{\beta}), where $\bm{\beta}$ are the weights we wish to extract from data, in our case $\beta_0$ and $\beta_1$. -===== The equations to solve ===== Our compact equations used a definition of a vector $\bm{y}$ with $n$ elements $y_i$, an $n\times p$ matrix $\bm{X}$ which contains the @@ -678,7 +674,6 @@ $p(y_i\vert x_i,\bm{\beta})(1-p(y_i\vert x_i,\bm{\beta})$, we can obtain a compa This defines what is called the Hessian matrix. -===== Solving using Newton-Raphson's method ===== If we can set up these equations, Newton-Raphson's iterative method is normally the method of choice. It requires however that we can compute in an efficient way the matrices that define the first and second derivatives. @@ -701,9 +696,6 @@ The right-hand side is computed with the old values of $\beta$. If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement. - -===== Brief reminder on Newton-Raphson's method ===== - Let us quickly remind ourselves how we derive the above method. Perhaps the most celebrated of all one-dimensional root-finding @@ -715,7 +707,6 @@ numerically and/or your function is not of the smooth type, we normally discourage the use of this method. -===== The equations ===== The Newton-Raphson formula consists geometrically of extending the tangent line at a current point until it crosses zero, then setting @@ -755,7 +746,6 @@ Having in mind an iterative procedure, it is natural to start iterating with !et -===== Simple geometric interpretation ===== The above is Newton-Raphson's method. It has a simple geometric interpretation, namely $x_{n+1}$ is the point where the tangent from @@ -771,7 +761,6 @@ vanishes, then Newton-Raphson may fail totally -===== Extending to more than one variable ===== Newton's method can be generalized to systems of several non-linear equations and variables. Consider the case with two equations @@ -848,8 +837,6 @@ For $\gamma_k$ small enough, then $F(\mathbf{x}_{k+1}) \leq F(\mathbf{x}_k)$. This means that for a sufficiently small $\gamma_k$ we are always moving towards smaller function values, i.e a minimum. - -===== More on Steepest descent ===== The previous observation is the basis of the method of steepest descent, which is also referred to as just gradient descent (GD). One @@ -865,8 +852,6 @@ computes new approximations according to The parameter $\gamma_k$ is often referred to as the step length or the learning rate within the context of Machine Learning. - -===== The ideal ===== Ideally the sequence $\{\mathbf{x}_k \}_{k=0}$ converges to a global minimum of the function $F$. In general we do not know if we are in a @@ -888,7 +873,6 @@ Note that the gradient is a function of $\mathbf{x} = -===== The sensitiveness of the gradient descent ===== The gradient descent method is sensitive to the choice of learning rate $\gamma_k$. This is due @@ -919,13 +903,10 @@ $\mathbb{R}$. Examples of convex sets of $\mathbb{R}^2$ are the regular polygons (triangles, rectangles, pentagons, etc...). -===== Convex function ===== _Convex function_: Let $X \subset \mathbb{R}^n$ be a convex set. Assume that the function $f: X \rightarrow \mathbb{R}$ is continuous, then $f$ is said to be convex if $$f(tx_1 + (1-t)x_2) \leq tf(x_1) + (1-t)f(x_2) $$ for all $x_1, x_2 \in X$ and for all $t \in [0,1]$. If $\leq$ is replaced with a strict inequaltiy in the definition, we demand $x_1 \neq x_2$ and $t\in(0,1)$ then $f$ is said to be strictly convex. For a single variable function, convexity means that if you draw a straight line connecting $f(x_1)$ and $f(x_2)$, the value of the function on the interval $[x_1,x_2]$ is always below the line as illustrated below. -===== Conditions on convex functions ===== - In the following we state first and second-order conditions which ensures convexity of a function $f$. We write $D_f$ to denote the domain of $f$, i.e the subset of $R^n$ where $f$ is defined. For more @@ -954,7 +935,6 @@ everywhere. This condition is particularly useful since it gives us an procedure for determining if the function under consideration is convex, apart from using the definition. -===== More on convex functions ===== The next result is of great importance to us and the reason why we are going on about convex functions. In machine learning we frequently @@ -975,7 +955,7 @@ $x^*$ that satisfies $\nabla f(x^*) = 0$ is a global minimum. This result means that if we know that the cost/loss function is convex and we are able to find a minimum, we are guaranteed that it is a global minimum. -===== Some simple problems ===== +=== Some simple problems === o Show that $f(x)=x^2$ is convex for $x \in \mathbb{R}$ using the definition of convexity. Hint: If you re-write the definition, $f$ is convex if the following holds for all $x,y \in D_f$ and any $\lambda \in [0,1]$ $\lambda f(x)+(1-\lambda)f(y)-f(\lambda x + (1-\lambda) y ) \geq 0$. @@ -992,13 +972,6 @@ o A norm is any function that satisfy the following properties Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this). - -===== Friday September 25 ===== - -"Video of Lecture":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureSeptember25.mp4?vrtx=view-as-webpage" and "link to handwritten notes":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/NotesSeptember25.pdf". - - - ===== Standard steepest descent ===== @@ -1029,7 +1002,6 @@ where $\bm{r}$ is the so-called residual or error in the iterative process. When we have found the exact solution, $\bm{r}=0$. -===== Gradient method ===== The residual is zero when we reach the minimum of the quadratic equation !bt @@ -1042,9 +1014,6 @@ with the constraint that the matrix $\bm{A}$ is positive definite and symmetric. This defines also the Hessian and we want it to be positive definite. - -===== Steepest descent method ===== - We denote the initial guess for $\bm{x}$ as $\bm{x}_0$. We can assume without loss of generality that !bt @@ -1061,9 +1030,6 @@ or consider the system instead. - -===== Steepest descent method ===== -!bblock One can show that the solution $\bm{x}$ is also the unique minimizer of the quadratic form !bt \begin{equation*} @@ -1081,11 +1047,7 @@ which equals and $\bm{x}_0=0$ it is equal $-\bm{b}$. -!eblock - -===== Final expressions ===== -!bblock We can compute the residual iteratively as !bt \begin{equation*} @@ -1117,12 +1079,6 @@ leading to the iterative scheme \bm{x}_{k+1}=\bm{x}_k-\alpha_k\bm{r}_{k}, \end{equation*} !et -!eblock - - - - -===== Steepest descent example ===== !bc pycod import numpy as np @@ -1177,7 +1133,6 @@ pt.plot(it_array.T[0], it_array.T[1], "x-") ===== Conjugate gradient method ===== -!bblock In the CG method we define so-called conjugate directions and two vectors $\bm{s}$ and $\bm{t}$ are said to be @@ -1196,11 +1151,7 @@ of our vectors $\bm{x}_i$ obeying the above criterion, namely !et Two vectors are conjugate if they are orthogonal with respect to this inner product. Being conjugate is a symmetric relation: if $\bm{s}$ is conjugate to $\bm{t}$, then $\bm{t}$ is conjugate to $\bm{s}$. -!eblock - -===== Conjugate gradient method ===== -!bblock An example is given by the eigenvectors of the matrix !bt \begin{equation*} @@ -1208,12 +1159,7 @@ An example is given by the eigenvectors of the matrix \end{equation*} !et which is zero unless $i=j$. -!eblock - - -===== Conjugate gradient method ===== -!bblock Assume now that we have a symmetric positive-definite matrix $\bm{A}$ of size $n\times n$. At each iteration $i+1$ we obtain the conjugate direction of a vector !bt @@ -1230,11 +1176,7 @@ $ \bm{A}\bm{x} = \bm{b}$ in this basis, namely \bm{x} = \sum^{n}_{i=1} \alpha_i \bm{p}_i. \end{equation*} !et -!eblock - -===== Conjugate gradient method ===== -!bblock The coefficients are given by !bt \begin{equation*} @@ -1255,11 +1197,6 @@ and we can define the coefficients $\alpha_k$ as \alpha_k = \frac{\bm{p}_k^T \bm{b}}{\bm{p}_k^T \bm{A} \bm{p}_k} \end{equation*} !et -!eblock - - -===== Conjugate gradient method and iterations ===== -!bblock If we choose the conjugate vectors $\bm{p}_k$ carefully, then we may not need all of them to obtain a good approximation to the solution @@ -1282,12 +1219,7 @@ or consider the system \end{equation*} !et instead. -!eblock - - -===== Conjugate gradient method ===== -!bblock One can show that the solution $\bm{x}$ is also the unique minimizer of the quadratic form !bt \begin{equation*} @@ -1306,12 +1238,7 @@ and $\bm{x}_0=0$ it is equal $-\bm{b}$. The other vectors in the basis will be conjugate to the gradient, hence the name conjugate gradient method. -!eblock - - -===== Conjugate gradient method ===== -!bblock Let $\bm{r}_k$ be the residual at the $k$-th step: !bt \begin{equation*} @@ -1330,11 +1257,7 @@ This gives the following expression \bm{p}_{k+1}=\bm{r}_k-\frac{\bm{p}_k^T \bm{A}\bm{r}_k}{\bm{p}_k^T\bm{A}\bm{p}_k} \bm{p}_k. \end{equation*} !et -!eblock - -===== Conjugate gradient method ===== -!bblock We can also compute the residual iteratively as !bt \begin{equation*} @@ -1360,21 +1283,16 @@ which gives \bm{r}_{k+1}=\bm{r}_k-\bm{A}\bm{p}_{k}, \end{equation*} !et -!eblock - - - - -===== Revisiting our first homework ===== +===== Revisiting our Linear Regression Solvers ===== We will use linear regression as a case study for the gradient descent methods. Linear regression is a great test case for the gradient descent methods discussed in the lectures since it has several desirable properties such as: -o An analytical solution (recall homework set 1). +o An analytical solution. o The gradient can be computed analytically. o The cost function is convex which guarantees that gradient descent converges for small enough learning rates @@ -1398,8 +1316,6 @@ such that \] !et - -===== Gradient descent example ===== Let $\mathbf{y} = (y_1,\cdots,y_n)^T$, $\mathbf{\bm{y}} = (\bm{y}_1,\cdots,\bm{y}_n)^T$ and $\beta = (\beta_0, \beta_1)^T$ @@ -1422,8 +1338,6 @@ C(\beta) = \frac{1}{n}||X\beta-\mathbf{y}||_{2}^{2} = \frac{1}{n}\sum_{i=1}^{100 and we want to find $\beta$ such that $C(\beta)$ is minimized. -===== The derivative of the cost/loss function ===== - Computing $\partial C(\beta) / \partial \beta_0$ and $\partial C(\beta) / \partial \beta_1$ we can show that the gradient can be written as !bt \[ @@ -1435,7 +1349,6 @@ Computing $\partial C(\beta) / \partial \beta_0$ and $\partial C(\beta) / \parti where $X$ is the design matrix defined above. -===== The Hessian matrix ===== The Hessian matrix of $C(\beta)$ is given by !bt \[ @@ -1449,10 +1362,6 @@ This result implies that $C(\beta)$ is a convex function since the matrix $X^T X - - -===== Simple program ===== - We can now write a program that minimizes $C(\beta)$ using the gradient descent method with a constant learning rate $\gamma$ according to !bt \[ @@ -1468,9 +1377,8 @@ And finally we can compare our solution for $\beta$ with the analytic result giv $\beta= (X^TX)^{-1} X^T \mathbf{y}$. -===== Gradient Descent Example ===== -Here our simple example +Here is our simple example !bc pycod # Importing various packages @@ -1521,8 +1429,7 @@ plt.show() !ec - -===== And a corresponding example using _scikit-learn_ ===== +Alternatively, we can use _Scikit-Learn_ as done here !bc pycod # Importing various packages @@ -1545,10 +1452,6 @@ print(sgdreg.intercept_, sgdreg.coef_) !ec - - -===== Gradient descent and Ridge ===== - We have also discussed Ridge regression where the loss function contains a regularized term given by the $L_2$ norm of $\beta$, !bt \[ @@ -1573,8 +1476,6 @@ We can easily extend our program to minimize $C_{\text{ridge}}(\beta)$ using gra !et - -===== Program example for gradient descent with Ridge Regression ===== !bc pycod from random import random, seed import numpy as np @@ -1658,8 +1559,6 @@ C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i, !et -===== Computation of gradients ===== - This in turn means that the gradient can be computed as a sum over $i$-gradients !bt @@ -1676,7 +1575,7 @@ minibatches. We denote these minibatches by $B_k$ where $k=1,\cdots,n/M$. -===== SGD example ===== + As an example, suppose we have $10$ data points $(\mathbf{x}_1,\cdots, \mathbf{x}_{10})$ and we choose to have $M=5$ minibathces, then each minibatch contains two data points. In particular we have @@ -1699,8 +1598,6 @@ c_i(\mathbf{x}_i, \mathbf{\beta}). !et -===== The gradient step ===== - Thus a gradient descent step now looks like !bt \[ @@ -1715,9 +1612,6 @@ minibathces (n/M) is commonly referred to as an epoch. Thus it is typical to choose a number of epochs and for each epoch iterate over the number of minibatches, as exemplified in the code below. - -===== Simple example code ===== - !bc pycod import numpy as np @@ -1744,7 +1638,6 @@ cheaper since we sum over the datapoints in the $k-th$ minibatch and not all $n$ datapoints. -===== When do we stop? ===== A natural question is when do we stop the search for a new minimum? One possibility is to compute the full gradient after a given number @@ -1758,7 +1651,6 @@ compare the values of the cost function and keep the $\beta$ that gave the lowest value. -===== Slightly different approach ===== Another approach is to let the step length $\gamma_j$ depend on the number of epochs in such a way that it becomes very small after a @@ -1800,11 +1692,7 @@ print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j)) !ec - - - - -===== Program for stochastic gradient ===== +=== Program for stochastic gradient === !bc pycod # Importing various packages @@ -1874,12 +1762,584 @@ plt.show() !ec -_Challenge_: try to write a similar code for a Logistic Regression case. - - - - - + +===== Momentum based GD ===== + +The stochastic gradient descent (SGD) is almost always used with a +*momentum* or inertia term that serves as a memory of the direction we +are moving in parameter space. This is typically implemented as +follows + +!bt +\begin{align} +\mathbf{v}_{t}&=\gamma \mathbf{v}_{t-1}+\eta_{t}\nabla_\theta E(\boldsymbol{\theta}_t) \nonumber \\ +\boldsymbol{\theta}_{t+1}&= \boldsymbol{\theta}_t -\mathbf{v}_{t}, +\end{align} +!et + +where we have introduced a momentum parameter $\gamma$, with +$0\le\gamma\le 1$, and for brevity we dropped the explicit notation to +indicate the gradient is to be taken over a different mini-batch at +each step. We call this algorithm gradient descent with momentum +(GDM). From these equations, it is clear that $\mathbf{v}_t$ is a +running average of recently encountered gradients and +$(1-\gamma)^{-1}$ sets the characteristic time scale for the memory +used in the averaging procedure. Consistent with this, when +$\gamma=0$, this just reduces down to ordinary SGD as discussed +earlier. An equivalent way of writing the updates is + +!bt +\[ +\Delta \boldsymbol{\theta}_{t+1} = \gamma \Delta \boldsymbol{\theta}_t -\ \eta_{t}\nabla_\theta E(\boldsymbol{\theta}_t), +\] +!et +where we have defined $\Delta \boldsymbol{\theta}_{t}= \boldsymbol{\theta}_t-\boldsymbol{\theta}_{t-1}$. + + + +Let us try to get more intuition from these equations. It is helpful +to consider a simple physical analogy with a particle of mass $m$ +moving in a viscous medium with drag coefficient $\mu$ and potential +$E(\mathbf{w})$. If we denote the particle's position by $\mathbf{w}$, +then its motion is described by + +!bt +\[ +m {d^2 \mathbf{w} \over dt^2} + \mu {d \mathbf{w} \over dt }= -\nabla_w E(\mathbf{w}). +\] +!et + +We can discretize this equation in the usual way to get + +!bt +\[ +m { \mathbf{w}_{t+\Delta t}-2 \mathbf{w}_{t} +\mathbf{w}_{t-\Delta t} \over (\Delta t)^2}+\mu {\mathbf{w}_{t+\Delta t}- \mathbf{w}_{t} \over \Delta t} = -\nabla_w E(\mathbf{w}). +\] +!et + +Rearranging this equation, we can rewrite this as + +!bt +\[ +\Delta \mathbf{w}_{t +\Delta t}= - { (\Delta t)^2 \over m +\mu \Delta t} \nabla_w E(\mathbf{w})+ {m \over m +\mu \Delta t} \Delta \mathbf{w}_t. +\] +!et + + + +Notice that this equation is identical to previous one if we identify +the position of the particle, $\mathbf{w}$, with the parameters +$\boldsymbol{\theta}$. This allows us to identify the momentum +parameter and learning rate with the mass of the particle and the +viscous drag as: + +!bt +\[ +\gamma= {m \over m +\mu \Delta t }, \qquad \eta = {(\Delta t)^2 \over m +\mu \Delta t}. +\] +!et + +Thus, as the name suggests, the momentum parameter is proportional to +the mass of the particle and effectively provides inertia. +Furthermore, in the large viscosity/small learning rate limit, our +memory time scales as $(1-\gamma)^{-1} \approx m/(\mu \Delta t)$. + +Why is momentum useful? SGD momentum helps the gradient descent +algorithm gain speed in directions with persistent but small gradients +even in the presence of stochasticity, while suppressing oscillations +in high-curvature directions. This becomes especially important in +situations where the landscape is shallow and flat in some directions +and narrow and steep in others. It has been argued that first-order +methods (with appropriate initial conditions) can perform comparable +to more expensive second order methods, especially in the context of +complex deep learning models. + +These beneficial properties of momentum can sometimes become even more +pronounced by using a slight modification of the classical momentum +algorithm called Nesterov Accelerated Gradient (NAG). + +In the NAG algorithm, rather than calculating the gradient at the +current parameters, $\nabla_\theta E(\boldsymbol{\theta}_t)$, one +calculates the gradient at the expected value of the parameters given +our current momentum, $\nabla_\theta E(\boldsymbol{\theta}_t +\gamma +\mathbf{v}_{t-1})$. This yields the NAG update rule + +!bt +\begin{align} +\mathbf{v}_{t}&=\gamma \mathbf{v}_{t-1}+\eta_{t}\nabla_\theta E(\boldsymbol{\theta}_t +\gamma \mathbf{v}_{t-1}) \nonumber \\ +\boldsymbol{\theta}_{t+1}&= \boldsymbol{\theta}_t -\mathbf{v}_{t}. +\end{align} +!et + +One of the major advantages of NAG is that it allows for the use of a larger learning rate than GDM for the same choice of $\gamma$. + + + +In stochastic gradient descent, with and without momentum, we still +have to specify a schedule for tuning the learning rates $\eta_t$ +as a function of time. As discussed in the context of Newton's +method, this presents a number of dilemmas. The learning rate is +limited by the steepest direction which can change depending on the +current position in the landscape. To circumvent this problem, ideally +our algorithm would keep track of curvature and take large steps in +shallow, flat directions and small steps in steep, narrow directions. +Second-order methods accomplish this by calculating or approximating +the Hessian and normalizing the learning rate by the +curvature. However, this is very computationally expensive for +extremely large models. Ideally, we would like to be able to +adaptively change the step size to match the landscape without paying +the steep computational price of calculating or approximating +Hessians. + +Recently, a number of methods have been introduced that accomplish +this by tracking not only the gradient, but also the second moment of +the gradient. These methods include AdaGrad, AdaDelta, RMS-Prop, and +ADAM. + + +=== RMS prop === + +In RMS prop, in addition to keeping a running average of the first +moment of the gradient, we also keep track of the second moment +denoted by $\mathbf{s}_t=\mathbb{E}[\mathbf{g}_t^2]$. The update rule +for RMS prop is given by + +!bt +\begin{align} +\mathbf{g}_t &= \nabla_\theta E(\boldsymbol{\theta}) \\ +\mathbf{s}_t &=\beta \mathbf{s}_{t-1} +(1-\beta)\mathbf{g}_t^2 \nonumber \\ +\boldsymbol{\theta}_{t+1}&=&\boldsymbol{\theta}_t - \eta_t { \mathbf{g}_t \over \sqrt{\mathbf{s}_t +\epsilon}}, \nonumber +\end{align} +!et + +where $\beta$ controls the averaging time of the second moment and is +typically taken to be about $\beta=0.9$, $\eta_t$ is a learning rate +typically chosen to be $10^{-3}$, and $\epsilon\sim 10^{-8} $ is a +small regularization constant to prevent divergences. Multiplication +and division by vectors is understood as an element-wise operation. It +is clear from this formula that the learning rate is reduced in +directions where the norm of the gradient is consistently large. This +greatly speeds up the convergence by allowing us to use a larger +learning rate for flat directions. + +=== ADAM optimizer === + +A related algorithm is the ADAM optimizer. In ADAM, we keep a running +average of both the first and second moment of the gradient and use +this information to adaptively change the learning rate for different +parameters. In addition to keeping a running average of the first and +second moments of the gradient +(i.e. $\mathbf{m}_t=\mathbb{E}[\mathbf{g}_t]$ and +$\mathbf{s}_t=\mathbb{E}[\mathbf{g}^2_t]$, respectively), ADAM +performs an additional bias correction to account for the fact that we +are estimating the first two moments of the gradient using a running +average (denoted by the hats in the update rule below). The update +rule for ADAM is given by (where multiplication and division are once +again understood to be element-wise operations below) + +!bt +\begin{align} +\mathbf{g}_t &= \nabla_\theta E(\boldsymbol{\theta}) \\ +\mathbf{m}_t &= \beta_1 \mathbf{m}_{t-1} + (1-\beta_1) \mathbf{g}_t \nonumber \\ +\mathbf{s}_t &=\beta_2 \mathbf{s}_{t-1} +(1-\beta_2)\mathbf{g}_t^2 \nonumber \\ +\bm{\mathbf{m}}_t&={\mathbf{m}_t \over 1-\beta_1^t} \nonumber \\ +\bm{\mathbf{s}}_t &={\mathbf{s}_t \over1-\beta_2^t} \nonumber \\ +\boldsymbol{\theta}_{t+1}&=\boldsymbol{\theta}_t - \eta_t { \bm{\mathbf{m}}_t \over \sqrt{\bm{\mathbf{s}}_t} +\epsilon}, \nonumber \\ +\end{align} +!et + +where $\beta_1$ and $\beta_2$ set the memory lifetime of the first and +second moment and are typically taken to be $0.9$ and $0.99$ +respectively, and $\eta$ and $\epsilon$ are identical to RMSprop. + +Like in RMSprop, the effective step size of a parameter depends on the +magnitude of its gradient squared. To understand this better, let us +rewrite this expression in terms of the variance +$\boldsymbol{\sigma}_t^2 = \bm{\mathbf{s}}_t - +(\bm{\mathbf{m}}_t)^2$. Consider a single parameter $\theta_t$. The +update rule for this parameter is given by + +!bt +\[ +\Delta \theta_{t+1}= -\eta_t { \bm{m}_t \over \sqrt{\sigma_t^2 + m_t^2 }+\epsilon}. +\] +!et + + + + + +===== Practical tips ===== + +* _Randomize the data when making mini-batches_. It is always important to randomly shuffle the data when forming mini-batches. Otherwise, the gradient descent method can fit spurious correlations resulting from the order in which data is presented. + +* _Transform your inputs_. Learning becomes difficult when our landscape has a mixture of steep and flat directions. One simple trick for minimizing these situations is to standardize the data by subtracting the mean and normalizing the variance of input variables. Whenever possible, also decorrelate the inputs. To understand why this is helpful, consider the case of linear regression. It is easy to show that for the squared error cost function, the Hessian of the cost function is just the correlation matrix between the inputs. Thus, by standardizing the inputs, we are ensuring that the landscape looks homogeneous in all directions in parameter space. Since most deep networks can be viewed as linear transformations followed by a non-linearity at each layer, we expect this intuition to hold beyond the linear case. + +* _Monitor the out-of-sample performance._ Always monitor the performance of your model on a validation set (a small portion of the training data that is held out of the training process to serve as a proxy for the test set. If the validation error starts increasing, then the model is beginning to overfit. Terminate the learning process. This *early stopping* significantly improves performance in many settings. + +* _Adaptive optimization methods don't always have good generalization._ Recent studies have shown that adaptive methods such as ADAM, RMSPorp, and AdaGrad tend to have poor generalization compared to SGD or SGD with momentum, particularly in the high-dimensional limit (i.e. the number of parameters exceeds the number of data points). Although it is not clear at this stage why these methods perform so well in training deep neural networks, simpler procedures like properly-tuned SGD may work as well or better in these applications. + + +===== Automatic differentiation ===== + +"Automatic differentiation (AD)":"https://en.wikipedia.org/wiki/Automatic_differentiation", +also called algorithmic +differentiation or computational differentiation,is a set of +techniques to numerically evaluate the derivative of a function +specified by a computer program. AD exploits the fact that every +computer program, no matter how complicated, executes a sequence of +elementary arithmetic operations (addition, subtraction, +multiplication, division, etc.) and elementary functions (exp, log, +sin, cos, etc.). By applying the chain rule repeatedly to these +operations, derivatives of arbitrary order can be computed +automatically, accurately to working precision, and using at most a +small constant factor more arithmetic operations than the original +program. + +Automatic differentiation is neither: + +* Symbolic differentiation, nor +* Numerical differentiation (the method of finite differences). + +Symbolic differentiation can lead to inefficient code and faces the +difficulty of converting a computer program into a single expression, +while numerical differentiation can introduce round-off errors in the +discretization process and cancellation + + + +Python has tools for so-called _automatic differentiation_. +Consider the following example +!bt +\[ +f(x) = \sin\left(2\pi x + x^2\right) +\] +!et +which has the following derivative +!bt +\[ +f'(x) = \cos\left(2\pi x + x^2\right)\left(2\pi + 2x\right) +\] +!et +Using _autograd_ we have + +!bc pycod +import autograd.numpy as np + +# To do elementwise differentiation: +from autograd import elementwise_grad as egrad + +# To plot: +import matplotlib.pyplot as plt + + +def f(x): + return np.sin(2*np.pi*x + x**2) + +def f_grad_analytic(x): + return np.cos(2*np.pi*x + x**2)*(2*np.pi + 2*x) + +# Do the comparison: +x = np.linspace(0,1,1000) + +f_grad = egrad(f) + +computed = f_grad(x) +analytic = f_grad_analytic(x) + +plt.title('Derivative computed from Autograd compared with the analytical derivative') +plt.plot(x,computed,label='autograd') +plt.plot(x,analytic,label='analytic') + +plt.xlabel('x') +plt.ylabel('y') +plt.legend() + +plt.show() + +print("The max absolute difference is: %g"%(np.max(np.abs(computed - analytic)))) +!ec + + +Here we +experiment with what kind of functions Autograd is capable +of finding the gradient of. The following Python functions are just +meant to illustrate what Autograd can do, but please feel free to +experiment with other, possibly more complicated, functions as well. + +!bc pycod +import autograd.numpy as np +from autograd import grad + +def f1(x): + return x**3 + 1 + +f1_grad = grad(f1) + +# Remember to send in float as argument to the computed gradient from Autograd! +a = 1.0 + +# See the evaluated gradient at a using autograd: +print("The gradient of f1 evaluated at a = %g using autograd is: %g"%(a,f1_grad(a))) + +# Compare with the analytical derivative, that is f1'(x) = 3*x**2 +grad_analytical = 3*a**2 +print("The gradient of f1 evaluated at a = %g by finding the analytic expression is: %g"%(a,grad_analytical)) +!ec + + + +To differentiate with respect to two (or more) arguments of a Python +function, Autograd need to know at which variable the function if +being differentiated with respect to. + +!bc pycod +import autograd.numpy as np +from autograd import grad +def f2(x1,x2): + return 3*x1**3 + x2*(x1 - 5) + 1 + +# By sending the argument 0, Autograd will compute the derivative w.r.t the first variable, in this case x1 +f2_grad_x1 = grad(f2,0) + +# ... and differentiate w.r.t x2 by sending 1 as an additional arugment to grad +f2_grad_x2 = grad(f2,1) + +x1 = 1.0 +x2 = 3.0 + +print("Evaluating at x1 = %g, x2 = %g"%(x1,x2)) +print("-"*30) + +# Compare with the analytical derivatives: + +# Derivative of f2 w.r.t x1 is: 9*x1**2 + x2: +f2_grad_x1_analytical = 9*x1**2 + x2 + +# Derivative of f2 w.r.t x2 is: x1 - 5: +f2_grad_x2_analytical = x1 - 5 + +# See the evaluated derivations: +print("The derivative of f2 w.r.t x1: %g"%( f2_grad_x1(x1,x2) )) +print("The analytical derivative of f2 w.r.t x1: %g"%( f2_grad_x1(x1,x2) )) + +print() + +print("The derivative of f2 w.r.t x2: %g"%( f2_grad_x2(x1,x2) )) +print("The analytical derivative of f2 w.r.t x2: %g"%( f2_grad_x2(x1,x2) )) +!ec + +Note that the grad function will not produce the true gradient of the function. The true gradient of a function with two or more variables will produce a vector, where each element is the function differentiated w.r.t a variable. + + + +!bc pycod +import autograd.numpy as np +from autograd import grad +def f3(x): # Assumes x is an array of length 5 or higher + return 2*x[0] + 3*x[1] + 5*x[2] + 7*x[3] + 11*x[4]**2 + +f3_grad = grad(f3) + +x = np.linspace(0,4,5) + +# Print the computed gradient: +print("The computed gradient of f3 is: ", f3_grad(x)) + +# The analytical gradient is: (2, 3, 5, 7, 22*x[4]) +f3_grad_analytical = np.array([2, 3, 5, 7, 22*x[4]]) + +# Print the analytical gradient: +print("The analytical gradient of f3 is: ", f3_grad_analytical) +!ec + +Note that in this case, when sending an array as input argument, the +output from Autograd is another array. This is the true gradient of +the function, as opposed to the function in the previous example. By +using arrays to represent the variables, the output from Autograd +might be easier to work with, as the output is closer to what one +could expect form a gradient-evaluting function. + + + +!bc pycod +import autograd.numpy as np +from autograd import grad +def f4(x): + return np.sqrt(1+x**2) + np.exp(x) + np.sin(2*np.pi*x) + +f4_grad = grad(f4) + +x = 2.7 + +# Print the computed derivative: +print("The computed derivative of f4 at x = %g is: %g"%(x,f4_grad(x))) + +# The analytical derivative is: x/sqrt(1 + x**2) + exp(x) + cos(2*pi*x)*2*pi +f4_grad_analytical = x/np.sqrt(1 + x**2) + np.exp(x) + np.cos(2*np.pi*x)*2*np.pi + +# Print the analytical gradient: +print("The analytical gradient of f4 at x = %g is: %g"%(x,f4_grad_analytical)) +!ec + + + +!bc pycod +import autograd.numpy as np +from autograd import grad +def f5(x): + if x >= 0: + return x**2 + else: + return -3*x + 1 + +f5_grad = grad(f5) + +x = 2.7 + +# Print the computed derivative: +print("The computed derivative of f5 at x = %g is: %g"%(x,f5_grad(x))) +!ec + + + +!bc pycod +import autograd.numpy as np +from autograd import grad +def f6_for(x): + val = 0 + for i in range(10): + val = val + x**i + return val + +def f6_while(x): + val = 0 + i = 0 + while i < 10: + val = val + x**i + i = i + 1 + return val + +f6_for_grad = grad(f6_for) +f6_while_grad = grad(f6_while) + +x = 0.5 + +# Print the computed derivaties of f6_for and f6_while +print("The computed derivative of f6_for at x = %g is: %g"%(x,f6_for_grad(x))) +print("The computed derivative of f6_while at x = %g is: %g"%(x,f6_while_grad(x))) +!ec +!bc pycod +import autograd.numpy as np +from autograd import grad +# Both of the functions are implementation of the sum: sum(x**i) for i = 0, ..., 9 +# The analytical derivative is: sum(i*x**(i-1)) +f6_grad_analytical = 0 +for i in range(10): + f6_grad_analytical += i*x**(i-1) + +print("The analytical derivative of f6 at x = %g is: %g"%(x,f6_grad_analytical)) +!ec + + +!bc pycod +import autograd.numpy as np +from autograd import grad + +def f7(n): # Assume that n is an integer + if n == 1 or n == 0: + return 1 + else: + return n*f7(n-1) + +f7_grad = grad(f7) + +n = 2.0 + +print("The computed derivative of f7 at n = %d is: %g"%(n,f7_grad(n))) + +# The function f7 is an implementation of the factorial of n. +# By using the product rule, one can find that the derivative is: + +f7_grad_analytical = 0 +for i in range(int(n)-1): + tmp = 1 + for k in range(int(n)-1): + if k != i: + tmp *= (n - k) + f7_grad_analytical += tmp + +print("The analytical derivative of f7 at n = %d is: %g"%(n,f7_grad_analytical)) + +!ec +Note that if n is equal to zero or one, Autograd will give an error message. This message appears when the output is independent on input. + + +Autograd supports many features. However, there are some functions that is not supported (yet) by Autograd. + +Assigning a value to the variable being differentiated with respect to +!bc pycod +import autograd.numpy as np +from autograd import grad +def f8(x): # Assume x is an array + x[2] = 3 + return x*2 + +f8_grad = grad(f8) + +x = 8.4 + +print("The derivative of f8 is:",f8_grad(x)) +!ec +Here, Autograd tells us that an 'ArrayBox' does not support item assignment. The item assignment is done when the program tries to assign x[2] to the value 3. However, Autograd has implemented the computation of the derivative such that this assignment is not possible. + + +!bc pycod +import autograd.numpy as np +from autograd import grad +def f9(a): # Assume a is an array with 2 elements + b = np.array([1.0,2.0]) + return a.dot(b) + +f9_grad = grad(f9) + +x = np.array([1.0,0.0]) + +print("The derivative of f9 is:",f9_grad(x)) +!ec + +Here we are told that the 'dot' function does not belong to Autograd's +version of a Numpy array. To overcome this, an alternative syntax +which also computed the dot product can be used: + +!bc pycod +import autograd.numpy as np +from autograd import grad +def f9_alternative(x): # Assume a is an array with 2 elements + b = np.array([1.0,2.0]) + return np.dot(x,b) # The same as x_1*b_1 + x_2*b_2 + +f9_alternative_grad = grad(f9_alternative) + +x = np.array([3.0,0.0]) + +print("The gradient of f9 is:",f9_alternative_grad(x)) + +# The analytical gradient of the dot product of vectors x and b with two elements (x_1,x_2) and (b_1, b_2) respectively +# w.r.t x is (b_1, b_2). +!ec + + +The documentation recommends to avoid inplace operations such as +!bc pycod +a += b +a -= b +a*= b +a /=b +!ec + + +More examples will be added, in particular how to compare autograd with own codes for the gradients. diff --git a/doc/LectureNotes/chapter1.ipynb b/doc/LectureNotes/chapter1.ipynb index c930e0a6b..9dab26cdc 100644 --- a/doc/LectureNotes/chapter1.ipynb +++ b/doc/LectureNotes/chapter1.ipynb @@ -13,38 +13,37 @@ "\n", "\n", "\n", - "Our emphasis throughout this series of lectures \n", - "is on understanding the mathematical aspects of\n", - "different algorithms used in the fields of data analysis and machine learning. \n", + "Our emphasis throughout this series of lectures is on understanding\n", + "the mathematical aspects of different algorithms used in the fields of\n", + "data analysis and machine learning.\n", "\n", - "However, where possible we will emphasize the\n", - "importance of using available software. We start thus with a hands-on\n", - "and top-down approach to machine learning. The aim is thus to start with\n", - "relevant data or data we have produced \n", - "and use these to introduce statistical data analysis\n", - "concepts and machine learning algorithms before we delve into the\n", - "algorithms themselves. The examples we will use in the beginning, start with simple\n", - "polynomials with random noise added. We will use the Python\n", - "software package [Scikit-Learn](http://scikit-learn.org/stable/) and\n", - "introduce various machine learning algorithms to make fits of\n", - "the data and predictions. We move thereafter to more interesting\n", - "cases such as data from say experiments (below we will look at experimental nuclear binding energies as an example).\n", - "These are examples where we can easily set up the data and\n", - "then use machine learning algorithms included in for example\n", - "**Scikit-Learn**. \n", + "However, where possible we will emphasize the importance of using\n", + "available software. We start thus with a hands-on and top-down\n", + "approach to machine learning. The aim is thus to start with relevant\n", + "data or data we have produced and use these to introduce statistical\n", + "data analysis concepts and machine learning algorithms before we delve\n", + "into the algorithms themselves. The examples we will use in the\n", + "beginning, start with simple polynomials with random noise added. We\n", + "will use the Python software package\n", + "[Scikit-Learn](http://scikit-learn.org/stable/) and introduce various\n", + "machine learning algorithms to make fits of the data and\n", + "predictions. We move thereafter to more interesting cases such as data\n", + "from say experiments (below we will look at experimental nuclear\n", + "binding energies as an example). These are examples where we can\n", + "easily set up the data and then use machine learning algorithms\n", + "included in for example **Scikit-Learn**.\n", "\n", "These examples will serve us the purpose of getting\n", "started. Furthermore, they allow us to catch more than two birds with\n", "a stone. They will allow us to bring in some programming specific\n", - "topics and tools as well as showing the power of various Python \n", - "libraries for machine learning and statistical data analysis. \n", + "topics and tools as well as showing the power of various Python\n", + "libraries for machine learning and statistical data analysis.\n", "\n", - "Here, we will mainly focus on two\n", - "specific Python packages for Machine Learning, Scikit-Learn and\n", - "Tensorflow (see below for links etc). Moreover, the examples we\n", - "introduce will serve as inputs to many of our discussions later, as\n", - "well as allowing you to set up models and produce your own data and\n", - "get started with programming.\n", + "Here, we will mainly focus on two specific Python packages for Machine\n", + "Learning, Scikit-Learn and Tensorflow (see below for links etc).\n", + "Moreover, the examples we introduce will serve as inputs to many of\n", + "our discussions later, as well as allowing you to set up models and\n", + "produce your own data and get started with programming.\n", "\n", "\n", "\n", @@ -286,7 +285,17 @@ "[gallery](https://matplotlib.org/gallery/index.html) of examples. In\n", "this example we plot our original values of $x$ and $y$ as well as the\n", "prediction **ypredict** ($\\tilde{y}$), which attempts at fitting our\n", - "data with a straight line.\n", + "data with a straight line. Note also that **Scikit-Learn** requires a\n", + "matrix as input for the input values $x$ and $y$. In the above code we\n", + "have solved this by declaring $x$ and $y$ as arrays of dimension\n", + "$n\\times 1$.\n", + "\n", + "In the code here we have also made a new array for $x\\in [0,1]$. Our\n", + "prediction is computed for these values, meaning that they were not\n", + "included in the data set used to *train* (or fit) the model.\n", + "This is a recurrring theme in machine learning and data analysis. We would like to train a model on a specific given data set.\n", + "Thereafter we wish to apply it to data which were not included in the training. Below we will encounter this again in the so-called *train-validate-test* spliting. We will typically split our data into different sets, oen for training, one for validation and finally, our data from the untouched test vault!\n", + "\n", "\n", "The Python code follows here." ] @@ -311,6 +320,7 @@ "y = 2*x+np.random.randn(100,1)\n", "linreg = LinearRegression()\n", "linreg.fit(x,y)\n", + "# This is our new x-array to which we test our model\n", "xnew = np.array([[0],[1]])\n", "ypredict = linreg.predict(xnew)\n", "\n", @@ -460,8 +470,7 @@ "Depending on the parameter in front of the normal distribution, we may\n", "have a small or larger relative error. Try to play around with\n", "different training data sets and study (graphically) the value of the\n", - "relative error. Note also that **Scikit-Learn** requires a matrix as input for the input values $x$ and $y$. In the above code we have\n", - "solved this by declaring $x$ and $y$ as arrays of dimension $n\\times 1$.\n", + "relative error.\n", "\n", "As mentioned above, **Scikit-Learn** has an impressive functionality.\n", "We can for example extract the values of $\\alpha$ and $\\beta$ and\n", @@ -1039,7 +1048,7 @@ "metadata": {}, "source": [ "The next step, and we will define this mathematically later, is to set up the so-called **design matrix**. We will throughout call this matrix $\\boldsymbol{X}$.\n", - "It has dimensionality $p\\times n$, where $n$ is the number of data points and $p$ are the so-called predictors. In our case here they are given by the number of polynomials in $A$ we wish to include in the fit." + "It has dimensionality $n\\times p$, where $n$ is the number of data points and $p$ are the so-called predictors. In our case here they are given by the number of polynomials in $A$ we wish to include in the fit." ] }, { @@ -1064,9 +1073,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note well that we have made life simple here. We perform a fit in terms of the number of nucleons only. A more sophisticated fit can be done by including an explicit dependence on the number of protons and neutrons in the asymmetry and Coulomb terms.\n", + "Note well that we have made life simple here. We perform a fit in\n", + "terms of the number of nucleons only. A more sophisticated fit can be\n", + "done by including an explicit dependence on the number of protons and\n", + "neutrons in the asymmetry and Coulomb terms. We leave this as an exercise to you the reader.\n", "\n", - "With **scikitlearn** we are now ready to use linear regression and fit our data." + "With **Scikit-Learn** we are now ready to use linear regression and fit our data." ] }, { @@ -1124,7 +1136,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "As a teaser, let us now see how we can do this with decision trees using **scikit-learn**. Later we will switch to so-called **random forests**!" + "As a teaser, let us now see how we can do this with decision trees using **Scikit-Learn**. Later we will switch to so-called **random forests**!" ] }, { @@ -1172,6 +1184,17 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "With a deeper and deeper tree level, we can almost reproduce every\n", + "single data point by increasing the max depth of the tree.\n", + "We can actually decide to make a decision tree which fits every single point.\n", + "As we will\n", + "see later, this has the benefit that we can really train a model which\n", + "traverses every single data point. However, the price we pay is that\n", + "we will easily overfit. That is, if we apply our model to unseen data,\n", + "we will most likely fail miserably in our attempt at making\n", + "predictions. As an exercise, try to make the tree level larger by adjusting the maximum depth variable. When printing out the predicition, you will note that the binding energy of every nucleus is accurately reproduced.\n", + "\n", + "\n", "The **seaborn** package allows us to visualize data in an efficient way. Note that we use **scikit-learn**'s multi-layer perceptron (or feed forward neural network) \n", "functionality." ] @@ -1977,6 +2000,47 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "We can then compute the second derivative of the cost function, which in our case is the second derivative\n", + "of the means squared error. This leads to" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\frac{\\partial^2 C(\\boldsymbol{\\beta})}{\\partial \\boldsymbol{\\beta}^T\\partial \\boldsymbol{\\beta}} =\\frac{2}{n}\\boldsymbol{X}^T\\boldsymbol{X}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This quantity defines was what is called the Hessian matrix (the second derivative of a function we want to optimize).\n", + "\n", + "The Hessian matrix plays an important role and is defined for the mean squared error as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{H}=\\boldsymbol{X}^T\\boldsymbol{X}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Hessian matrix for ordinary least squares is also proportional to\n", + "the covariance matrix. As we will see in the chapter on Ridge and Lasso regression, This means that we can use the Singular Value Decomposition of a matrix to find\n", + "the eigenvalues of the covariance matrix and the Hessian matrix in\n", + "terms of the singular values.\n", + "\n", + "\n", "The residuals $\\boldsymbol{\\epsilon}$ are in turn given by" ] }, @@ -2529,9 +2593,7 @@ "**pandas** again, rather extensively in order to organize our data.\n", "\n", "The difference now is that we use **Scikit-Learn's** regression tools\n", - "instead of our own matrix inversion implementation. Furthermore, we\n", - "sneak in **Ridge** regression (to be discussed below) which includes a\n", - "hyperparameter $\\lambda$, also to be explained below." + "instead of our own matrix inversion implementation." ] }, { @@ -2603,18 +2665,6 @@ "print('Mean absolute error: %.2f' % mean_absolute_error(Energies, ytilde))\n", "print(clf.coef_, clf.intercept_)\n", "\n", - "# The Ridge regression with a hyperparameter lambda = 0.1\n", - "_lambda = 0.1\n", - "clf_ridge = skl.Ridge(alpha=_lambda).fit(X, Energies)\n", - "yridge = clf_ridge.predict(X)\n", - "EoS['Eridge'] = yridge\n", - "# The mean squared error \n", - "print(\"Mean squared error: %.2f\" % mean_squared_error(Energies, yridge))\n", - "# Explained variance score: 1 is perfect prediction \n", - "print('Variance score: %.2f' % r2_score(Energies, yridge))\n", - "# Mean absolute error \n", - "print('Mean absolute error: %.2f' % mean_absolute_error(Energies, yridge))\n", - "print(clf_ridge.coef_, clf_ridge.intercept_)\n", "\n", "fig, ax = plt.subplots()\n", "ax.set_xlabel(r'$\\rho[\\mathrm{fm}^{-3}]$')\n", @@ -2623,8 +2673,6 @@ " label='Theoretical data')\n", "ax.plot(EoS['Density'], EoS['Eols'], alpha=0.7, lw=2, c='m',\n", " label='OLS')\n", - "ax.plot(EoS['Density'], EoS['Eridge'], alpha=0.7, lw=2, c='g',\n", - " label='Ridge $\\lambda = 0.1$')\n", "ax.legend()\n", "save_fig(\"EoSfitting\")\n", "plt.show()" @@ -2637,13 +2685,11 @@ "The above simple polynomial in density $\\rho$ gives an excellent fit\n", "to the data. \n", "\n", - "We note also that there is a small deviation between the\n", - "standard OLS and the Ridge regression at higher densities. We discuss this in more detail\n", - "below.\n", "\n", "\n", "## Splitting our Data in Training and Test data\n", "\n", + "\n", "It is normal in essentially all Machine Learning studies to split the\n", "data in a training set and a test set (sometimes also an additional\n", "validation set). **Scikit-Learn** has an own function for this. There\n", @@ -2653,6 +2699,107 @@ "postpone a discussion of this splitting to the end of these notes and\n", "our discussion of the so-called **bias-variance** tradeoff. Here we\n", "limit ourselves to repeat the above equation of state fitting example\n", + "but now splitting the data into a training set and a test set.\n", + "\n", + "Let us study some examples. The first code here takes a simple\n", + "one-dimensional second-order polynomial and we fit it to a\n", + "second-order polynomial. Depending on the strength of the added noise,\n", + "the various measures like the $R2$ score or the mean-squared error,\n", + "the fit becomes better or worse." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "from sklearn.model_selection import train_test_split\n", + "\n", + "\n", + "def R2(y_data, y_model):\n", + " return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n", + "def MSE(y_data,y_model):\n", + " n = np.size(y_model)\n", + " return np.sum((y_data-y_model)**2)/n\n", + "\n", + "x = np.random.rand(100)\n", + "y = 2.0+5*x*x+0.1*np.random.randn(100)\n", + "\n", + "\n", + "# The design matrix now as function of a given polynomial\n", + "X = np.zeros((len(x),3))\n", + "X[:,0] = 1.0\n", + "X[:,1] = x\n", + "X[:,2] = x**2\n", + "# We split the data in test and training data\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", + "# matrix inversion to find beta\n", + "beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n", + "print(beta)\n", + "# and then make the prediction\n", + "ytilde = X_train @ beta\n", + "print(\"Training R2\")\n", + "print(R2(y_train,ytilde))\n", + "print(\"Training MSE\")\n", + "print(MSE(y_train,ytilde))\n", + "ypredict = X_test @ beta\n", + "print(\"Test R2\")\n", + "print(R2(y_test,ypredict))\n", + "print(\"Test MSE\")\n", + "print(MSE(y_test,ypredict))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Alternatively, you could write your own test-train splitting function as shown here." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "# equivalently in numpy\n", + "def train_test_split_numpy(inputs, labels, train_size, test_size):\n", + " n_inputs = len(inputs)\n", + " inputs_shuffled = inputs.copy()\n", + " labels_shuffled = labels.copy()\n", + "\n", + " np.random.shuffle(inputs_shuffled)\n", + " np.random.shuffle(labels_shuffled)\n", + "\n", + " train_end = int(n_inputs*train_size)\n", + " X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:]\n", + " Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:]\n", + "\n", + " return X_train, X_test, Y_train, Y_test" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "But since **scikit-learn** has its own function for doing this and since\n", + "it interfaces easily with **tensorflow** and other libraries, we\n", + "normally recommend using the latter functionality.\n", + "\n", + "\n", + "As another example, we apply the training and testing split to \n", + "to the above equation of state fitting example\n", "but now splitting the data into a training set and a test set." ] }, @@ -3050,119 +3197,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Splitting our Data in Training and Test data\n", - "\n", - "\n", - "It is normal in essentially all Machine Learning studies to split the\n", - "data in a training set and a test set (sometimes also an additional\n", - "validation set). **Scikit-Learn** has an own function for this. There\n", - "is no explicit recipe for how much data should be included as training\n", - "data and say test data. An accepted rule of thumb is to use\n", - "approximately $2/3$ to $4/5$ of the data as training data. We will\n", - "postpone a discussion of this splitting to the end of these notes and\n", - "our discussion of the so-called **bias-variance** tradeoff. Here we\n", - "limit ourselves to repeat the above equation of state fitting example\n", - "but now splitting the data into a training set and a test set.\n", - "\n", - "Let us study some examples. The first code here takes a simple\n", - "one-dimensional second-order polynomial and we fit it to a\n", - "second-order polynomial. Depending on the strength of the added noise,\n", - "the various measures like the $R2$ score or the mean-squared error,\n", - "the fit becomes better or worse." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], - "source": [ - "import os\n", - "import numpy as np\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "from sklearn.model_selection import train_test_split\n", - "\n", - "\n", - "def R2(y_data, y_model):\n", - " return 1 - np.sum((y_data - y_model) ** 2) / np.sum((y_data - np.mean(y_data)) ** 2)\n", - "def MSE(y_data,y_model):\n", - " n = np.size(y_model)\n", - " return np.sum((y_data-y_model)**2)/n\n", - "\n", - "x = np.random.rand(100)\n", - "y = 2.0+5*x*x+0.1*np.random.randn(100)\n", - "\n", - "\n", - "# The design matrix now as function of a given polynomial\n", - "X = np.zeros((len(x),3))\n", - "X[:,0] = 1.0\n", - "X[:,1] = x\n", - "X[:,2] = x**2\n", - "# We split the data in test and training data\n", - "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)\n", - "# matrix inversion to find beta\n", - "beta = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train\n", - "print(beta)\n", - "# and then make the prediction\n", - "ytilde = X_train @ beta\n", - "print(\"Training R2\")\n", - "print(R2(y_train,ytilde))\n", - "print(\"Training MSE\")\n", - "print(MSE(y_train,ytilde))\n", - "ypredict = X_test @ beta\n", - "print(\"Test R2\")\n", - "print(R2(y_test,ypredict))\n", - "print(\"Test MSE\")\n", - "print(MSE(y_test,ypredict))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Alternatively, you could write your own test-train splitting function as shown here." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], - "source": [ - "# equivalently in numpy\n", - "def train_test_split_numpy(inputs, labels, train_size, test_size):\n", - " n_inputs = len(inputs)\n", - " inputs_shuffled = inputs.copy()\n", - " labels_shuffled = labels.copy()\n", - "\n", - " np.random.shuffle(inputs_shuffled)\n", - " np.random.shuffle(labels_shuffled)\n", - "\n", - " train_end = int(n_inputs*train_size)\n", - " X_train, X_test = inputs_shuffled[:train_end], inputs_shuffled[train_end:]\n", - " Y_train, Y_test = labels_shuffled[:train_end], labels_shuffled[train_end:]\n", - "\n", - " return X_train, X_test, Y_train, Y_test" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "But since **scikit-learn** has its own function for doing this and since\n", - "it interfaces easily with **tensorflow** and other libraries, we\n", - "normally recommend using the latter functionality.\n", - "\n", - "\n", - "\n", - "\n", "## Reducing the number of degrees of freedom, overarching view\n", "\n", "Many Machine Learning problems involve thousands or even millions of\n", diff --git a/doc/LectureNotes/chapter2.ipynb b/doc/LectureNotes/chapter2.ipynb index 651eadf4a..3049b1030 100644 --- a/doc/LectureNotes/chapter2.ipynb +++ b/doc/LectureNotes/chapter2.ipynb @@ -13,7 +13,7 @@ "What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD). \n", "\n", "\n", - "We have shown that in ordinary least squares the optimal parameters $\\beta$ are given by" + "We have shown that in ordinary least squares (OLS) the optimal parameters $\\beta$ are given by" ] }, { @@ -21,7 +21,7 @@ "metadata": {}, "source": [ "$$\n", - "\\hat{\\boldsymbol{\\beta}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}.\n", + "\\hat{\\boldsymbol{\\beta}}_{\\mathrm{OLS}} = \\left(\\boldsymbol{X}^T\\boldsymbol{X}\\right)^{-1}\\boldsymbol{X}^T\\boldsymbol{y}.\n", "$$" ] }, @@ -79,7 +79,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The matrix $\\boldsymbol{A}$ has the important property that $\\boldsymbol{A}^2=\\boldsymbol{A}$. This is the definition of a projection matrix.\n", + "The matrix $\\boldsymbol{A}$ has the important property that $\\boldsymbol{A}^2=\\boldsymbol{A}$. This is the definition of a [projection matrix](https://en.wikipedia.org/wiki/Projection_matrix).\n", "We can then interpret our optimal model $\\tilde{\\boldsymbol{y}}$ as being represented by an orthogonal projection of $\\boldsymbol{y}$ onto a space defined by the column vectors of $\\boldsymbol{X}$. In our case here the matrix $\\boldsymbol{A}$ is a square matrix. If it is a general rectangular matrix we have an oblique projection matrix.\n", "\n", "\n", @@ -545,7 +545,7 @@ "\n", "## Code for SVD and Inversion of Matrices\n", "\n", - "How do we use the SVD to invert a matrix $\\boldsymbol{X}^\\boldsymbol{X}$ which is singular or near singular?\n", + "How do we use the SVD to invert a matrix $\\boldsymbol{X}^T\\boldsymbol{X}$ which is singular or near singular?\n", "The simple answer is to use the linear algebra function for the pseudoinverse, that is" ] }, @@ -3258,6 +3258,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "**More text will be added to this example.**\n", + "\n", "## Linking the regression analysis with a statistical interpretation\n", "\n", "We will now couple the discussions of ordinary least squares, Ridge\n", @@ -3624,7 +3626,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "It is a conditional probability (see below) and reads as the likelihood of a domain of events $\\boldsymbol{D}$ given a set of parameters $\\boldsymbol{\\beta}$.\n", + "It is a conditional probability (see below) and reads as the\n", + "likelihood of a domain of events $\\boldsymbol{D}$ given a set of parameters\n", + "$\\boldsymbol{\\beta}$.\n", "\n", "\n", "In statistics, maximum likelihood estimation (MLE) is a method of\n", @@ -3728,7 +3732,7 @@ "We define also the likelihood for $X$ and $Y$ as $p(X)$ and $p(Y)$ respectively.\n", "The likelihood of a specific event $x_i$ (or $y_i$) is then written as $p(X=x_i)$ or just $p(x_i)=p_i$. \n", "\n", - "**Union of events is given by.**" + "The union of events is given by" ] }, { @@ -3744,7 +3748,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "**The product rule (aka joint probability) is given by.**" + "The product rule (aka joint probability) is given by" ] }, { @@ -3762,8 +3766,6 @@ "source": [ "where we read $p(X\\vert Y)$ as the likelihood of obtaining $X$ given $Y$.\n", "\n", - "\n", - "\n", "If we have independent events then $p(X,Y)=p(X)p(Y)$.\n", "\n", "\n", @@ -3840,21 +3842,27 @@ "the parameter space represented by $X$. This function is not\n", "necesseraly normalized and is normally called the likelihood function.\n", "\n", - "The function $p(X)$ on the right hand side is called the prior while the function on the left hand side is the called the posterior probability. The denominator on the right hand side serves as a normalization factor for the posterior distribution.\n", + "The function $p(X)$ on the right hand side is called the prior while\n", + "the function on the left hand side is the called the posterior\n", + "probability. The denominator on the right hand side serves as a\n", + "normalization factor for the posterior distribution.\n", "\n", "Let us try to illustrate Bayes' theorem through an example.\n", "\n", "\n", - "Let us suppose that you are undergoing a series of mammography scans in\n", - "order to rule out possible breast cancer cases. We define the\n", + "Let us suppose that you are undergoing a series of mammography scans\n", + "in order to rule out possible breast cancer cases. We define the\n", "sensitivity for a positive event by the variable $X$. It takes binary\n", "values with $X=1$ representing a positive event and $X=0$ being a\n", "negative event. We reserve $Y$ as a classification parameter for\n", - "either a negative or a positive breast cancer confirmation. (Short note on wordings: positive here means having breast cancer, although none of us would consider this being a positive thing).\n", + "either a negative or a positive breast cancer confirmation. (Short\n", + "note on wordings: positive here means having breast cancer, although\n", + "none of us would consider this being a positive thing).\n", "\n", "We let $Y=1$ represent the the case of having breast cancer and $Y=0$ as not.\n", "\n", - "Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$, that is we have" + "Let us assume that if you have breast cancer, the test will be positive with a probability of $0.8$ (the numbers here are all made up),\n", + "that is we have" ] }, { @@ -3870,12 +3878,32 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This obviously sounds scary since many would conclude that if the test is positive, there is a likelihood of $80\\%$ for having cancer.\n", - "It is however not correct, as the following Bayesian analysis shows.\n", + "This obviously sounds scary since many would conclude that if the test\n", + "is positive, there is a likelihood of $80\\%$ for having cancer. It is\n", + "however not correct, as the following Bayesian analysis shows. The correct question to pose is *what is the probability of having breast cancer in case of a positive test?*\n", + "We are thus interested in" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "p(Y=1\\vert X=1),\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "instead of $p(X=1\\vert Y=1)$.\n", "\n", "\n", - "If we look at various national surveys on breast cancer, the general likelihood of developing breast cancer is a very small number.\n", - "Let us assume that the prior probability in the population as a whole is" + "\n", + "If we look at various national surveys on breast cancer, the general\n", + "likelihood of developing breast cancer is a very small number. Let us\n", + "assume that the prior probability in the population as a whole is" ] }, { @@ -3891,7 +3919,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We need also to account for the fact that the test may produce a false positive result (false alarm). Let us here assume that we have" + "We need also to account for the fact that the test may produce a false\n", + "positive result (false alarm). Let us here assume that we have" ] }, { @@ -3907,15 +3936,38 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Using Bayes' theorem we can then find the posterior probability that the person has breast cancer in case of a positive test, that is we can compute" + "Using Bayes' theorem we can then find the posterior probability that\n", + "the person has breast cancer in case of a positive test, that is we\n", + "can compute" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ + "\n", + "
\n", + "\n", "$$\n", - "p(Y=1\\vert X=1)=\\frac{p(X=1\\vert Y=1)p(Y=1)}{p(X=1\\vert Y=1)p(Y=1)+p(X=1\\vert Y=0)p(Y=0)}=\\frac{0.8\\times 0.004}{0.8\\times 0.004+0.1\\times 0.996}=0.031.\n", + "\\begin{equation}\n", + "p(Y=1\\vert X=1)=\\frac{p(X=1\\vert Y=1)p(Y=1)}{p(X=1\\vert Y=1)p(Y=1)+p(X=1\\vert Y=0)p(Y=0)}= \n", + "\\label{_auto2} \\tag{2}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + " \\frac{0.8\\times 0.004}{0.8\\times 0.004+0.1\\times 0.996}=0.031.\n", + "\\label{_auto3} \\tag{3}\n", + "\\end{equation}\n", "$$" ] }, @@ -3941,9 +3993,12 @@ "parameters $\\boldsymbol{\\beta}$ using OLS, Ridge and Lasso regression. For\n", "OLS, you will notice as function of the noise and polynomial degree,\n", "that the parameters $\\beta$ will fluctuate from order to order in the\n", - "polynomial fit and that for larger and larger polynomial degrees of freedom, the parameters will tend to increase in value for OLS.\n", + "polynomial fit and that for larger and larger polynomial degrees of\n", + "freedom, the parameters will tend to increase in value for OLS.\n", "\n", - "For Ridge and Lasso regression, the higher order parameters will typically be reduced, providing thereby less fluctuations from one order to another one." + "For Ridge and Lasso regression, the higher order parameters will\n", + "typically be reduced, providing thereby less fluctuations from one\n", + "order to another one." ] }, { @@ -4008,8 +4063,8 @@ " # Compute the MSE and print it\n", " MSERidgePredict[i] = MSE(y_test,ypredictRidge)\n", " MSELassoPredict[i] = MSE(y_test,ypredictLasso)\n", - " print(lmb,RegRidge.coef_)\n", - " print(lmb,RegLasso.coef_)\n", + "# print(lmb,RegRidge.coef_)\n", + "# print(lmb,RegLasso.coef_)\n", "# Now plot the results\n", "plt.figure()\n", "plt.plot(np.log10(lambdas), MSERidgePredict, 'b', label = 'MSE Ridge Test')\n", @@ -4024,10 +4079,13 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "How can we understand this? \n", + "How can we understand this? **More text to be added**.\n", "\n", "\n", "\n", + "## Linking Bayes' Theorem with Ridge and Lasso Regression\n", + "\n", + "\n", "Using Bayes' theorem we can gain a better intuition about Ridge and Lasso regression. \n", "\n", "For ordinary least squares we postulated that the maximum likelihood for the doamin of events $\\boldsymbol{D}$ (one-dimensional case)" diff --git a/doc/LectureNotes/chapter3.ipynb b/doc/LectureNotes/chapter3.ipynb index d7f8222d5..b2933e098 100644 --- a/doc/LectureNotes/chapter3.ipynb +++ b/doc/LectureNotes/chapter3.ipynb @@ -575,196 +575,6 @@ "plt.show()" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Various steps in cross-validation\n", - "\n", - "When the repetitive splitting of the data set is done randomly,\n", - "samples may accidently end up in a fast majority of the splits in\n", - "either training or test set. Such samples may have an unbalanced\n", - "influence on either model building or prediction evaluation. To avoid\n", - "this $k$-fold cross-validation structures the data splitting. The\n", - "samples are divided into $k$ more or less equally sized exhaustive and\n", - "mutually exclusive subsets. In turn (at each split) one of these\n", - "subsets plays the role of the test set while the union of the\n", - "remaining subsets constitutes the training set. Such a splitting\n", - "warrants a balanced representation of each sample in both training and\n", - "test set over the splits. Still the division into the $k$ subsets\n", - "involves a degree of randomness. This may be fully excluded when\n", - "choosing $k=n$. This particular case is referred to as leave-one-out\n", - "cross-validation (LOOCV). \n", - "\n", - "\n", - "* Define a range of interest for the penalty parameter.\n", - "\n", - "* Divide the data set into training and test set comprising samples $\\{1, \\ldots, n\\} \\setminus i$ and $\\{ i \\}$, respectively.\n", - "\n", - "* Fit the linear regression model by means of ridge estimation for each $\\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\\boldsymbol{\\sigma}_{-i}^2(\\lambda)$, as" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "$$\n", - "\\begin{align*}\n", - "\\boldsymbol{\\beta}_{-i}(\\lambda) & = ( \\boldsymbol{X}_{-i, \\ast}^{T}\n", - "\\boldsymbol{X}_{-i, \\ast} + \\lambda \\boldsymbol{I}_{pp})^{-1}\n", - "\\boldsymbol{X}_{-i, \\ast}^{T} \\boldsymbol{y}_{-i}\n", - "\\end{align*}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "* Evaluate the prediction performance of these models on the test set by $\\log\\{L[y_i, \\boldsymbol{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}$. Or, by the prediction error $|y_i - \\boldsymbol{X}_{i, \\ast} \\boldsymbol{\\beta}_{-i}(\\lambda)|$, the relative error, the error squared or the R2 score function.\n", - "\n", - "* Repeat the first three steps such that each sample plays the role of the test set once.\n", - "\n", - "* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "$$\n", - "\\begin{align*}\n", - "\\frac{1}{n} \\sum_{i = 1}^n \\log\\{L[y_i, \\mathbf{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}.\n", - "\\end{align*}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For the various values of $k$\n", - "\n", - "1. shuffle the dataset randomly.\n", - "\n", - "2. Split the dataset into $k$ groups.\n", - "\n", - "3. For each unique group:\n", - "\n", - "a. Decide which group to use as set for test data\n", - "\n", - "b. Take the remaining groups as a training data set\n", - "\n", - "c. Fit a model on the training set and evaluate it on the test set\n", - "\n", - "d. Retain the evaluation score and discard the model\n", - "\n", - "\n", - "5. Summarize the model using the sample of model evaluation scores\n", - "\n", - "The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "from sklearn.model_selection import KFold\n", - "from sklearn.linear_model import Ridge\n", - "from sklearn.model_selection import cross_val_score\n", - "from sklearn.preprocessing import PolynomialFeatures\n", - "\n", - "# A seed just to ensure that the random numbers are the same for every run.\n", - "# Useful for eventual debugging.\n", - "np.random.seed(3155)\n", - "\n", - "# Generate the data.\n", - "nsamples = 100\n", - "x = np.random.randn(nsamples)\n", - "y = 3*x**2 + np.random.randn(nsamples)\n", - "\n", - "## Cross-validation on Ridge regression using KFold only\n", - "\n", - "# Decide degree on polynomial to fit\n", - "poly = PolynomialFeatures(degree = 6)\n", - "\n", - "# Decide which values of lambda to use\n", - "nlambdas = 500\n", - "lambdas = np.logspace(-3, 5, nlambdas)\n", - "\n", - "# Initialize a KFold instance\n", - "k = 5\n", - "kfold = KFold(n_splits = k)\n", - "\n", - "# Perform the cross-validation to estimate MSE\n", - "scores_KFold = np.zeros((nlambdas, k))\n", - "\n", - "i = 0\n", - "for lmb in lambdas:\n", - " ridge = Ridge(alpha = lmb)\n", - " j = 0\n", - " for train_inds, test_inds in kfold.split(x):\n", - " xtrain = x[train_inds]\n", - " ytrain = y[train_inds]\n", - "\n", - " xtest = x[test_inds]\n", - " ytest = y[test_inds]\n", - "\n", - " Xtrain = poly.fit_transform(xtrain[:, np.newaxis])\n", - " ridge.fit(Xtrain, ytrain[:, np.newaxis])\n", - "\n", - " Xtest = poly.fit_transform(xtest[:, np.newaxis])\n", - " ypred = ridge.predict(Xtest)\n", - "\n", - " scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)\n", - "\n", - " j += 1\n", - " i += 1\n", - "\n", - "\n", - "estimated_mse_KFold = np.mean(scores_KFold, axis = 1)\n", - "\n", - "## Cross-validation using cross_val_score from sklearn along with KFold\n", - "\n", - "# kfold is an instance initialized above as:\n", - "# kfold = KFold(n_splits = k)\n", - "\n", - "estimated_mse_sklearn = np.zeros(nlambdas)\n", - "i = 0\n", - "for lmb in lambdas:\n", - " ridge = Ridge(alpha = lmb)\n", - "\n", - " X = poly.fit_transform(x[:, np.newaxis])\n", - " estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)\n", - "\n", - " # cross_val_score return an array containing the estimated negative mse for every fold.\n", - " # we have to the the mean of every array in order to get an estimate of the mse of the model\n", - " estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n", - "\n", - " i += 1\n", - "\n", - "## Plot and compare the slightly different ways to perform cross-validation\n", - "\n", - "plt.figure()\n", - "\n", - "plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n", - "plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')\n", - "\n", - "plt.xlabel('log10(lambda)')\n", - "plt.ylabel('mse')\n", - "\n", - "plt.legend()\n", - "\n", - "plt.show()" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -1221,6 +1031,203 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Cross-validation\n", + "\n", + "When the repetitive splitting of the data set is done randomly,\n", + "samples may accidently end up in a fast majority of the splits in\n", + "either training or test set. Such samples may have an unbalanced\n", + "influence on either model building or prediction evaluation. To avoid\n", + "this $k$-fold cross-validation structures the data splitting. The\n", + "samples are divided into $k$ more or less equally sized exhaustive and\n", + "mutually exclusive subsets. In turn (at each split) one of these\n", + "subsets plays the role of the test set while the union of the\n", + "remaining subsets constitutes the training set. Such a splitting\n", + "warrants a balanced representation of each sample in both training and\n", + "test set over the splits. Still the division into the $k$ subsets\n", + "involves a degree of randomness. This may be fully excluded when\n", + "choosing $k=n$. This particular case is referred to as leave-one-out\n", + "cross-validation (LOOCV). \n", + "\n", + "\n", + "* Define a range of interest for the penalty parameter.\n", + "\n", + "* Divide the data set into training and test set comprising samples $\\{1, \\ldots, n\\} \\setminus i$ and $\\{ i \\}$, respectively.\n", + "\n", + "* Fit the linear regression model by means of ridge estimation for each $\\lambda$ in the grid using the training set, and the corresponding estimate of the error variance $\\boldsymbol{\\sigma}_{-i}^2(\\lambda)$, as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "\\boldsymbol{\\beta}_{-i}(\\lambda) & = ( \\boldsymbol{X}_{-i, \\ast}^{T}\n", + "\\boldsymbol{X}_{-i, \\ast} + \\lambda \\boldsymbol{I}_{pp})^{-1}\n", + "\\boldsymbol{X}_{-i, \\ast}^{T} \\boldsymbol{y}_{-i}\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* Evaluate the prediction performance of these models on the test set by $\\log\\{L[y_i, \\boldsymbol{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}$. Or, by the prediction error $|y_i - \\boldsymbol{X}_{i, \\ast} \\boldsymbol{\\beta}_{-i}(\\lambda)|$, the relative error, the error squared or the R2 score function.\n", + "\n", + "* Repeat the first three steps such that each sample plays the role of the test set once.\n", + "\n", + "* Average the prediction performances of the test sets at each grid point of the penalty bias/parameter. It is an estimate of the prediction performance of the model corresponding to this value of the penalty parameter on novel data. It is defined as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\begin{align*}\n", + "\\frac{1}{n} \\sum_{i = 1}^n \\log\\{L[y_i, \\mathbf{X}_{i, \\ast}; \\boldsymbol{\\beta}_{-i}(\\lambda), \\boldsymbol{\\sigma}_{-i}^2(\\lambda)]\\}.\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the various values of $k$\n", + "\n", + "1. shuffle the dataset randomly.\n", + "\n", + "2. Split the dataset into $k$ groups.\n", + "\n", + "3. For each unique group:\n", + "\n", + "a. Decide which group to use as set for test data\n", + "\n", + "b. Take the remaining groups as a training data set\n", + "\n", + "c. Fit a model on the training set and evaluate it on the test set\n", + "\n", + "d. Retain the evaluation score and discard the model\n", + "\n", + "\n", + "5. Summarize the model using the sample of model evaluation scores\n", + "\n", + "The code here uses Ridge regression with cross-validation (CV) resampling and $k$-fold CV in order to fit a specific polynomial." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from sklearn.model_selection import KFold\n", + "from sklearn.linear_model import Ridge\n", + "from sklearn.model_selection import cross_val_score\n", + "from sklearn.preprocessing import PolynomialFeatures\n", + "\n", + "# A seed just to ensure that the random numbers are the same for every run.\n", + "# Useful for eventual debugging.\n", + "np.random.seed(3155)\n", + "\n", + "# Generate the data.\n", + "nsamples = 100\n", + "x = np.random.randn(nsamples)\n", + "y = 3*x**2 + np.random.randn(nsamples)\n", + "\n", + "## Cross-validation on Ridge regression using KFold only\n", + "\n", + "# Decide degree on polynomial to fit\n", + "poly = PolynomialFeatures(degree = 6)\n", + "\n", + "# Decide which values of lambda to use\n", + "nlambdas = 500\n", + "lambdas = np.logspace(-3, 5, nlambdas)\n", + "\n", + "# Initialize a KFold instance\n", + "k = 5\n", + "kfold = KFold(n_splits = k)\n", + "\n", + "# Perform the cross-validation to estimate MSE\n", + "scores_KFold = np.zeros((nlambdas, k))\n", + "\n", + "i = 0\n", + "for lmb in lambdas:\n", + " ridge = Ridge(alpha = lmb)\n", + " j = 0\n", + " for train_inds, test_inds in kfold.split(x):\n", + " xtrain = x[train_inds]\n", + " ytrain = y[train_inds]\n", + "\n", + " xtest = x[test_inds]\n", + " ytest = y[test_inds]\n", + "\n", + " Xtrain = poly.fit_transform(xtrain[:, np.newaxis])\n", + " ridge.fit(Xtrain, ytrain[:, np.newaxis])\n", + "\n", + " Xtest = poly.fit_transform(xtest[:, np.newaxis])\n", + " ypred = ridge.predict(Xtest)\n", + "\n", + " scores_KFold[i,j] = np.sum((ypred - ytest[:, np.newaxis])**2)/np.size(ypred)\n", + "\n", + " j += 1\n", + " i += 1\n", + "\n", + "\n", + "estimated_mse_KFold = np.mean(scores_KFold, axis = 1)\n", + "\n", + "## Cross-validation using cross_val_score from sklearn along with KFold\n", + "\n", + "# kfold is an instance initialized above as:\n", + "# kfold = KFold(n_splits = k)\n", + "\n", + "estimated_mse_sklearn = np.zeros(nlambdas)\n", + "i = 0\n", + "for lmb in lambdas:\n", + " ridge = Ridge(alpha = lmb)\n", + "\n", + " X = poly.fit_transform(x[:, np.newaxis])\n", + " estimated_mse_folds = cross_val_score(ridge, X, y[:, np.newaxis], scoring='neg_mean_squared_error', cv=kfold)\n", + "\n", + " # cross_val_score return an array containing the estimated negative mse for every fold.\n", + " # we have to the the mean of every array in order to get an estimate of the mse of the model\n", + " estimated_mse_sklearn[i] = np.mean(-estimated_mse_folds)\n", + "\n", + " i += 1\n", + "\n", + "## Plot and compare the slightly different ways to perform cross-validation\n", + "\n", + "plt.figure()\n", + "\n", + "plt.plot(np.log10(lambdas), estimated_mse_sklearn, label = 'cross_val_score')\n", + "plt.plot(np.log10(lambdas), estimated_mse_KFold, 'r--', label = 'KFold')\n", + "\n", + "plt.xlabel('log10(lambda)')\n", + "plt.ylabel('mse')\n", + "\n", + "plt.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "More examples of the application of cross-validation follow here." + ] + }, { "cell_type": "code", "execution_count": null, diff --git a/doc/LectureNotes/chapter4.ipynb b/doc/LectureNotes/chapter4.ipynb index 7cd60a191..81b9d76c8 100644 --- a/doc/LectureNotes/chapter4.ipynb +++ b/doc/LectureNotes/chapter4.ipynb @@ -7,8 +7,6 @@ "# Logistic Regression\n", "\n", "\n", - "[Video of Lecture](https://www.uio.no/studier/emner/matnat/fys/FYS-STK3155/h20/forelesningsvideoer/LectureSeptember18.mp4?vrtx=view-as-webpage)\n", - "\n", "\n", "## Logistic Regression\n", "\n", @@ -911,7 +909,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Optimization, the central part of any Machine Learning algortithm\n", + "# Optimization, the central part of any Machine Learning algortithm\n", "\n", "Almost every problem in machine learning and data science starts with\n", "a dataset $X$, a model $g(\\beta)$, which is a function of the\n", @@ -924,7 +922,6 @@ "\n", "\n", "\n", - "## Revisiting our Logistic Regression case\n", "\n", "In our discussion on Logistic Regression we studied the \n", "case of\n", @@ -953,7 +950,6 @@ "where $\\boldsymbol{\\beta}$ are the weights we wish to extract from data, in our case $\\beta_0$ and $\\beta_1$. \n", "\n", "\n", - "## The equations to solve\n", "\n", "Our compact equations used a definition of a vector $\\boldsymbol{y}$ with $n$\n", "elements $y_i$, an $n\\times p$ matrix $\\boldsymbol{X}$ which contains the\n", @@ -995,7 +991,6 @@ "This defines what is called the Hessian matrix.\n", "\n", "\n", - "## Solving using Newton-Raphson's method\n", "\n", "If we can set up these equations, Newton-Raphson's iterative method is normally the method of choice. It requires however that we can compute in an efficient way the matrices that define the first and second derivatives. \n", "\n", @@ -1036,9 +1031,6 @@ "If we can compute these matrices, in particular the Hessian, the above is often the easiest method to implement. \n", "\n", "\n", - "\n", - "## Brief reminder on Newton-Raphson's method\n", - "\n", "Let us quickly remind ourselves how we derive the above method.\n", "\n", "Perhaps the most celebrated of all one-dimensional root-finding\n", @@ -1050,7 +1042,6 @@ "normally discourage the use of this method.\n", "\n", "\n", - "## The equations\n", "\n", "The Newton-Raphson formula consists geometrically of extending the\n", "tangent line at a current point until it crosses zero, then setting\n", @@ -1125,8 +1116,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Simple geometric interpretation\n", - "\n", "The above is Newton-Raphson's method. It has a simple geometric\n", "interpretation, namely $x_{n+1}$ is the point where the tangent from\n", "$(x_n,f(x_n))$ crosses the $x$-axis. Close to the solution,\n", @@ -1141,7 +1130,6 @@ "\n", "\n", "\n", - "## Extending to more than one variable\n", "\n", "Newton's method can be generalized to systems of several non-linear equations\n", "and variables. Consider the case with two equations" @@ -1278,8 +1266,6 @@ "we are always moving towards smaller function values, i.e a minimum.\n", "\n", "\n", - "## More on Steepest descent\n", - "\n", "The previous observation is the basis of the method of steepest\n", "descent, which is also referred to as just gradient descent (GD). One\n", "starts with an initial guess $\\mathbf{x}_0$ for a minimum of $F$ and\n", @@ -1303,8 +1289,6 @@ "the learning rate within the context of Machine Learning.\n", "\n", "\n", - "## The ideal\n", - "\n", "Ideally the sequence $\\{\\mathbf{x}_k \\}_{k=0}$ converges to a global\n", "minimum of the function $F$. In general we do not know if we are in a\n", "global or local minimum. In the special case when $F$ is a convex\n", @@ -1325,7 +1309,6 @@ "\n", "\n", "\n", - "## The sensitiveness of the gradient descent\n", "\n", "The gradient descent method \n", "is sensitive to the choice of learning rate $\\gamma_k$. This is due\n", @@ -1356,13 +1339,10 @@ "regular polygons (triangles, rectangles, pentagons, etc...).\n", "\n", "\n", - "## Convex function\n", "\n", "**Convex function**: Let $X \\subset \\mathbb{R}^n$ be a convex set. Assume that the function $f: X \\rightarrow \\mathbb{R}$ is continuous, then $f$ is said to be convex if $$f(tx_1 + (1-t)x_2) \\leq tf(x_1) + (1-t)f(x_2) $$ for all $x_1, x_2 \\in X$ and for all $t \\in [0,1]$. If $\\leq$ is replaced with a strict inequaltiy in the definition, we demand $x_1 \\neq x_2$ and $t\\in(0,1)$ then $f$ is said to be strictly convex. For a single variable function, convexity means that if you draw a straight line connecting $f(x_1)$ and $f(x_2)$, the value of the function on the interval $[x_1,x_2]$ is always below the line as illustrated below.\n", "\n", "\n", - "## Conditions on convex functions\n", - "\n", "In the following we state first and second-order conditions which\n", "ensures convexity of a function $f$. We write $D_f$ to denote the\n", "domain of $f$, i.e the subset of $R^n$ where $f$ is defined. For more\n", @@ -1395,7 +1375,6 @@ "This condition is particularly useful since it gives us an procedure for determining if the function under consideration is convex, apart from using the definition.\n", "\n", "\n", - "## More on convex functions\n", "\n", "The next result is of great importance to us and the reason why we are\n", "going on about convex functions. In machine learning we frequently\n", @@ -1418,7 +1397,7 @@ "This result means that if we know that the cost/loss function is convex and we are able to find a minimum, we are guaranteed that it is a global minimum.\n", "\n", "\n", - "## Some simple problems\n", + "### Some simple problems\n", "\n", "1. Show that $f(x)=x^2$ is convex for $x \\in \\mathbb{R}$ using the definition of convexity. Hint: If you re-write the definition, $f$ is convex if the following holds for all $x,y \\in D_f$ and any $\\lambda \\in [0,1]$ $\\lambda f(x)+(1-\\lambda)f(y)-f(\\lambda x + (1-\\lambda) y ) \\geq 0$.\n", "\n", @@ -1443,13 +1422,6 @@ "Using the definition of convexity, try to show that a function satisfying the properties above is convex (the third condition is not needed to show this).\n", "\n", "\n", - "\n", - "## Friday September 25\n", - "\n", - "[Video of Lecture](https://www.uio.no/studier/emner/matnat/fys/FYS-STK4155/h20/forelesningsvideoer/LectureSeptember25.mp4?vrtx=view-as-webpage) and [link to handwritten notes](https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/NotesSeptember25.pdf).\n", - "\n", - "\n", - "\n", "## Standard steepest descent\n", "\n", "\n", @@ -1498,7 +1470,6 @@ "When we have found the exact solution, $\\boldsymbol{r}=0$.\n", "\n", "\n", - "## Gradient method\n", "\n", "The residual is zero when we reach the minimum of the quadratic equation" ] @@ -1520,9 +1491,6 @@ "symmetric. This defines also the Hessian and we want it to be positive definite. \n", "\n", "\n", - "\n", - "## Steepest descent method\n", - "\n", "We denote the initial guess for $\\boldsymbol{x}$ as $\\boldsymbol{x}_0$. \n", "We can assume without loss of generality that" ] @@ -1559,8 +1527,6 @@ "instead.\n", "\n", "\n", - "\n", - "## Steepest descent method\n", "One can show that the solution $\\boldsymbol{x}$ is also the unique minimizer of the quadratic form" ] }, @@ -1599,9 +1565,6 @@ "$\\boldsymbol{x}_0=0$ it is equal $-\\boldsymbol{b}$.\n", "\n", "\n", - "\n", - "\n", - "## Final expressions\n", "We can compute the residual iteratively as" ] }, @@ -1678,13 +1641,6 @@ "$$" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Steepest descent example" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1848,10 +1804,6 @@ "Two vectors are conjugate if they are orthogonal with respect to \n", "this inner product. Being conjugate is a symmetric relation: if $\\boldsymbol{s}$ is conjugate to $\\boldsymbol{t}$, then $\\boldsymbol{t}$ is conjugate to $\\boldsymbol{s}$.\n", "\n", - "\n", - "\n", - "\n", - "## Conjugate gradient method\n", "An example is given by the eigenvectors of the matrix" ] }, @@ -1868,13 +1820,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "which is zero unless $i=j$.\n", + "which is zero unless $i=j$. \n", "\n", - "\n", - "\n", - "\n", - "\n", - "## Conjugate gradient method\n", "Assume now that we have a symmetric positive-definite matrix $\\boldsymbol{A}$ of size\n", "$n\\times n$. At each iteration $i+1$ we obtain the conjugate direction of a vector" ] @@ -1910,7 +1857,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Conjugate gradient method\n", "The coefficients are given by" ] }, @@ -1959,8 +1905,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Conjugate gradient method and iterations\n", - "\n", "If we choose the conjugate vectors $\\boldsymbol{p}_k$ carefully, \n", "then we may not need all of them to obtain a good approximation to the solution \n", "$\\boldsymbol{x}$. \n", @@ -2003,11 +1947,6 @@ "source": [ "instead.\n", "\n", - "\n", - "\n", - "\n", - "\n", - "## Conjugate gradient method\n", "One can show that the solution $\\boldsymbol{x}$ is also the unique minimizer of the quadratic form" ] }, @@ -2047,11 +1986,6 @@ "The other vectors in the basis will be conjugate to the gradient, \n", "hence the name conjugate gradient method.\n", "\n", - "\n", - "\n", - "\n", - "\n", - "## Conjugate gradient method\n", "Let $\\boldsymbol{r}_k$ be the residual at the $k$-th step:" ] }, @@ -2090,7 +2024,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Conjugate gradient method\n", "We can also compute the residual iteratively as" ] }, @@ -2155,14 +2088,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Revisiting our first homework\n", + "## Revisiting our Linear Regression Solvers\n", "\n", "We will use linear regression as a case study for the gradient descent\n", "methods. Linear regression is a great test case for the gradient\n", "descent methods discussed in the lectures since it has several\n", "desirable properties such as:\n", "\n", - "1. An analytical solution (recall homework set 1).\n", + "1. An analytical solution.\n", "\n", "2. The gradient can be computed analytically.\n", "\n", @@ -2221,8 +2154,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Gradient descent example\n", - "\n", "Let $\\mathbf{y} = (y_1,\\cdots,y_n)^T$, $\\mathbf{\\boldsymbol{y}} = (\\boldsymbol{y}_1,\\cdots,\\boldsymbol{y}_n)^T$ and $\\beta = (\\beta_0, \\beta_1)^T$\n", "\n", "It is convenient to write $\\mathbf{\\boldsymbol{y}} = X\\beta$ where $X \\in \\mathbb{R}^{100 \\times 2} $ is the design matrix given by (we keep the intercept here)" @@ -2264,8 +2195,6 @@ "and we want to find $\\beta$ such that $C(\\beta)$ is minimized.\n", "\n", "\n", - "## The derivative of the cost/loss function\n", - "\n", "Computing $\\partial C(\\beta) / \\partial \\beta_0$ and $\\partial C(\\beta) / \\partial \\beta_1$ we can show that the gradient can be written as" ] }, @@ -2287,7 +2216,6 @@ "where $X$ is the design matrix defined above.\n", "\n", "\n", - "## The Hessian matrix\n", "The Hessian matrix of $C(\\beta)$ is given by" ] }, @@ -2311,10 +2239,6 @@ "\n", "\n", "\n", - "\n", - "\n", - "## Simple program\n", - "\n", "We can now write a program that minimizes $C(\\beta)$ using the gradient descent method with a constant learning rate $\\gamma$ according to" ] }, @@ -2339,9 +2263,8 @@ "$\\beta= (X^TX)^{-1} X^T \\mathbf{y}$.\n", "\n", "\n", - "## Gradient Descent Example\n", "\n", - "Here our simple example" + "Here is our simple example" ] }, { @@ -2405,7 +2328,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## And a corresponding example using **scikit-learn**" + "Alternatively, we can use **Scikit-Learn** as done here" ] }, { @@ -2439,8 +2362,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Gradient descent and Ridge\n", - "\n", "We have also discussed Ridge regression where the loss function contains a regularized term given by the $L_2$ norm of $\\beta$," ] }, @@ -2487,13 +2408,6 @@ "$$" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Program example for gradient descent with Ridge Regression" - ] - }, { "cell_type": "code", "execution_count": null, @@ -2590,8 +2504,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Computation of gradients\n", - "\n", "This in turn means that the gradient can be\n", "computed as a sum over $i$-gradients" ] @@ -2617,7 +2529,7 @@ "$k=1,\\cdots,n/M$.\n", "\n", "\n", - "## SGD example\n", + "\n", "As an example, suppose we have $10$ data points $(\\mathbf{x}_1,\\cdots, \\mathbf{x}_{10})$ \n", "and we choose to have $M=5$ minibathces,\n", "then each minibatch contains two data points. In particular we have\n", @@ -2648,8 +2560,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## The gradient step\n", - "\n", "Thus a gradient descent step now looks like" ] }, @@ -2671,10 +2581,7 @@ "probability from $[1,n/M]$. An iteration over the number of\n", "minibathces (n/M) is commonly referred to as an epoch. Thus it is\n", "typical to choose a number of epochs and for each epoch iterate over\n", - "the number of minibatches, as exemplified in the code below.\n", - "\n", - "\n", - "## Simple example code" + "the number of minibatches, as exemplified in the code below." ] }, { @@ -2715,7 +2622,6 @@ "all $n$ datapoints.\n", "\n", "\n", - "## When do we stop?\n", "\n", "A natural question is when do we stop the search for a new minimum?\n", "One possibility is to compute the full gradient after a given number\n", @@ -2729,7 +2635,6 @@ "gave the lowest value.\n", "\n", "\n", - "## Slightly different approach\n", "\n", "Another approach is to let the step length $\\gamma_j$ depend on the\n", "number of epochs in such a way that it becomes very small after a\n", @@ -2783,7 +2688,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Program for stochastic gradient" + "### Program for stochastic gradient" ] }, { @@ -2865,7 +2770,922 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "**Challenge**: try to write a similar code for a Logistic Regression case." + "## Momentum based GD\n", + "\n", + "The stochastic gradient descent (SGD) is almost always used with a\n", + "*momentum* or inertia term that serves as a memory of the direction we\n", + "are moving in parameter space. This is typically implemented as\n", + "follows" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathbf{v}_{t}=\\gamma \\mathbf{v}_{t-1}+\\eta_{t}\\nabla_\\theta E(\\boldsymbol{\\theta}_t) \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\boldsymbol{\\theta}_{t+1}= \\boldsymbol{\\theta}_t -\\mathbf{v}_{t},\n", + "\\label{_auto2} \\tag{3}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where we have introduced a momentum parameter $\\gamma$, with\n", + "$0\\le\\gamma\\le 1$, and for brevity we dropped the explicit notation to\n", + "indicate the gradient is to be taken over a different mini-batch at\n", + "each step. We call this algorithm gradient descent with momentum\n", + "(GDM). From these equations, it is clear that $\\mathbf{v}_t$ is a\n", + "running average of recently encountered gradients and\n", + "$(1-\\gamma)^{-1}$ sets the characteristic time scale for the memory\n", + "used in the averaging procedure. Consistent with this, when\n", + "$\\gamma=0$, this just reduces down to ordinary SGD as discussed\n", + "earlier. An equivalent way of writing the updates is" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\Delta \\boldsymbol{\\theta}_{t+1} = \\gamma \\Delta \\boldsymbol{\\theta}_t -\\ \\eta_{t}\\nabla_\\theta E(\\boldsymbol{\\theta}_t),\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where we have defined $\\Delta \\boldsymbol{\\theta}_{t}= \\boldsymbol{\\theta}_t-\\boldsymbol{\\theta}_{t-1}$.\n", + "\n", + "\n", + "\n", + "Let us try to get more intuition from these equations. It is helpful\n", + "to consider a simple physical analogy with a particle of mass $m$\n", + "moving in a viscous medium with drag coefficient $\\mu$ and potential\n", + "$E(\\mathbf{w})$. If we denote the particle's position by $\\mathbf{w}$,\n", + "then its motion is described by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "m {d^2 \\mathbf{w} \\over dt^2} + \\mu {d \\mathbf{w} \\over dt }= -\\nabla_w E(\\mathbf{w}).\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can discretize this equation in the usual way to get" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "m { \\mathbf{w}_{t+\\Delta t}-2 \\mathbf{w}_{t} +\\mathbf{w}_{t-\\Delta t} \\over (\\Delta t)^2}+\\mu {\\mathbf{w}_{t+\\Delta t}- \\mathbf{w}_{t} \\over \\Delta t} = -\\nabla_w E(\\mathbf{w}).\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Rearranging this equation, we can rewrite this as" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\Delta \\mathbf{w}_{t +\\Delta t}= - { (\\Delta t)^2 \\over m +\\mu \\Delta t} \\nabla_w E(\\mathbf{w})+ {m \\over m +\\mu \\Delta t} \\Delta \\mathbf{w}_t.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that this equation is identical to previous one if we identify\n", + "the position of the particle, $\\mathbf{w}$, with the parameters\n", + "$\\boldsymbol{\\theta}$. This allows us to identify the momentum\n", + "parameter and learning rate with the mass of the particle and the\n", + "viscous drag as:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\gamma= {m \\over m +\\mu \\Delta t }, \\qquad \\eta = {(\\Delta t)^2 \\over m +\\mu \\Delta t}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Thus, as the name suggests, the momentum parameter is proportional to\n", + "the mass of the particle and effectively provides inertia.\n", + "Furthermore, in the large viscosity/small learning rate limit, our\n", + "memory time scales as $(1-\\gamma)^{-1} \\approx m/(\\mu \\Delta t)$.\n", + "\n", + "Why is momentum useful? SGD momentum helps the gradient descent\n", + "algorithm gain speed in directions with persistent but small gradients\n", + "even in the presence of stochasticity, while suppressing oscillations\n", + "in high-curvature directions. This becomes especially important in\n", + "situations where the landscape is shallow and flat in some directions\n", + "and narrow and steep in others. It has been argued that first-order\n", + "methods (with appropriate initial conditions) can perform comparable\n", + "to more expensive second order methods, especially in the context of\n", + "complex deep learning models.\n", + "\n", + "These beneficial properties of momentum can sometimes become even more\n", + "pronounced by using a slight modification of the classical momentum\n", + "algorithm called Nesterov Accelerated Gradient (NAG).\n", + "\n", + "In the NAG algorithm, rather than calculating the gradient at the\n", + "current parameters, $\\nabla_\\theta E(\\boldsymbol{\\theta}_t)$, one\n", + "calculates the gradient at the expected value of the parameters given\n", + "our current momentum, $\\nabla_\\theta E(\\boldsymbol{\\theta}_t +\\gamma\n", + "\\mathbf{v}_{t-1})$. This yields the NAG update rule" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathbf{v}_{t}=\\gamma \\mathbf{v}_{t-1}+\\eta_{t}\\nabla_\\theta E(\\boldsymbol{\\theta}_t +\\gamma \\mathbf{v}_{t-1}) \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\boldsymbol{\\theta}_{t+1}= \\boldsymbol{\\theta}_t -\\mathbf{v}_{t}.\n", + "\\label{_auto3} \\tag{4}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "One of the major advantages of NAG is that it allows for the use of a larger learning rate than GDM for the same choice of $\\gamma$.\n", + "\n", + "\n", + "\n", + "In stochastic gradient descent, with and without momentum, we still\n", + "have to specify a schedule for tuning the learning rates $\\eta_t$\n", + "as a function of time. As discussed in the context of Newton's\n", + "method, this presents a number of dilemmas. The learning rate is\n", + "limited by the steepest direction which can change depending on the\n", + "current position in the landscape. To circumvent this problem, ideally\n", + "our algorithm would keep track of curvature and take large steps in\n", + "shallow, flat directions and small steps in steep, narrow directions.\n", + "Second-order methods accomplish this by calculating or approximating\n", + "the Hessian and normalizing the learning rate by the\n", + "curvature. However, this is very computationally expensive for\n", + "extremely large models. Ideally, we would like to be able to\n", + "adaptively change the step size to match the landscape without paying\n", + "the steep computational price of calculating or approximating\n", + "Hessians.\n", + "\n", + "Recently, a number of methods have been introduced that accomplish\n", + "this by tracking not only the gradient, but also the second moment of\n", + "the gradient. These methods include AdaGrad, AdaDelta, RMS-Prop, and\n", + "ADAM.\n", + "\n", + "\n", + "### RMS prop\n", + "\n", + "In RMS prop, in addition to keeping a running average of the first\n", + "moment of the gradient, we also keep track of the second moment\n", + "denoted by $\\mathbf{s}_t=\\mathbb{E}[\\mathbf{g}_t^2]$. The update rule\n", + "for RMS prop is given by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\mathbf{g}_t = \\nabla_\\theta E(\\boldsymbol{\\theta}) \n", + "\\label{_auto4} \\tag{5}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathbf{s}_t =\\beta \\mathbf{s}_{t-1} +(1-\\beta)\\mathbf{g}_t^2 \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{\\theta}_{t+1}=\\boldsymbol{\\theta}_t - \\eta_t { \\mathbf{g}_t \\over \\sqrt{\\mathbf{s}_t +\\epsilon}}, \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $\\beta$ controls the averaging time of the second moment and is\n", + "typically taken to be about $\\beta=0.9$, $\\eta_t$ is a learning rate\n", + "typically chosen to be $10^{-3}$, and $\\epsilon\\sim 10^{-8} $ is a\n", + "small regularization constant to prevent divergences. Multiplication\n", + "and division by vectors is understood as an element-wise operation. It\n", + "is clear from this formula that the learning rate is reduced in\n", + "directions where the norm of the gradient is consistently large. This\n", + "greatly speeds up the convergence by allowing us to use a larger\n", + "learning rate for flat directions.\n", + "\n", + "### ADAM optimizer\n", + "\n", + "A related algorithm is the ADAM optimizer. In ADAM, we keep a running\n", + "average of both the first and second moment of the gradient and use\n", + "this information to adaptively change the learning rate for different\n", + "parameters. In addition to keeping a running average of the first and\n", + "second moments of the gradient\n", + "(i.e. $\\mathbf{m}_t=\\mathbb{E}[\\mathbf{g}_t]$ and\n", + "$\\mathbf{s}_t=\\mathbb{E}[\\mathbf{g}^2_t]$, respectively), ADAM\n", + "performs an additional bias correction to account for the fact that we\n", + "are estimating the first two moments of the gradient using a running\n", + "average (denoted by the hats in the update rule below). The update\n", + "rule for ADAM is given by (where multiplication and division are once\n", + "again understood to be element-wise operations below)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation}\n", + "\\mathbf{g}_t = \\nabla_\\theta E(\\boldsymbol{\\theta}) \n", + "\\label{_auto5} \\tag{6}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathbf{m}_t = \\beta_1 \\mathbf{m}_{t-1} + (1-\\beta_1) \\mathbf{g}_t \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\mathbf{s}_t =\\beta_2 \\mathbf{s}_{t-1} +(1-\\beta_2)\\mathbf{g}_t^2 \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{\\mathbf{m}}_t={\\mathbf{m}_t \\over 1-\\beta_1^t} \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{\\mathbf{s}}_t ={\\mathbf{s}_t \\over1-\\beta_2^t} \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\boldsymbol{\\theta}_{t+1}=\\boldsymbol{\\theta}_t - \\eta_t { \\boldsymbol{\\mathbf{m}}_t \\over \\sqrt{\\boldsymbol{\\mathbf{s}}_t} +\\epsilon}, \\nonumber\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "
\n", + "\n", + "$$\n", + "\\begin{equation} \n", + "\\label{_auto6} \\tag{7}\n", + "\\end{equation}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where $\\beta_1$ and $\\beta_2$ set the memory lifetime of the first and\n", + "second moment and are typically taken to be $0.9$ and $0.99$\n", + "respectively, and $\\eta$ and $\\epsilon$ are identical to RMSprop.\n", + "\n", + "Like in RMSprop, the effective step size of a parameter depends on the\n", + "magnitude of its gradient squared. To understand this better, let us\n", + "rewrite this expression in terms of the variance\n", + "$\\boldsymbol{\\sigma}_t^2 = \\boldsymbol{\\mathbf{s}}_t -\n", + "(\\boldsymbol{\\mathbf{m}}_t)^2$. Consider a single parameter $\\theta_t$. The\n", + "update rule for this parameter is given by" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\Delta \\theta_{t+1}= -\\eta_t { \\boldsymbol{m}_t \\over \\sqrt{\\sigma_t^2 + m_t^2 }+\\epsilon}.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Practical tips\n", + "\n", + "* **Randomize the data when making mini-batches**. It is always important to randomly shuffle the data when forming mini-batches. Otherwise, the gradient descent method can fit spurious correlations resulting from the order in which data is presented.\n", + "\n", + "* **Transform your inputs**. Learning becomes difficult when our landscape has a mixture of steep and flat directions. One simple trick for minimizing these situations is to standardize the data by subtracting the mean and normalizing the variance of input variables. Whenever possible, also decorrelate the inputs. To understand why this is helpful, consider the case of linear regression. It is easy to show that for the squared error cost function, the Hessian of the cost function is just the correlation matrix between the inputs. Thus, by standardizing the inputs, we are ensuring that the landscape looks homogeneous in all directions in parameter space. Since most deep networks can be viewed as linear transformations followed by a non-linearity at each layer, we expect this intuition to hold beyond the linear case.\n", + "\n", + "* **Monitor the out-of-sample performance.** Always monitor the performance of your model on a validation set (a small portion of the training data that is held out of the training process to serve as a proxy for the test set. If the validation error starts increasing, then the model is beginning to overfit. Terminate the learning process. This *early stopping* significantly improves performance in many settings.\n", + "\n", + "* **Adaptive optimization methods don't always have good generalization.** Recent studies have shown that adaptive methods such as ADAM, RMSPorp, and AdaGrad tend to have poor generalization compared to SGD or SGD with momentum, particularly in the high-dimensional limit (i.e. the number of parameters exceeds the number of data points). Although it is not clear at this stage why these methods perform so well in training deep neural networks, simpler procedures like properly-tuned SGD may work as well or better in these applications.\n", + "\n", + "## Automatic differentiation\n", + "\n", + "[Automatic differentiation (AD)](https://en.wikipedia.org/wiki/Automatic_differentiation), \n", + "also called algorithmic\n", + "differentiation or computational differentiation,is a set of\n", + "techniques to numerically evaluate the derivative of a function\n", + "specified by a computer program. AD exploits the fact that every\n", + "computer program, no matter how complicated, executes a sequence of\n", + "elementary arithmetic operations (addition, subtraction,\n", + "multiplication, division, etc.) and elementary functions (exp, log,\n", + "sin, cos, etc.). By applying the chain rule repeatedly to these\n", + "operations, derivatives of arbitrary order can be computed\n", + "automatically, accurately to working precision, and using at most a\n", + "small constant factor more arithmetic operations than the original\n", + "program.\n", + "\n", + "Automatic differentiation is neither:\n", + "\n", + "* Symbolic differentiation, nor\n", + "\n", + "* Numerical differentiation (the method of finite differences).\n", + "\n", + "Symbolic differentiation can lead to inefficient code and faces the\n", + "difficulty of converting a computer program into a single expression,\n", + "while numerical differentiation can introduce round-off errors in the\n", + "discretization process and cancellation\n", + "\n", + "\n", + "\n", + "Python has tools for so-called **automatic differentiation**.\n", + "Consider the following example" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "f(x) = \\sin\\left(2\\pi x + x^2\\right)\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "which has the following derivative" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "f'(x) = \\cos\\left(2\\pi x + x^2\\right)\\left(2\\pi + 2x\\right)\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using **autograd** we have" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "\n", + "# To do elementwise differentiation:\n", + "from autograd import elementwise_grad as egrad \n", + "\n", + "# To plot:\n", + "import matplotlib.pyplot as plt \n", + "\n", + "\n", + "def f(x):\n", + " return np.sin(2*np.pi*x + x**2)\n", + "\n", + "def f_grad_analytic(x):\n", + " return np.cos(2*np.pi*x + x**2)*(2*np.pi + 2*x)\n", + "\n", + "# Do the comparison:\n", + "x = np.linspace(0,1,1000)\n", + "\n", + "f_grad = egrad(f)\n", + "\n", + "computed = f_grad(x)\n", + "analytic = f_grad_analytic(x)\n", + "\n", + "plt.title('Derivative computed from Autograd compared with the analytical derivative')\n", + "plt.plot(x,computed,label='autograd')\n", + "plt.plot(x,analytic,label='analytic')\n", + "\n", + "plt.xlabel('x')\n", + "plt.ylabel('y')\n", + "plt.legend()\n", + "\n", + "plt.show()\n", + "\n", + "print(\"The max absolute difference is: %g\"%(np.max(np.abs(computed - analytic))))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we\n", + "experiment with what kind of functions Autograd is capable\n", + "of finding the gradient of. The following Python functions are just\n", + "meant to illustrate what Autograd can do, but please feel free to\n", + "experiment with other, possibly more complicated, functions as well." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "\n", + "def f1(x):\n", + " return x**3 + 1\n", + "\n", + "f1_grad = grad(f1)\n", + "\n", + "# Remember to send in float as argument to the computed gradient from Autograd!\n", + "a = 1.0\n", + "\n", + "# See the evaluated gradient at a using autograd:\n", + "print(\"The gradient of f1 evaluated at a = %g using autograd is: %g\"%(a,f1_grad(a)))\n", + "\n", + "# Compare with the analytical derivative, that is f1'(x) = 3*x**2 \n", + "grad_analytical = 3*a**2\n", + "print(\"The gradient of f1 evaluated at a = %g by finding the analytic expression is: %g\"%(a,grad_analytical))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To differentiate with respect to two (or more) arguments of a Python\n", + "function, Autograd need to know at which variable the function if\n", + "being differentiated with respect to." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "def f2(x1,x2):\n", + " return 3*x1**3 + x2*(x1 - 5) + 1\n", + "\n", + "# By sending the argument 0, Autograd will compute the derivative w.r.t the first variable, in this case x1\n", + "f2_grad_x1 = grad(f2,0)\n", + "\n", + "# ... and differentiate w.r.t x2 by sending 1 as an additional arugment to grad\n", + "f2_grad_x2 = grad(f2,1)\n", + "\n", + "x1 = 1.0\n", + "x2 = 3.0 \n", + "\n", + "print(\"Evaluating at x1 = %g, x2 = %g\"%(x1,x2))\n", + "print(\"-\"*30)\n", + "\n", + "# Compare with the analytical derivatives:\n", + "\n", + "# Derivative of f2 w.r.t x1 is: 9*x1**2 + x2:\n", + "f2_grad_x1_analytical = 9*x1**2 + x2\n", + "\n", + "# Derivative of f2 w.r.t x2 is: x1 - 5:\n", + "f2_grad_x2_analytical = x1 - 5\n", + "\n", + "# See the evaluated derivations:\n", + "print(\"The derivative of f2 w.r.t x1: %g\"%( f2_grad_x1(x1,x2) ))\n", + "print(\"The analytical derivative of f2 w.r.t x1: %g\"%( f2_grad_x1(x1,x2) ))\n", + "\n", + "print()\n", + "\n", + "print(\"The derivative of f2 w.r.t x2: %g\"%( f2_grad_x2(x1,x2) ))\n", + "print(\"The analytical derivative of f2 w.r.t x2: %g\"%( f2_grad_x2(x1,x2) ))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the grad function will not produce the true gradient of the function. The true gradient of a function with two or more variables will produce a vector, where each element is the function differentiated w.r.t a variable." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "def f3(x): # Assumes x is an array of length 5 or higher\n", + " return 2*x[0] + 3*x[1] + 5*x[2] + 7*x[3] + 11*x[4]**2\n", + "\n", + "f3_grad = grad(f3)\n", + "\n", + "x = np.linspace(0,4,5)\n", + "\n", + "# Print the computed gradient:\n", + "print(\"The computed gradient of f3 is: \", f3_grad(x))\n", + "\n", + "# The analytical gradient is: (2, 3, 5, 7, 22*x[4])\n", + "f3_grad_analytical = np.array([2, 3, 5, 7, 22*x[4]])\n", + "\n", + "# Print the analytical gradient:\n", + "print(\"The analytical gradient of f3 is: \", f3_grad_analytical)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that in this case, when sending an array as input argument, the\n", + "output from Autograd is another array. This is the true gradient of\n", + "the function, as opposed to the function in the previous example. By\n", + "using arrays to represent the variables, the output from Autograd\n", + "might be easier to work with, as the output is closer to what one\n", + "could expect form a gradient-evaluting function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "def f4(x):\n", + " return np.sqrt(1+x**2) + np.exp(x) + np.sin(2*np.pi*x)\n", + "\n", + "f4_grad = grad(f4)\n", + "\n", + "x = 2.7\n", + "\n", + "# Print the computed derivative:\n", + "print(\"The computed derivative of f4 at x = %g is: %g\"%(x,f4_grad(x)))\n", + "\n", + "# The analytical derivative is: x/sqrt(1 + x**2) + exp(x) + cos(2*pi*x)*2*pi\n", + "f4_grad_analytical = x/np.sqrt(1 + x**2) + np.exp(x) + np.cos(2*np.pi*x)*2*np.pi\n", + "\n", + "# Print the analytical gradient:\n", + "print(\"The analytical gradient of f4 at x = %g is: %g\"%(x,f4_grad_analytical))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "def f5(x):\n", + " if x >= 0:\n", + " return x**2\n", + " else:\n", + " return -3*x + 1\n", + "\n", + "f5_grad = grad(f5)\n", + "\n", + "x = 2.7\n", + "\n", + "# Print the computed derivative:\n", + "print(\"The computed derivative of f5 at x = %g is: %g\"%(x,f5_grad(x)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "2\n", + "6\n", + " \n", + "<\n", + "<\n", + "<\n", + "!\n", + "!\n", + "C\n", + "O\n", + "D\n", + "E\n", + "_\n", + "B\n", + "L\n", + "O\n", + "C\n", + "K\n", + " \n", + " \n", + "p\n", + "y\n", + "c\n", + "o\n", + "d" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "# Both of the functions are implementation of the sum: sum(x**i) for i = 0, ..., 9\n", + "# The analytical derivative is: sum(i*x**(i-1)) \n", + "f6_grad_analytical = 0\n", + "for i in range(10):\n", + " f6_grad_analytical += i*x**(i-1)\n", + "\n", + "print(\"The analytical derivative of f6 at x = %g is: %g\"%(x,f6_grad_analytical))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "\n", + "def f7(n): # Assume that n is an integer\n", + " if n == 1 or n == 0:\n", + " return 1\n", + " else:\n", + " return n*f7(n-1)\n", + "\n", + "f7_grad = grad(f7)\n", + "\n", + "n = 2.0\n", + "\n", + "print(\"The computed derivative of f7 at n = %d is: %g\"%(n,f7_grad(n)))\n", + "\n", + "# The function f7 is an implementation of the factorial of n.\n", + "# By using the product rule, one can find that the derivative is:\n", + "\n", + "f7_grad_analytical = 0\n", + "for i in range(int(n)-1):\n", + " tmp = 1\n", + " for k in range(int(n)-1):\n", + " if k != i:\n", + " tmp *= (n - k)\n", + " f7_grad_analytical += tmp\n", + "\n", + "print(\"The analytical derivative of f7 at n = %d is: %g\"%(n,f7_grad_analytical))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that if n is equal to zero or one, Autograd will give an error message. This message appears when the output is independent on input.\n", + "\n", + "\n", + "Autograd supports many features. However, there are some functions that is not supported (yet) by Autograd.\n", + "\n", + "Assigning a value to the variable being differentiated with respect to" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "def f8(x): # Assume x is an array\n", + " x[2] = 3\n", + " return x*2\n", + "\n", + "f8_grad = grad(f8)\n", + "\n", + "x = 8.4\n", + "\n", + "print(\"The derivative of f8 is:\",f8_grad(x))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, Autograd tells us that an 'ArrayBox' does not support item assignment. The item assignment is done when the program tries to assign x[2] to the value 3. However, Autograd has implemented the computation of the derivative such that this assignment is not possible." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "def f9(a): # Assume a is an array with 2 elements\n", + " b = np.array([1.0,2.0])\n", + " return a.dot(b)\n", + "\n", + "f9_grad = grad(f9)\n", + "\n", + "x = np.array([1.0,0.0])\n", + "\n", + "print(\"The derivative of f9 is:\",f9_grad(x))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we are told that the 'dot' function does not belong to Autograd's\n", + "version of a Numpy array. To overcome this, an alternative syntax\n", + "which also computed the dot product can be used:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import autograd.numpy as np\n", + "from autograd import grad\n", + "def f9_alternative(x): # Assume a is an array with 2 elements\n", + " b = np.array([1.0,2.0])\n", + " return np.dot(x,b) # The same as x_1*b_1 + x_2*b_2\n", + "\n", + "f9_alternative_grad = grad(f9_alternative)\n", + "\n", + "x = np.array([3.0,0.0])\n", + "\n", + "print(\"The gradient of f9 is:\",f9_alternative_grad(x))\n", + "\n", + "# The analytical gradient of the dot product of vectors x and b with two elements (x_1,x_2) and (b_1, b_2) respectively\n", + "# w.r.t x is (b_1, b_2)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The documentation recommends to avoid inplace operations such as" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "a += b\n", + "a -= b\n", + "a*= b\n", + "a /=b" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "More examples will be added, in particular how to compare autograd with own codes for the gradients." ] } ],