typos
This commit is contained in:
Morten Hjorth-Jensen
@@ -67,3 +67,13 @@ found info about 5 exercises
|
||||
|
||||
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
|
||||
output in chapter1.ipynb
|
||||
Translating doconce text in chapter1.do.txt to ipynb
|
||||
*** replacing \bm{...} by \boldsymbol{...} (\bm is not supported by MathJax)
|
||||
found info about 5 exercises
|
||||
|
||||
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
|
||||
|
||||
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
|
||||
|
||||
*** warning: latex envir \begin{bmatrix} does not work well in Markdown. Stick to \[ ... \], equation, equation*, align, or align* environments in math environments.
|
||||
output in chapter1.ipynb
|
||||
|
||||
@@ -88,7 +88,7 @@ role when we develop a specific machine learning algorithm.
|
||||
Machine learning is an extremely rich field, in spite of its young
|
||||
age. The increases we have seen during the last three decades in
|
||||
computational capabilities have been followed by developments of
|
||||
methods and techniques for analyzing and handling large date sets,
|
||||
methods and techniques for analyzing and handling large data sets,
|
||||
relying heavily on statistics, computer science and mathematics. The
|
||||
field is rather new and developing rapidly. Popular software packages
|
||||
written in Python for machine learning like
|
||||
@@ -113,7 +113,7 @@ two main categories. In *supervised learning* we know the answer to a
|
||||
problem, and let the computer deduce the logic behind it. On the other
|
||||
hand, *unsupervised learning* is a method for finding patterns and
|
||||
relationship in data sets without any prior knowledge of the system.
|
||||
Some authours also operate with a third category, namely
|
||||
Some authors also operate with a third category, namely
|
||||
*reinforcement learning*. This is a paradigm of learning inspired by
|
||||
behavioral psychology, where learning is achieved by trial-and-error,
|
||||
solely from rewards and punishment.
|
||||
@@ -176,14 +176,14 @@ what is the likelihood of finding $B$.
|
||||
In science and engineering we often end up in situations where we want to infer (or learn) a
|
||||
quantitative model $M$ for a given set of sample points $\bm{X} \in [x_1, x_2,\dots x_N]$.
|
||||
|
||||
As we will see repeatedely in these lectures, we could try to fit these data points to a model given by a
|
||||
As we will see repeatedly in these lectures, we could try to fit these data points to a model given by a
|
||||
straight line, or if we wish to be more sophisticated to a more complex
|
||||
function.
|
||||
|
||||
The reason for inferring such a model is that it
|
||||
serves many useful purposes. On the one hand, the model can reveal information
|
||||
encoded in the data or underlying mechanisms from which the data were generated. For instance, we could discover important
|
||||
corelations that relate interesting physics interpretations.
|
||||
correlations that relate interesting physics interpretations.
|
||||
|
||||
In addition, it can simplify the representation of the given data set and help
|
||||
us in making predictions about future data samples.
|
||||
@@ -328,7 +328,7 @@ y = 10x+0.01 \times N(0,1),
|
||||
where $x$ is defined as before. Does the fit look better? Indeed, by
|
||||
reducing the role of the noise given by the normal distribution we see immediately that
|
||||
our linear prediction seemingly reproduces better the training
|
||||
set. However, this testing 'by the eye' is obviouly not satisfactory in the
|
||||
set. However, this testing 'by the eye' is obviously not satisfactory in the
|
||||
long run. Here we have only defined the training data and our model, and
|
||||
have not discussed a more rigorous approach to the _cost_ function.
|
||||
|
||||
|
||||
@@ -104,7 +104,7 @@ later shrinkage methods like Ridge and Lasso regressions.
|
||||
This is given by the _Singular Value Decomposition_ (SVD) algorithm,
|
||||
perhaps the most powerful linear algebra algorithm. The SVD provides
|
||||
a numerically stable matrix decomposition that is used in a large
|
||||
swath oc applications and the decomposition is always stable
|
||||
swath of applications and the decomposition is always stable
|
||||
numerically.
|
||||
|
||||
In machine learning it plays a central role in dealing with for
|
||||
@@ -123,7 +123,7 @@ when the matrix $\bm{X}$ (our so-called design matrix) is high-dimensional,
|
||||
are problems with near singular or singular matrices. The column vectors of $\bm{X}$
|
||||
may be linearly dependent, normally referred to as super-collinearity.
|
||||
This means that the matrix may be rank deficient and it is basically impossible to
|
||||
to model the data using linear regression. As an example, consider the matrix
|
||||
model the data using linear regression. As an example, consider the matrix
|
||||
!bt
|
||||
\begin{align*}
|
||||
\mathbf{X} & = \left[
|
||||
@@ -143,7 +143,7 @@ The columns of $\bm{X}$ are linearly dependent. We see this easily since the
|
||||
the first column is the row-wise sum of the other two columns. The rank (more correct,
|
||||
the column rank) of a matrix is the dimension of the space spanned by the
|
||||
column vectors. Hence, the rank of $\mathbf{X}$ is equal to the number
|
||||
of linearly independent columns. In this particular case the matrix has rank 2.
|
||||
of linearly independent columns. In this particular case the matrix has rank 1.
|
||||
|
||||
Super-collinearity of an $(n \times p)$-dimensional design matrix $\mathbf{X}$ implies
|
||||
that the inverse of the matrix $\bm{X}^T\bm{X}$ (the matrix we need to invert to solve the linear regression equations) is non-invertible. If we have a square matrix that does not have an inverse, we say this matrix singular. The example here demonstrates this
|
||||
@@ -170,7 +170,7 @@ If our design matrix $\bm{X}$ which enters the linear regression problem
|
||||
!et
|
||||
has linearly dependent column vectors, we will not be able to compute the inverse
|
||||
of $\bm{X}^T\bm{X}$ and we cannot find the parameters (estimators) $\beta_i$.
|
||||
The estimators are only well-defined if $(\bm{X}^{T}\bm{X})^{-1}$ exits.
|
||||
The estimators are only well-defined if $(\bm{X}^{T}\bm{X})$ can be inverted.
|
||||
This is more likely to happen when the matrix $\bm{X}$ is high-dimensional. In this case it is likely to encounter a situation where
|
||||
the regression parameters $\beta_i$ cannot be estimated.
|
||||
|
||||
@@ -188,7 +188,7 @@ where $\bm{I}$ is the identity matrix. When we discuss _Ridge_ regression this
|
||||
===== Basic math of the SVD =====
|
||||
|
||||
|
||||
From standard linear algebra we know that a square matrix $\bm{X}$ can be diagonalized if and only it is
|
||||
From standard linear algebra we know that a square matrix $\bm{X}$ can be diagonalized if and only if it is
|
||||
a so-called "normal matrix":"https://en.wikipedia.org/wiki/Normal_matrix", that is if $\bm{X}\in {\mathbb{R}}^{n\times n}$
|
||||
we have $\bm{X}\bm{X}^T=\bm{X}^T\bm{X}$ or if $\bm{X}\in {\mathbb{C}}^{n\times n}$ we have $\bm{X}\bm{X}^{\dagger}=\bm{X}^{\dagger}\bm{X}$.
|
||||
The matrix has then a set of eigenpairs
|
||||
@@ -381,7 +381,6 @@ def SVDinv(A):
|
||||
return np.matmul(V,np.matmul(invD,UT))
|
||||
|
||||
|
||||
#X = np.array([ [1.0, -1.0, 2.0], [1.0, 0.0, 1.0], [1.0, 2.0, -1.0], [1.0, 1.0, 0.0] ])
|
||||
# Non-singular square matrix
|
||||
X = np.array( [ [1,2,3],[2,4,5],[3,5,6]])
|
||||
print(X)
|
||||
@@ -402,7 +401,7 @@ rectangular matrices where the number of rows and columns are not equal.
|
||||
It is also called the the Moore-Penrose Inverse after two independent discoverers of the method or the Generalized Inverse.
|
||||
It is used for the calculation of the inverse for singular or near singular matrices and for rectangular matrices.
|
||||
|
||||
Using the SVD we can obtain the pseudoinverse of a matrix $\bm{A}$ (labeled here as $\bm{A}_{\mathrm{PI}}$
|
||||
Using the SVD we can obtain the pseudoinverse (PI) of a matrix $\bm{A}$ (labeled here as $\bm{A}_{\mathrm{PI}}$
|
||||
!bt
|
||||
\[
|
||||
\bm{A}_{\mathrm{PI}}= \bm{V}\bm{D}_{\mathrm{PI}}\bm{U}^T,
|
||||
@@ -469,7 +468,7 @@ We can SVD decompose our matrix as
|
||||
!et
|
||||
where $\bm{U}$ is an orthogonal matrix of dimension $n\times n$, meaning that $\bm{U}\bm{U}^T=\bm{U}^T\bm{U}=\bm{I}_n$. Here $\bm{I}_n$ is the unit matrix of dimension $n \times n$.
|
||||
|
||||
Similarly, $\bm{V}$ is an orthogonal matrix of dimension $p\times p$, meaning that $\bm{V}\bm{V}^T=\bm{V}^T\bm{V}=\bm{I}_p$. Here $\bm{I}_n$ is the unit matrix of dimension $p \times p$.
|
||||
Similarly, $\bm{V}$ is an orthogonal matrix of dimension $p\times p$, meaning that $\bm{V}\bm{V}^T=\bm{V}^T\bm{V}=\bm{I}_p$. Here $\bm{I}_p$ is the unit matrix of dimension $p \times p$.
|
||||
|
||||
Finally $\bm{\Sigma}$ contains the singular values $\sigma_i$. This matrix has dimension $n\times p$ and the singular values $\sigma_i$ are all positive. The non-zero values are ordered in descending order, that is
|
||||
|
||||
@@ -656,7 +655,7 @@ function, that is we have
|
||||
\frac{\partial^2 C(\bm{\beta})}{\partial \bm{\beta}^T\partial \bm{\beta}} =\frac{2}{n}\bm{X}^T\bm{X}.
|
||||
\]
|
||||
!et
|
||||
This quantity defines was what is called the Hessian matrix (the second derivative of a function we want to optimize).
|
||||
This quantity defines what is called the Hessian matrix (the second derivative of the cost function we want to optimize).
|
||||
|
||||
The Hessian matrix plays an important role and is defined in this course as
|
||||
|
||||
@@ -777,7 +776,7 @@ with a given vector
|
||||
!et
|
||||
|
||||
With these definitions, we can now rewrite our $2\times 2$
|
||||
correlation/covariance matrix in terms of a moe general design/feature
|
||||
correlation/covariance matrix in terms of a more general design/feature
|
||||
matrix $\bm{X}\in {\mathbb{R}}^{n\times p}$. This leads to a $p\times p$
|
||||
covariance matrix for the vectors $\bm{x}_i$ with $i=0,1,\dots,p-1$
|
||||
|
||||
@@ -883,7 +882,7 @@ The above procedure with _numpy_ can be made more compact if we use _pandas_.
|
||||
|
||||
|
||||
|
||||
We whow here how we can set up the correlation matrix using _pandas_, as done in this simple code
|
||||
We know here how we can set up the correlation matrix using _pandas_, as done in this simple code
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
@@ -1344,7 +1343,7 @@ This equation does not lead to a nice analytical equation as in Ridge regression
|
||||
Let us assume that our design matrix is given by unit (identity) matrix, that is a square diagonal matrix with ones only along the
|
||||
diagonal. In this case we have an equal number of rows and columns $n=p$.
|
||||
|
||||
Our model approximation is just $\tilde{\bm{y}}=\bm{\beta}$ and the mean squared error and thereby the cost function for ordinary least sqquares (OLS) is then (we drop the term $1/n$)
|
||||
Our model approximation is just $\tilde{\bm{y}}=\bm{\beta}$ and the mean squared error and thereby the cost function for ordinary least squares (OLS) is then (we drop the term $1/n$)
|
||||
!bt
|
||||
\[
|
||||
C(\bm{\beta})=\sum_{i=0}^{p-1}(y_i-\beta_i)^2,
|
||||
@@ -1396,7 +1395,7 @@ which leads to
|
||||
Plotting these results ("figure in handwritten notes for week 36":"https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2021/NotesSeptember9.pdf") shows clearly that Lasso regression suppresses (sets to zero) values of $\beta_i$ for specific values of $\lambda$. Ridge regression reduces on the other hand the values of $\beta_i$ as function of $\lambda$.
|
||||
|
||||
|
||||
As another examples,
|
||||
As another example,
|
||||
let us assume we have a data set with outputs/targets given by the vector
|
||||
|
||||
!bt
|
||||
@@ -1600,13 +1599,13 @@ plt.show()
|
||||
|
||||
!ec
|
||||
|
||||
We see here that we reach a plateau for the Ridge results. Writing out the coefficients $\bm{\beta}$, we that they are getting smaller and smaller and our error stabilizes since the predicted values of $\tilde{\bm{y}}$ approach zero.
|
||||
We see here that we reach a plateau for the Ridge results. Writing out the coefficients $\bm{\beta}$, we observe that they are getting smaller and smaller and our error stabilizes since the predicted values of $\tilde{\bm{y}}$ approach zero.
|
||||
|
||||
This happens also for Lasso regression, as seen from the next code
|
||||
output. The difference is that Lasso shrinks the values of $\beta$ to
|
||||
zero at a much earlier stage and the results flatten out. We see that
|
||||
Lasso gives also an excellent fit for small values of $\lambda$ and
|
||||
shows rthe best performance of the three regression methods.
|
||||
shows the best performance of the three regression methods.
|
||||
|
||||
!bc pycod
|
||||
import os
|
||||
@@ -1967,7 +1966,7 @@ p(y_i, \bm{X}\vert\bm{\beta})=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left[-\frac{(y_
|
||||
!et
|
||||
which reads as finding the likelihood of an event $y_i$ with the input variables $\bm{X}$ given the parameters (to be determined) $\bm{\beta}$.
|
||||
|
||||
Since these events are assumed to be independent and identicall distributed we can build the probability distribution function (PDF) for all possible event $\bm{y}$ as the product of the single events, that is we have
|
||||
Since these events are assumed to be independent and identically distributed we can build the probability distribution function (PDF) for all possible event $\bm{y}$ as the product of the single events, that is we have
|
||||
|
||||
!bt
|
||||
\[
|
||||
@@ -2072,7 +2071,7 @@ p(X \cup Y)= p(X)+p(Y)-p(X \cap Y).
|
||||
The product rule (aka joint probability) is given by
|
||||
!bt
|
||||
\[
|
||||
p(X \cup Y)= p(X,Y)= p(X\vert Y)p(Y)=p(Y\vert X)p(X),
|
||||
p(X \cap Y)= p(X,Y)= p(X\vert Y)p(Y)=p(Y\vert X)p(X),
|
||||
\]
|
||||
!et
|
||||
where we read $p(X\vert Y)$ as the likelihood of obtaining $X$ given $Y$.
|
||||
@@ -2352,7 +2351,7 @@ parameters $\beta_j$ as function of polynomial order and of the added
|
||||
noise. Here we recommend to use $\sigma^2=1$ as variance for the
|
||||
added noise (which follows a normal distribution with mean value zero).
|
||||
Comment your results. If you have a large noise term, do the parameters $\beta_j$ vary more as function
|
||||
model complexity? And what about their variance?
|
||||
of model complexity? And what about their variance?
|
||||
|
||||
|
||||
|
||||
@@ -2393,7 +2392,7 @@ p(\bm{\beta}\vert\bm{D})\propto p(\bm{D}\vert\bm{\beta})p(\bm{\beta}).
|
||||
\]
|
||||
!et
|
||||
|
||||
We have a model for $p(\bm{D}\vert\bm{\beta})$ but need one for the _prior_ $p(\bm{\beta}$!
|
||||
We have a model for $p(\bm{D}\vert\bm{\beta})$ but need one for the _prior_ $p(\bm{\beta})$!
|
||||
|
||||
|
||||
|
||||
@@ -2464,14 +2463,14 @@ constants terms that do not depend on $\beta$, we have
|
||||
|
||||
!bt
|
||||
\[
|
||||
C(\bm{\beta}=\frac{\vert\vert (\bm{y}-\bm{X}\bm{\beta})\vert\vert_2^2}{2\sigma^2}+\frac{1}{\tau}\vert\vert\bm{\beta}\vert\vert_1,
|
||||
C(\bm{\beta})=\frac{\vert\vert (\bm{y}-\bm{X}\bm{\beta})\vert\vert_2^2}{2\sigma^2}+\frac{1}{\tau}\vert\vert\bm{\beta}\vert\vert_1,
|
||||
\]
|
||||
!et
|
||||
and replacing $1/\tau$ with $\lambda$ we have
|
||||
|
||||
!bt
|
||||
\[
|
||||
C(\bm{\beta}=\frac{\vert\vert (\bm{y}-\bm{X}\bm{\beta})\vert\vert_2^2}{2\sigma^2}+\lambda\vert\vert\bm{\beta}\vert\vert_1,
|
||||
C(\bm{\beta})=\frac{\vert\vert (\bm{y}-\bm{X}\bm{\beta})\vert\vert_2^2}{2\sigma^2}+\lambda\vert\vert\bm{\beta}\vert\vert_1,
|
||||
\]
|
||||
!et
|
||||
which is our Lasso cost function!
|
||||
|
||||
@@ -280,9 +280,19 @@ The convex subsets of $\mathbb{R}$ are the intervals of
|
||||
$\mathbb{R}$. Examples of convex sets of $\mathbb{R}^2$ are the
|
||||
regular polygons (triangles, rectangles, pentagons, etc...).
|
||||
|
||||
_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]$.
|
||||
|
||||
|
||||
_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.
|
||||
If $\leq$ is replaced with a strict inequality 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 discussed below.
|
||||
|
||||
|
||||
In the following we state first and second-order conditions which
|
||||
@@ -296,7 +306,7 @@ all $x$ in the domain of $f$). Then $f$ is convex if and only if $D_f$
|
||||
is a convex set and $$f(y) \geq f(x) + \nabla f(x)^T (y-x) $$ holds
|
||||
for all $x,y \in D_f$. This condition means that for a convex function
|
||||
the first order Taylor expansion (right hand side above) at any point
|
||||
a global under estimator of the function. To convince yourself you can
|
||||
is a global under estimator of the function. To convince yourself you can
|
||||
make a drawing of $f(x) = x^2+1$ and draw the tangent line to $f(x)$ and
|
||||
note that it is always below the graph.
|
||||
!eblock
|
||||
@@ -1756,7 +1766,14 @@ a /=b
|
||||
|
||||
|
||||
|
||||
===== Using Autograd with OLS =====
|
||||
!split
|
||||
===== Replace or not =====
|
||||
|
||||
In the above code, we have use replacement in setting up the
|
||||
mini-batches. The discussion
|
||||
"here":"https://sebastianraschka.com/faq/docs/sgd-methods.html" may be
|
||||
useful.
|
||||
===== Using Autograd =====
|
||||
|
||||
We conclude the part on optmization by showing how we can make codes
|
||||
for linear regression and logistic regression using _autograd_. The
|
||||
@@ -1816,8 +1833,110 @@ plt.show()
|
||||
!ec
|
||||
|
||||
|
||||
===== Same code but now with momentum gradient descent =====
|
||||
!bc pycod
|
||||
# Using Autograd to calculate gradients for OLS
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import autograd.numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from autograd import grad
|
||||
|
||||
=== Including Stochastic Gradient Descent with Autograd ===
|
||||
def CostOLS(beta):
|
||||
return (1.0/n)*np.sum((y-X @ beta)**2)
|
||||
|
||||
n = 100
|
||||
x = 2*np.random.rand(n,1)
|
||||
y = 4+3*x#+np.random.randn(n,1)
|
||||
|
||||
X = np.c_[np.ones((n,1)), x]
|
||||
XT_X = X.T @ X
|
||||
theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
|
||||
print("Own inversion")
|
||||
print(theta_linreg)
|
||||
# Hessian matrix
|
||||
H = (2.0/n)* XT_X
|
||||
EigValues, EigVectors = np.linalg.eig(H)
|
||||
print(f"Eigenvalues of Hessian Matrix:{EigValues}")
|
||||
|
||||
theta = np.random.randn(2,1)
|
||||
eta = 1.0/np.max(EigValues)
|
||||
Niterations = 30
|
||||
|
||||
# define the gradient
|
||||
training_gradient = grad(CostOLS)
|
||||
|
||||
for iter in range(Niterations):
|
||||
gradients = training_gradient(theta)
|
||||
theta -= eta*gradients
|
||||
print(iter,gradients[0],gradients[1])
|
||||
print("theta from own gd")
|
||||
print(theta)
|
||||
|
||||
# Now improve with momentum gradient descent
|
||||
change = 0.0
|
||||
delta_momentum = 0.3
|
||||
for iter in range(Niterations):
|
||||
# calculate gradient
|
||||
gradients = training_gradient(theta)
|
||||
# calculate update
|
||||
new_change = eta*gradients+delta_momentum*change
|
||||
# take a step
|
||||
theta -= new_change
|
||||
# save the change
|
||||
change = new_change
|
||||
print(iter,gradients[0],gradients[1])
|
||||
print("theta from own gd wth momentum")
|
||||
print(theta)
|
||||
|
||||
!ec
|
||||
We note indeed a considerable increase in efficiency here, we less iterations needed.
|
||||
However, if we can invert the Hessian matrix, this is the preferred approach, as shown in the example here.
|
||||
|
||||
!bc pycod
|
||||
# Using Newton's method
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import autograd.numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from autograd import grad
|
||||
|
||||
def CostOLS(beta):
|
||||
return (1.0/n)*np.sum((y-X @ beta)**2)
|
||||
|
||||
n = 100
|
||||
x = 2*np.random.rand(n,1)
|
||||
y = 4+3*x+np.random.randn(n,1)
|
||||
|
||||
X = np.c_[np.ones((n,1)), x]
|
||||
XT_X = X.T @ X
|
||||
beta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
|
||||
print("Own inversion")
|
||||
print(beta_linreg)
|
||||
# Hessian matrix
|
||||
H = (2.0/n)* XT_X
|
||||
# Note that here the Hessian does not depend on the parameters beta
|
||||
invH = np.linalg.pinv(H)
|
||||
EigValues, EigVectors = np.linalg.eig(H)
|
||||
print(f"Eigenvalues of Hessian Matrix:{EigValues}")
|
||||
|
||||
beta = np.random.randn(2,1)
|
||||
Niterations = 5
|
||||
|
||||
# define the gradient
|
||||
training_gradient = grad(CostOLS)
|
||||
|
||||
for iter in range(Niterations):
|
||||
gradients = training_gradient(beta)
|
||||
beta -= invH @ gradients
|
||||
print(iter,gradients[0],gradients[1])
|
||||
print("beta from own Newton code")
|
||||
print(beta)
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
===== Including Stochastic Gradient Descent with Autograd =====
|
||||
In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using _autograd_.
|
||||
|
||||
!bc pycod
|
||||
@@ -1898,43 +2017,223 @@ print(theta)
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
=== And Logistic Regression ===
|
||||
|
||||
Here we include momentum in the standard gradient descent approach.
|
||||
!bc pycod
|
||||
# Using Autograd to calculate gradients using SGD
|
||||
# OLS example
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import autograd.numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from autograd import grad
|
||||
|
||||
def sigmoid(x):
|
||||
return 0.5 * (np.tanh(x / 2.) + 1)
|
||||
# Note change from previous example
|
||||
def CostOLS(y,X,theta):
|
||||
return np.sum((y-X @ theta)**2)
|
||||
|
||||
def logistic_predictions(weights, inputs):
|
||||
# Outputs probability of a label being true according to logistic model.
|
||||
return sigmoid(np.dot(inputs, weights))
|
||||
n = 100
|
||||
x = 2*np.random.rand(n,1)
|
||||
y = 4+3*x+np.random.randn(n,1)
|
||||
|
||||
def training_loss(weights):
|
||||
# Training loss is the negative log-likelihood of the training labels.
|
||||
preds = logistic_predictions(weights, inputs)
|
||||
label_probabilities = preds * targets + (1 - preds) * (1 - targets)
|
||||
return -np.sum(np.log(label_probabilities))
|
||||
X = np.c_[np.ones((n,1)), x]
|
||||
XT_X = X.T @ X
|
||||
theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
|
||||
print("Own inversion")
|
||||
print(theta_linreg)
|
||||
# Hessian matrix
|
||||
H = (2.0/n)* XT_X
|
||||
EigValues, EigVectors = np.linalg.eig(H)
|
||||
print(f"Eigenvalues of Hessian Matrix:{EigValues}")
|
||||
|
||||
# Build a toy dataset.
|
||||
inputs = np.array([[0.52, 1.12, 0.77],
|
||||
[0.88, -1.08, 0.15],
|
||||
[0.52, 0.06, -1.30],
|
||||
[0.74, -2.49, 1.39]])
|
||||
targets = np.array([True, True, False, True])
|
||||
theta = np.random.randn(2,1)
|
||||
eta = 1.0/np.max(EigValues)
|
||||
Niterations = 100
|
||||
|
||||
# Define a function that returns gradients of training loss using Autograd.
|
||||
training_gradient_fun = grad(training_loss)
|
||||
# Note that we request the derivative wrt third argument (theta, 2 here)
|
||||
training_gradient = grad(CostOLS,2)
|
||||
|
||||
# Optimize weights using gradient descent.
|
||||
weights = np.array([0.0, 0.0, 0.0])
|
||||
print("Initial loss:", training_loss(weights))
|
||||
for i in range(100):
|
||||
weights -= training_gradient_fun(weights) * 0.01
|
||||
for iter in range(Niterations):
|
||||
gradients = (1.0/n)*training_gradient(y, X, theta)
|
||||
theta -= eta*gradients
|
||||
print("theta from own gd")
|
||||
print(theta)
|
||||
|
||||
print("Trained loss:", training_loss(weights))
|
||||
|
||||
n_epochs = 50
|
||||
M = 5 #size of each minibatch
|
||||
m = int(n/M) #number of minibatches
|
||||
t0, t1 = 5, 50
|
||||
def learning_schedule(t):
|
||||
return t0/(t+t1)
|
||||
|
||||
theta = np.random.randn(2,1)
|
||||
|
||||
change = 0.0
|
||||
delta_momentum = 0.3
|
||||
|
||||
for epoch in range(n_epochs):
|
||||
for i in range(m):
|
||||
random_index = M*np.random.randint(m)
|
||||
xi = X[random_index:random_index+M]
|
||||
yi = y[random_index:random_index+M]
|
||||
gradients = (1.0/M)*training_gradient(yi, xi, theta)
|
||||
eta = learning_schedule(epoch*m+i)
|
||||
# calculate update
|
||||
new_change = eta*gradients+delta_momentum*change
|
||||
# take a step
|
||||
theta -= new_change
|
||||
# save the change
|
||||
change = new_change
|
||||
print("theta from own sdg with momentum")
|
||||
print(theta)
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
=== Similar (second order function now) problem but now with AdaGrad ===
|
||||
!bc pycod
|
||||
# Using Autograd to calculate gradients using AdaGrad and Stochastic Gradient descent
|
||||
# OLS example
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import autograd.numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from autograd import grad
|
||||
|
||||
# Note change from previous example
|
||||
def CostOLS(y,X,theta):
|
||||
return np.sum((y-X @ theta)**2)
|
||||
|
||||
n = 10000
|
||||
x = np.random.rand(n,1)
|
||||
y = 2.0+3*x +4*x*x# +np.random.randn(n,1)
|
||||
|
||||
X = np.c_[np.ones((n,1)), x, x*x]
|
||||
XT_X = X.T @ X
|
||||
theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
|
||||
print("Own inversion")
|
||||
print(theta_linreg)
|
||||
|
||||
|
||||
# Note that we request the derivative wrt third argument (theta, 2 here)
|
||||
training_gradient = grad(CostOLS,2)
|
||||
# Define parameters for Stochastic Gradient Descent
|
||||
n_epochs = 50
|
||||
M = 5 #size of each minibatch
|
||||
m = int(n/M) #number of minibatches
|
||||
# Guess for unknown parameters theta
|
||||
theta = np.random.randn(3,1)
|
||||
|
||||
# Value for learning rate
|
||||
eta = 0.01
|
||||
# Including AdaGrad parameter to avoid possible division by zero
|
||||
delta = 1e-8
|
||||
for epoch in range(n_epochs):
|
||||
# The outer product is calculated from scratch for each epoch
|
||||
Giter = np.zeros(shape=(3,3))
|
||||
for i in range(m):
|
||||
random_index = M*np.random.randint(m)
|
||||
xi = X[random_index:random_index+M]
|
||||
yi = y[random_index:random_index+M]
|
||||
gradients = (1.0/M)*training_gradient(yi, xi, theta)
|
||||
# Calculate the outer product of the gradients
|
||||
Giter +=gradients @ gradients.T
|
||||
# Simpler algorithm with only diagonal elements
|
||||
Ginverse = np.c_[eta/(delta+np.sqrt(np.diagonal(Giter)))]
|
||||
# compute update
|
||||
update = np.multiply(Ginverse,gradients)
|
||||
theta -= update
|
||||
print("theta from own AdaGrad")
|
||||
print(theta)
|
||||
!ec
|
||||
|
||||
Running this code we note an almost perfect agreement with the results from matrix inversion.
|
||||
|
||||
Similarly, here is our implementation of RMSprop.
|
||||
!bc pycod
|
||||
# Using Autograd to calculate gradients using RMSprop and Stochastic Gradient descent
|
||||
# OLS example
|
||||
from random import random, seed
|
||||
import numpy as np
|
||||
import autograd.numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from autograd import grad
|
||||
|
||||
# Note change from previous example
|
||||
def CostOLS(y,X,theta):
|
||||
return np.sum((y-X @ theta)**2)
|
||||
|
||||
n = 10000
|
||||
x = np.random.rand(n,1)
|
||||
y = 2.0+3*x +4*x*x# +np.random.randn(n,1)
|
||||
|
||||
X = np.c_[np.ones((n,1)), x, x*x]
|
||||
XT_X = X.T @ X
|
||||
theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
|
||||
print("Own inversion")
|
||||
print(theta_linreg)
|
||||
|
||||
|
||||
# Note that we request the derivative wrt third argument (theta, 2 here)
|
||||
training_gradient = grad(CostOLS,2)
|
||||
# Define parameters for Stochastic Gradient Descent
|
||||
n_epochs = 50
|
||||
M = 5 #size of each minibatch
|
||||
m = int(n/M) #number of minibatches
|
||||
# Guess for unknown parameters theta
|
||||
theta = np.random.randn(3,1)
|
||||
|
||||
# Value for learning rate
|
||||
eta = 0.01
|
||||
# Value for parameter rho
|
||||
rho = 0.99
|
||||
# Including AdaGrad parameter to avoid possible division by zero
|
||||
delta = 1e-8
|
||||
for epoch in range(n_epochs):
|
||||
Giter = np.zeros(shape=(3,3))
|
||||
for i in range(m):
|
||||
random_index = M*np.random.randint(m)
|
||||
xi = X[random_index:random_index+M]
|
||||
yi = y[random_index:random_index+M]
|
||||
gradients = (1.0/M)*training_gradient(yi, xi, theta)
|
||||
# Previous value for the outer product of gradients
|
||||
Previous = Giter
|
||||
# Accumulated gradient
|
||||
Giter +=gradients @ gradients.T
|
||||
# Scaling with rho the new and the previous results
|
||||
Gnew = (rho*Previous+(1-rho)*Giter)
|
||||
# Taking the diagonal only and inverting
|
||||
Ginverse = np.c_[eta/(delta+np.sqrt(np.diagonal(Gnew)))]
|
||||
# Hadamard product
|
||||
update = np.multiply(Ginverse,gradients)
|
||||
theta -= update
|
||||
print("theta from own RMSprop")
|
||||
print(theta)
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
|
||||
===== Introducing "JAX":"https://jax.readthedocs.io/en/latest/" =====
|
||||
|
||||
Presently, instead of using _autograd_, we recommend using "JAX":"https://jax.readthedocs.io/en/latest/"
|
||||
|
||||
_JAX_ is Autograd and "XLA (Accelerated Linear Algebra))":"https://www.tensorflow.org/xla",
|
||||
brought together for high-performance numerical computing and machine learning research.
|
||||
It provides composable transformations of Python+NumPy programs: differentiate, vectorize, parallelize, Just-In-Time compile to GPU/TPU, and more.
|
||||
|
||||
Here's a simple example on how you can use _JAX_ to compute the derivate of the logistic function.
|
||||
|
||||
!bc pycod
|
||||
import jax.numpy as jnp
|
||||
from jax import grad, jit, vmap
|
||||
|
||||
def sum_logistic(x):
|
||||
return jnp.sum(1.0 / (1.0 + jnp.exp(-x)))
|
||||
|
||||
x_small = jnp.arange(3.)
|
||||
derivative_fn = grad(sum_logistic)
|
||||
print(derivative_fn(x_small))
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
+299
-299
File diff suppressed because it is too large
Load Diff
+371
-372
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user