updating codes
This commit is contained in:
Morten Hjorth-Jensen
@@ -0,0 +1,67 @@
|
||||
# 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 = 100
|
||||
|
||||
# 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)
|
||||
print(np.size(gradients))
|
||||
|
||||
n_epochs = 50
|
||||
M = 5 #size of each minibatch
|
||||
m = int(n/M) #number of minibatches
|
||||
theta = np.random.randn(2,1)
|
||||
|
||||
|
||||
# Including AdaGrad
|
||||
delta = 0.000001
|
||||
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)
|
||||
# calculate squared gradient by Hadamard multiplication
|
||||
r -= gradients*gradients
|
||||
# compute update
|
||||
update = (1.0/delta+np.sqrt(r))*gradients
|
||||
theta = eta*update
|
||||
print("theta from own AdaGrad")
|
||||
print(theta)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,38 @@
|
||||
# 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)
|
||||
@@ -1,293 +0,0 @@
|
||||
TITLE: Codes
|
||||
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University
|
||||
DATE: today
|
||||
!bc pycod
|
||||
from numpy import asarray
|
||||
from numpy import arange
|
||||
from numpy.random import rand
|
||||
from numpy.random import seed
|
||||
from matplotlib import pyplot
|
||||
|
||||
# objective function
|
||||
def objective(x):
|
||||
return x**2.0
|
||||
|
||||
# derivative of objective function
|
||||
def derivative(x):
|
||||
return x * 2.0
|
||||
|
||||
# gradient descent algorithm
|
||||
def gradient_descent(objective, derivative, bounds, n_iter, step_size):
|
||||
# track all solutions
|
||||
solutions, scores = list(), list()
|
||||
# generate an initial point
|
||||
solution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
|
||||
# run the gradient descent
|
||||
for i in range(n_iter):
|
||||
# calculate gradient
|
||||
gradient = derivative(solution)
|
||||
# take a step
|
||||
solution = solution - step_size * gradient
|
||||
# evaluate candidate point
|
||||
solution_eval = objective(solution)
|
||||
# store solution
|
||||
solutions.append(solution)
|
||||
scores.append(solution_eval)
|
||||
# report progress
|
||||
print('>%d f(%s) = %.5f' % (i, solution, solution_eval))
|
||||
return [solutions, scores]
|
||||
|
||||
# seed the pseudo random number generator
|
||||
seed(4)
|
||||
# define range for input
|
||||
bounds = asarray([[-1.0, 1.0]])
|
||||
# define the total iterations
|
||||
n_iter = 30
|
||||
# define the step size
|
||||
step_size = 0.1
|
||||
# perform the gradient descent search
|
||||
solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size)
|
||||
# sample input range uniformly at 0.1 increments
|
||||
inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1)
|
||||
# compute targets
|
||||
results = objective(inputs)
|
||||
# create a line plot of input vs result
|
||||
pyplot.plot(inputs, results)
|
||||
# plot the solutions found
|
||||
pyplot.plot(solutions, scores, '.-', color='red')
|
||||
# show the plot
|
||||
pyplot.show()
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Same code but now with momentum gradient descent =====
|
||||
|
||||
!bc pycod
|
||||
from numpy import asarray
|
||||
from numpy import arange
|
||||
from numpy.random import rand
|
||||
from numpy.random import seed
|
||||
from matplotlib import pyplot
|
||||
|
||||
# objective function
|
||||
def objective(x):
|
||||
return x**2.0
|
||||
|
||||
# derivative of objective function
|
||||
def derivative(x):
|
||||
return x * 2.0
|
||||
|
||||
# gradient descent algorithm
|
||||
def gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum):
|
||||
# track all solutions
|
||||
solutions, scores = list(), list()
|
||||
# generate an initial point
|
||||
solution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])
|
||||
# keep track of the change
|
||||
change = 0.0
|
||||
# run the gradient descent
|
||||
for i in range(n_iter):
|
||||
# calculate gradient
|
||||
gradient = derivative(solution)
|
||||
# calculate update
|
||||
new_change = step_size * gradient + momentum * change
|
||||
# take a step
|
||||
solution = solution - new_change
|
||||
# save the change
|
||||
change = new_change
|
||||
# evaluate candidate point
|
||||
solution_eval = objective(solution)
|
||||
# store solution
|
||||
solutions.append(solution)
|
||||
scores.append(solution_eval)
|
||||
# report progress
|
||||
print('>%d f(%s) = %.5f' % (i, solution, solution_eval))
|
||||
return [solutions, scores]
|
||||
|
||||
# seed the pseudo random number generator
|
||||
seed(4)
|
||||
# define range for input
|
||||
bounds = asarray([[-1.0, 1.0]])
|
||||
# define the total iterations
|
||||
n_iter = 30
|
||||
# define the step size
|
||||
step_size = 0.1
|
||||
# define momentum
|
||||
momentum = 0.3
|
||||
# perform the gradient descent search with momentum
|
||||
solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum)
|
||||
# sample input range uniformly at 0.1 increments
|
||||
inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1)
|
||||
# compute targets
|
||||
results = objective(inputs)
|
||||
# create a line plot of input vs result
|
||||
pyplot.plot(inputs, results)
|
||||
# plot the solutions found
|
||||
pyplot.plot(solutions, scores, '.-', color='red')
|
||||
# show the plot
|
||||
pyplot.show()
|
||||
!ec
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== 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
|
||||
|
||||
|
||||
!split
|
||||
===== 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 = 100
|
||||
|
||||
# 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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!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
|
||||
@@ -2228,11 +2228,53 @@ for iter in range(Niterations):
|
||||
print("theta from own gd wth momentum")
|
||||
print(theta)
|
||||
|
||||
|
||||
|
||||
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== But noen of these can compete with Newton's method =====
|
||||
|
||||
!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
|
||||
|
||||
|
||||
!split
|
||||
===== 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_.
|
||||
@@ -2389,11 +2431,72 @@ print(theta)
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Same problem but now with AdaGrad =====
|
||||
!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 = 100
|
||||
|
||||
# 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)
|
||||
print(np.size(gradients))
|
||||
|
||||
n_epochs = 50
|
||||
M = 5 #size of each minibatch
|
||||
m = int(n/M) #number of minibatches
|
||||
theta = np.random.randn(2,1)
|
||||
|
||||
|
||||
# Including AdaGrad
|
||||
delta = 0.000001
|
||||
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)
|
||||
# calculate squared gradient by Hadamard multiplication
|
||||
r -= gradients*gradients
|
||||
# compute update
|
||||
update = (1.0/delta+np.sqrt(r))*gradients
|
||||
theta = eta*update
|
||||
print("theta from own AdaGrad")
|
||||
print(theta)
|
||||
|
||||
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
@@ -2459,3 +2562,10 @@ derivative_fn = grad(sum_logistic)
|
||||
print(derivative_fn(x_small))
|
||||
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== Weekend challenge =====
|
||||
|
||||
* Try to run the above codes and implement the stochastic gradient descent with RMSprop and ADAM.
|
||||
* Add a more complicated function and study the rate of convergence for the derivatives as function of the different methods
|
||||
* Extend from linear regression to logistic regression.
|
||||
|
||||
Reference in New Issue
Block a user