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

update

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2021-12-08 06:54:30 +01:00
parent 9fdff62ca4
commit d760cdfa85
10 changed files with 7133 additions and 4030 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/BookChapters/chapteroptimization.do.txt
+252 -39
View File
@@ -922,7 +922,31 @@ plt.show()
===== Stochastic Gradient Descent =====
===== Stochastic Gradient Descent (SGD) =====
In stochastic gradient descent, the extreme case is the case where we
have only one batch, that is we include the whole data set.
This process is called Stochastic Gradient
Descent (SGD) (or also sometimes on-line gradient descent). This is
relatively less common to see because in practice due to vectorized
code optimizations it can be computationally much more efficient to
evaluate the gradient for 100 examples, than the gradient for one
example 100 times. Even though SGD technically refers to using a
single example at a time to evaluate the gradient, you will hear
people use the term SGD even when referring to mini-batch gradient
descent (i.e. mentions of MGD for “Minibatch Gradient Descent”, or BGD
for “Batch gradient descent” are rare to see), where it is usually
assumed that mini-batches are used. The size of the mini-batch is a
hyperparameter but it is not very common to cross-validate or bootstrap it. It is
usually based on memory constraints (if any), or set to some value,
e.g. 32, 64 or 128. We use powers of 2 in practice because many
vectorized operation implementations work faster when their inputs are
sized in powers of 2.
In our notes with SGD we mean stochastic gradient descent with mini-batches.
Stochastic gradient descent (SGD) and variants thereof address some of
the shortcomings of the Gradient descent method discussed above.
@@ -954,7 +978,6 @@ minibatches. We denote these minibatches by $B_k$ where
$k=1,\cdots,n/M$.
As an example, suppose we have $10$ data points $(\mathbf{x}_1,\cdots, \mathbf{x}_{10})$
and we choose to have $M=5$ minibathces,
then each minibatch contains two data points. In particular we have
@@ -991,11 +1014,12 @@ minibathces (n/M) is commonly referred to as an epoch. Thus it is
typical to choose a number of epochs and for each epoch iterate over
the number of minibatches, as exemplified in the code below.
!bc pycod
import numpy as np
n = 100 #100 datapoints
M = 5 #size of each minibatch
M = 5 #size of each mini-batche
m = int(n/M) #number of minibatches
n_epochs = 10 #number of epochs
@@ -1017,7 +1041,6 @@ cheaper since we sum over the datapoints in the $k-th$ minibatch and not
all $n$ datapoints.
A natural question is when do we stop the search for a new minimum?
One possibility is to compute the full gradient after a given number
of epochs and check if the norm of the gradient is smaller than some
@@ -1030,7 +1053,6 @@ compare the values of the cost function and keep the $\beta$ that
gave the lowest value.
Another approach is to let the step length $\gamma_j$ depend on the
number of epochs in such a way that it becomes very small after a
reasonable time such that we do not move at all.
@@ -1071,37 +1093,41 @@ print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
!ec
We note that we have defined several hyperparameters. These are now the number of epochs, the number of mini-batches and the parameters $t_0$ and $t_1$.
=== Program for stochastic gradient ===
!bc pycod
# Importing various packages
# Importing various packages
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import SGDRegressor
m = 100
x = 2*np.random.rand(m,1)
y = 4+3*x+np.random.randn(m,1)
n = 100
x = 2*np.random.rand(n,1)
y = 4+3*x+np.random.randn(n,1)
X = np.c_[np.ones((m,1)), x]
X = np.c_[np.ones((n,1)), x]
XT_X = X.T @ X
theta_linreg = np.linalg.inv(X.T @ X) @ (X.T @ y)
print("Own inversion")
print(theta_linreg)
sgdreg = SGDRegressor(max_iter = 50, penalty=None, eta0=0.1)
sgdreg.fit(x,y.ravel())
print("sgdreg from scikit")
print(sgdreg.intercept_, sgdreg.coef_)
# 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 = 0.1
eta = 1.0/np.max(EigValues)
Niterations = 1000
for iter in range(Niterations):
gradients = 2.0/m*X.T @ ((X @ theta)-y)
gradients = 2.0/n*X.T @ ((X @ theta)-y)
theta -= eta*gradients
print("theta from own gd")
print(theta)
@@ -1111,8 +1137,9 @@ Xnew = np.c_[np.ones((2,1)), xnew]
ypredict = Xnew.dot(theta)
ypredict2 = Xnew.dot(theta_linreg)
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)
@@ -1120,16 +1147,20 @@ def learning_schedule(t):
theta = np.random.randn(2,1)
for epoch in range(n_epochs):
# Can you figure out a better way of setting up the contributions to each batch?
for i in range(m):
random_index = np.random.randint(m)
xi = X[random_index:random_index+1]
yi = y[random_index:random_index+1]
gradients = 2 * xi.T @ ((xi @ theta)-yi)
random_index = M*np.random.randint(m)
xi = X[random_index:random_index+M]
yi = y[random_index:random_index+M]
gradients = (2.0/M)* xi.T @ ((xi @ theta)-yi)
eta = learning_schedule(epoch*m+i)
theta = theta - eta*gradients
print("theta from own sdg")
print(theta)
plt.plot(xnew, ypredict, "r-")
plt.plot(xnew, ypredict2, "b-")
plt.plot(x, y ,'ro')
@@ -1142,6 +1173,13 @@ plt.show()
!ec
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. More material will be added later.
===== Momentum based GD =====
The stochastic gradient descent (SGD) is almost always used with a
@@ -1175,7 +1213,6 @@ earlier. An equivalent way of writing the updates is
where we have defined $\Delta \boldsymbol{\theta}_{t}= \boldsymbol{\theta}_t-\boldsymbol{\theta}_{t-1}$.
Let us try to get more intuition from these equations. It is helpful
to consider a simple physical analogy with a particle of mass $m$
moving in a viscous medium with drag coefficient $\mu$ and potential
@@ -1205,7 +1242,6 @@ Rearranging this equation, we can rewrite this as
!et
Notice that this equation is identical to previous one if we identify
the position of the particle, $\mathbf{w}$, with the parameters
$\boldsymbol{\theta}$. This allows us to identify the momentum
@@ -1254,6 +1290,7 @@ One of the major advantages of NAG is that it allows for the use of a larger lea
In stochastic gradient descent, with and without momentum, we still
have to specify a schedule for tuning the learning rates $\eta_t$
as a function of time. As discussed in the context of Newton's
@@ -1272,10 +1309,9 @@ Hessians.
Recently, a number of methods have been introduced that accomplish
this by tracking not only the gradient, but also the second moment of
the gradient. These methods include AdaGrad, AdaDelta, RMS-Prop, and
the gradient. These methods include AdaGrad, AdaDelta, Root Mean Squared Propagation (RMS-Prop), and
ADAM.
=== RMS prop ===
In RMS prop, in addition to keeping a running average of the first
@@ -1301,6 +1337,8 @@ directions where the norm of the gradient is consistently large. This
greatly speeds up the convergence by allowing us to use a larger
learning rate for flat directions.
=== ADAM optimizer ===
A related algorithm is the ADAM optimizer. In ADAM, we keep a running
@@ -1359,6 +1397,7 @@ update rule for this parameter is given by
* _Adaptive optimization methods don't always have good generalization._ Recent studies have shown that adaptive methods such as ADAM, RMSPorp, and AdaGrad tend to have poor generalization compared to SGD or SGD with momentum, particularly in the high-dimensional limit (i.e. the number of parameters exceeds the number of data points). Although it is not clear at this stage why these methods perform so well in training deep neural networks, simpler procedures like properly-tuned SGD may work as well or better in these applications.
===== Automatic differentiation =====
"Automatic differentiation (AD)":"https://en.wikipedia.org/wiki/Automatic_differentiation",
@@ -1540,7 +1579,6 @@ might be easier to work with, as the output is closer to what one
could expect form a gradient-evaluting function.
!bc pycod
import autograd.numpy as np
from autograd import grad
@@ -1581,7 +1619,6 @@ print("The computed derivative of f5 at x = %g is: %g"%(x,f5_grad(x)))
!ec
!bc pycod
import autograd.numpy as np
from autograd import grad
@@ -1654,21 +1691,21 @@ print("The analytical derivative of f7 at n = %d is: %g"%(n,f7_grad_analytical))
Note that if n is equal to zero or one, Autograd will give an error message. This message appears when the output is independent on input.
Autograd supports many features. However, there are some functions that are not supported (yet) by Autograd.
Autograd supports many features. However, there are some functions that is not supported (yet) by Autograd.
Assigning a value to the variable being differentiated with respect to is an example thereof.
Assigning a value to the variable being differentiated with respect to
!bc pycod
#import autograd.numpy as np
#from autograd import grad
#def f8(x): # Assume x is an array
# x[2] = 3
# return x*2
import autograd.numpy as np
from autograd import grad
def f8(x): # Assume x is an array
x[2] = 3
return x*2
#f8_grad = grad(f8)
f8_grad = grad(f8)
#x = 8.4
x = 8.4
#print("The derivative of f8 is:",f8_grad(x))
print("The derivative of f8 is:",f8_grad(x))
!ec
Here, Autograd tells us that an 'ArrayBox' does not support item assignment. The item assignment is done when the program tries to assign x[2] to the value 3. However, Autograd has implemented the computation of the derivative such that this assignment is not possible.
@@ -1718,10 +1755,186 @@ a /=b
!ec
More examples will be added, in particular how to compare autograd with own codes for the gradients.
===== Using Autograd with OLS =====
We conclude the part on optmization by showing how we can make codes
for linear regression and logistic regression using _autograd_. The
first example shows results with ordinary leats squares.
!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
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 = 1000
# define the gradient
training_gradient = grad(CostOLS)
for iter in range(Niterations):
gradients = training_gradient(theta)
theta -= eta*gradients
print("theta from own gd")
print(theta)
xnew = np.array([[0],[2]])
Xnew = np.c_[np.ones((2,1)), xnew]
ypredict = Xnew.dot(theta)
ypredict2 = Xnew.dot(theta_linreg)
plt.plot(xnew, ypredict, "r-")
plt.plot(xnew, ypredict2, "b-")
plt.plot(x, y ,'ro')
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Random numbers ')
plt.show()
!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
# 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
# Note change from previous example
def CostOLS(y,X,theta):
return np.sum((y-X @ theta)**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 = 1000
# Note that we request the derivative wrt third argument (theta, 2 here)
training_gradient = grad(CostOLS,2)
for iter in range(Niterations):
gradients = (1.0/n)*training_gradient(y, X, theta)
theta -= eta*gradients
print("theta from own gd")
print(theta)
xnew = np.array([[0],[2]])
Xnew = np.c_[np.ones((2,1)), xnew]
ypredict = Xnew.dot(theta)
ypredict2 = Xnew.dot(theta_linreg)
plt.plot(xnew, ypredict, "r-")
plt.plot(xnew, ypredict2, "b-")
plt.plot(x, y ,'ro')
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Random numbers ')
plt.show()
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)
for epoch in range(n_epochs):
# Can you figure out a better way of setting up the contributions to each batch?
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)
theta = theta - eta*gradients
print("theta from own sdg")
print(theta)
!ec
=== And Logistic Regression ===
!bc pycod
import autograd.numpy as np
from autograd import grad
def sigmoid(x):
return 0.5 * (np.tanh(x / 2.) + 1)
def logistic_predictions(weights, inputs):
# Outputs probability of a label being true according to logistic model.
return sigmoid(np.dot(inputs, weights))
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))
# 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])
# Define a function that returns gradients of training loss using Autograd.
training_gradient_fun = grad(training_loss)
# 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
print("Trained loss:", training_loss(weights))
!ec