11. Clustering Analysis¶
+12. Clustering Analysis¶
In this chapter we will concern ourselves with the study of cluster analysis. In general terms cluster analysis, or clustering, is the task of grouping a data-set into different distinct categories based on some measure of equality of @@ -343,7 +343,7 @@ exists different categories of clustering algorithms. A few of the most common are: centroid-, distribution-, density- and hierarchical- clustering. We will concern ourselves primarily with the first one.
-
11.1. Basic Idea of the K-means Clustering Algorithm¶
+12.1. Basic Idea of the K-means Clustering Algorithm¶
The simplest of all clustering algorithms is the aptly named k-means algorithm , sometimes also referred to as Lloyds algorithm. It is the simplest and also the most common. From its simplicity it obtains both strengths and weaknesses. @@ -436,7 +436,7 @@ no-where, don’t fret, the method does get a whole lot simpler once we start programming.
-
11.2. The K-means Clustering Algorithm¶
+12.2. The K-means Clustering Algorithm¶
The k-means clustering algorithm goes as follows (note in my opinion this description is a bit complicated and is lifted directly out of ESL HASTIE for deeper understanding purposes)
@@ -460,12 +460,12 @@ understood but to help you along the way we summarize the algorithm as followsThat’s it, nothing magical happening.
-
11.3. Writing Our Own Code¶
+12.3. Writing Our Own Code¶
In the following section we will work to develop a deeper understanding of the previously discussed mathematics through developing codes to do k-means cluster analysis.
-
11.3.1. Basic Python¶
+12.3.1. Basic Python¶
Let us now program the most basic version of the algorithm using nothing but Python with numpy arrays. This code is kept intentionally simple to gradually progress our understanding. There is no vectorization of any kind, and even most @@ -718,7 +718,7 @@ two improvements.
-
- 7. Support Vector Machines, overarching aims
+ 8. Support Vector Machines, overarching aims
@@ -171,12 +171,12 @@
@@ -188,12 +188,12 @@
@@ -205,12 +205,12 @@
@@ -357,39 +357,34 @@
3.7. Housing data, the code
-
-
- 3.8. Splitting our Data in Training and Test data
-
-
- 3.9. Reducing the number of degrees of freedom, overarching view
+ 3.8. Reducing the number of degrees of freedom, overarching view
- 3.10. Testing the Means Squared Error as function of Complexity
+ 3.9. Testing the Means Squared Error as function of Complexity
- 3.11. Exercises
+ 3.10. Exercises
@@ -409,35 +404,34 @@
@@ -893,7 +895,7 @@ a linear \(x\)-dependence we s
-
@@ -1117,7 +1119,7 @@ A
-
-
+
+
+
+
@@ -1878,6 +1884,22 @@ K
\[
\frac{\partial\log{\vert\boldsymbol{A}\vert}}{\partial \boldsymbol{A}}=(\boldsymbol{A}^{-1})^T.
\]
+
-
+
-
+
11.4. Towards a More Numpythonic Code¶
+12.4. Towards a More Numpythonic Code¶
def get_distances_to_clusters(data, centroids):
@@ -1000,8 +1000,8 @@ Runtime: 0.004456043243408203 seconds
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_11_0.png b/doc/LectureNotes/_build/html/_images/chapter1_11_0.png
index 8585b7491..65d331dc5 100644
Binary files a/doc/LectureNotes/_build/html/_images/chapter1_11_0.png and b/doc/LectureNotes/_build/html/_images/chapter1_11_0.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_136_0.png b/doc/LectureNotes/_build/html/_images/chapter1_136_0.png
new file mode 100644
index 000000000..29a1bd652
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter1_136_0.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_13_1.png b/doc/LectureNotes/_build/html/_images/chapter1_13_1.png
index b145ce82a..d70ab26ee 100644
Binary files a/doc/LectureNotes/_build/html/_images/chapter1_13_1.png and b/doc/LectureNotes/_build/html/_images/chapter1_13_1.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_183_1.png b/doc/LectureNotes/_build/html/_images/chapter1_183_1.png
new file mode 100644
index 000000000..0821feaf6
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter1_183_1.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_199_1.png b/doc/LectureNotes/_build/html/_images/chapter1_199_1.png
new file mode 100644
index 000000000..f19e88b16
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter1_199_1.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_201_1.png b/doc/LectureNotes/_build/html/_images/chapter1_201_1.png
new file mode 100644
index 000000000..13fd31349
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter1_201_1.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_203_0.png b/doc/LectureNotes/_build/html/_images/chapter1_203_0.png
new file mode 100644
index 000000000..dfe1e253c
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter1_203_0.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_210_0.png b/doc/LectureNotes/_build/html/_images/chapter1_210_0.png
new file mode 100644
index 000000000..a3eac788b
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter1_210_0.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_218_0.png b/doc/LectureNotes/_build/html/_images/chapter1_218_0.png
new file mode 100644
index 000000000..cbac4fd70
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter1_218_0.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_27_0.png b/doc/LectureNotes/_build/html/_images/chapter1_27_0.png
index bef273ce3..889056124 100644
Binary files a/doc/LectureNotes/_build/html/_images/chapter1_27_0.png and b/doc/LectureNotes/_build/html/_images/chapter1_27_0.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_3_0.png b/doc/LectureNotes/_build/html/_images/chapter1_3_0.png
index 0be2e475d..656724cd5 100644
Binary files a/doc/LectureNotes/_build/html/_images/chapter1_3_0.png and b/doc/LectureNotes/_build/html/_images/chapter1_3_0.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter1_61_10.png b/doc/LectureNotes/_build/html/_images/chapter1_61_10.png
index 16b5dca78..2c196f8b3 100644
Binary files a/doc/LectureNotes/_build/html/_images/chapter1_61_10.png and b/doc/LectureNotes/_build/html/_images/chapter1_61_10.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter2_245_1.png b/doc/LectureNotes/_build/html/_images/chapter2_245_1.png
new file mode 100644
index 000000000..dc3a3ce3e
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter2_245_1.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter2_316_1.png b/doc/LectureNotes/_build/html/_images/chapter2_316_1.png
new file mode 100644
index 000000000..07584f6bb
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter2_316_1.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter3_25_2.png b/doc/LectureNotes/_build/html/_images/chapter3_25_2.png
index b77a1febb..dada82b44 100644
Binary files a/doc/LectureNotes/_build/html/_images/chapter3_25_2.png and b/doc/LectureNotes/_build/html/_images/chapter3_25_2.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter4_114_1.png b/doc/LectureNotes/_build/html/_images/chapter4_114_1.png
new file mode 100644
index 000000000..6c642a6a1
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter4_114_1.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter4_116_0.png b/doc/LectureNotes/_build/html/_images/chapter4_116_0.png
new file mode 100644
index 000000000..dfd44cb86
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter4_116_0.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter4_122_1.png b/doc/LectureNotes/_build/html/_images/chapter4_122_1.png
new file mode 100644
index 000000000..59590fe49
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter4_122_1.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter4_57_2.png b/doc/LectureNotes/_build/html/_images/chapter4_57_2.png
new file mode 100644
index 000000000..fa199fba3
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter4_57_2.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter4_57_3.png b/doc/LectureNotes/_build/html/_images/chapter4_57_3.png
new file mode 100644
index 000000000..88bffd58a
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter4_57_3.png differ
diff --git a/doc/LectureNotes/_build/html/_images/chapter4_57_4.png b/doc/LectureNotes/_build/html/_images/chapter4_57_4.png
new file mode 100644
index 000000000..e242a4058
Binary files /dev/null and b/doc/LectureNotes/_build/html/_images/chapter4_57_4.png differ
diff --git a/doc/LectureNotes/_build/html/_sources/chapter1.ipynb b/doc/LectureNotes/_build/html/_sources/chapter1.ipynb
index c930e0a6b..9dab26cdc 100644
--- a/doc/LectureNotes/_build/html/_sources/chapter1.ipynb
+++ b/doc/LectureNotes/_build/html/_sources/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/_build/html/_sources/chapter2.ipynb b/doc/LectureNotes/_build/html/_sources/chapter2.ipynb
index 651eadf4a..3049b1030 100644
--- a/doc/LectureNotes/_build/html/_sources/chapter2.ipynb
+++ b/doc/LectureNotes/_build/html/_sources/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/_build/html/_sources/chapter3.ipynb b/doc/LectureNotes/_build/html/_sources/chapter3.ipynb
index d7f8222d5..b2933e098 100644
--- a/doc/LectureNotes/_build/html/_sources/chapter3.ipynb
+++ b/doc/LectureNotes/_build/html/_sources/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/_build/html/_sources/chapter4.ipynb b/doc/LectureNotes/_build/html/_sources/chapter4.ipynb
index 7cd60a191..81b9d76c8 100644
--- a/doc/LectureNotes/_build/html/_sources/chapter4.ipynb
+++ b/doc/LectureNotes/_build/html/_sources/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."
]
}
],
diff --git a/doc/LectureNotes/_build/html/chapter1.html b/doc/LectureNotes/_build/html/chapter1.html
index 3e61ffcbd..2e7c0a670 100644
--- a/doc/LectureNotes/_build/html/chapter1.html
+++ b/doc/LectureNotes/_build/html/chapter1.html
@@ -159,7 +159,7 @@
3. Linear Regression¶
3.1. Introduction¶
-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 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.
+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 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.
-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.
3.2. What is Machine Learning?¶
@@ -625,7 +619,15 @@ quality figures. Feel free to explore the extensive gallery 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.
@@ -640,6 +642,7 @@ data with a straight line.
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)
@@ -745,8 +748,7 @@ to be dominated by outliers.
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 their error estimates, or the variance and standard deviation and many @@ -787,13 +789,13 @@ example of the functionality of Scikit-Learn.
The intercept alpha:
- [1.84860939]
+ [1.89992521]
Coefficient beta :
- [[5.2571699]]
-Mean squared error: 0.25
-Variance score: 0.90
+ [[5.12908963]]
+Mean squared error: 0.20
+Variance score: 0.92
Mean squared log error: 0.01
-Mean absolute error: 0.40
+Mean absolute error: 0.36
@@ -893,7 +895,7 @@ a linear \(x\)-dependence we s
-0.0050000000000000044
+0.005000000000000007
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}\). -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.# Now we set up the design matrix X
@@ -1131,8 +1133,11 @@ It has dimensionality \(p\times n\)
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.
-With scikitlearn we are now ready to use linear regression and fit our data.
+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 Scikit-Learn we are now ready to use linear regression and fit our data.
clf = skl.LinearRegression().fit(X, Energies)
@@ -1176,7 +1181,7 @@ Mean absolute error: 0.05
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!
#Decision Tree Regression
@@ -1228,11 +1233,20 @@ A
270 3344 160 110 270 Ds 7.253775 7.253775
[267 rows x 6 columns]
-0.009883615646716186
+0.009883615646716188
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,14 +1284,6 @@ functionality.
-
-
+
+
-
-
-
-
@@ -1360,7 +1362,11 @@ functionality.
warnings.warn(
-/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
- warnings.warn(
-/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
- warnings.warn(
-/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
- warnings.warn(
-/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
@@ -1334,16 +1340,12 @@ functionality.
warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
warnings.warn(
+/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
+ warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
warnings.warn(
-/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
- warnings.warn(
-/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
- warnings.warn(
/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
warnings.warn(
+/Users/mhjensen/opt/anaconda3/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:582: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (100) reached and the optimization hasn't converged yet.
+ warnings.warn(
+
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
+
+\[
+\frac{\partial^2 C(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}^T\partial \boldsymbol{\beta}} =\frac{2}{n}\boldsymbol{X}^T\boldsymbol{X}.
+\]
+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
+
+\[
+\boldsymbol{H}=\boldsymbol{X}^T\boldsymbol{X}.
+\]
+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 \(\boldsymbol{\epsilon}\) are in turn given by
\[
@@ -1935,7 +1957,7 @@ write
-
+
+
We can easily test our fit by computing the \(R2\) score that we discussed in connection with the functionality of Scikit-Learn in the introductory slides. @@ -2138,9 +2160,7 @@ time the file is presented as a standard csv file.
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.# Common imports
@@ -2203,18 +2223,6 @@ hyperparameter \(\lambda\), al
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}]$')
@@ -2223,8 +2231,6 @@ hyperparameter \(\lambda\), al
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()
@@ -2236,20 +2242,13 @@ hyperparameter \(\lambda\), al
Variance score: 1.00
Mean absolute error: 2.83
[ 0. 618.32047562 -861.13519106 1404.91549644] -11.057088709963637
-Mean squared error: 197.93
-Variance score: 1.00
-Mean absolute error: 11.69
-[ 0. 28.18220995 282.79902342 842.30879705] 12.946893955211749
+
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.
@@ -2264,6 +2263,95 @@ 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.
+
@@ -2460,7 +2548,7 @@ dtype: int64
warnings.warn(msg, FutureWarning)
-
+
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.
+
+
+
+
+
+
+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))
+
+
+
+
+[ 2.00247083e+00 -3.54710530e-03 5.02096350e+00]
+Training R2
+0.99620088364924
+Training MSE
+0.008171076530362356
+Test R2
+0.9968966158779216
+Test MSE
+0.006530396683779109
+Alternatively, you could write your own test-train splitting function as shown here.
+
+
+
+
+
+
+# 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
+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.
import os
@@ -2335,13 +2423,13 @@ but now splitting the data into a training set and a test set.
Training R2
-0.9999869956119286
+0.9999850282256434
Training MSE
-5.745136489050356
+6.653725417896576
Test R2
-0.9999787681537219
+0.9999887274726137
Test MSE
-9.712199063818309
+5.071447584032141
+
It is now useful to look at the correlation matrix
@@ -2478,7 +2566,7 @@ dtype: int64<AxesSubplot:>
+
From the above coorelation plot we can see that MEDV is strongly correlated to LSTAT and RM. We see also that RAD and TAX are stronly correlated, but we don’t include this in our features together to avoid multi-colinearity
@@ -2501,7 +2589,7 @@ dtype: int64
+
+
+
+
+
+
+
+