1544 lines
46 KiB
Plaintext
1544 lines
46 KiB
Plaintext
TITLE: Week 40: Gradient descent methods (continued) and start Neural networks
|
|
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo, Norway
|
|
DATE: September 29-October 3, 2025
|
|
|
|
|
|
|
|
!split
|
|
===== Lecture Monday September 30, 2024 =====
|
|
!bblock
|
|
o Stochastic Gradient descent with examples and automatic differentiation
|
|
o If we get time, we start with the basics of Neural Networks, setting up the basic steps, from the simple perceptron model to the multi-layer perceptron model
|
|
o "Video of lecture":"https://youtu.be/jdJoOrCIdII"
|
|
o Whiteboard notes at URL:"https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2024/NotesSeptember30.pdf"
|
|
!eblock
|
|
|
|
!split
|
|
===== Suggested readings and videos =====
|
|
!bblock Readings and Videos:
|
|
o The lecture notes for week 40 (these notes)
|
|
o For a good discussion on gradient methods, we would like to recommend Goodfellow et al section 4.3-4.5 and sections 8.3-8.6. We will come back to the latter chapter in our discussion of Neural networks as well.
|
|
o For neural networks we recommend Goodfellow et al chapter 6 and Raschka et al chapter 2 (contains also material about gradient descent) and chapter 11 (we will use this next week)
|
|
o Video on gradient descent at URL:"https://www.youtube.com/watch?v=sDv4f4s2SB8"
|
|
o Video on stochastic gradient descent at URL:"https://www.youtube.com/watch?v=vMh0zPT0tLI"
|
|
o Neural Networks demystified at URL:"https://www.youtube.com/watch?v=bxe2T-V8XRs&list=PLiaHhY2iBX9hdHaRr6b7XevZtgZRa1PoU&ab_channel=WelchLabs"
|
|
o Building Neural Networks from scratch at URL:https://www.youtube.com/watch?v=Wo5dMEP_BbI&list=PLQVvvaa0QuDcjD5BAw2DxE6OF2tius3V3&ab_channel=sentdex"
|
|
!eblock
|
|
|
|
!split
|
|
===== Lab sessions Tuesday and Wednesday =====
|
|
!bblock Material for the active learning sessions on Tuesday and Wednesday
|
|
* Work on project 1 and discussions on how to structure your report
|
|
* No weekly exercises for week 40, project work only
|
|
* Video on how to write scientific reports recorded during one of the lab sessions at URL:"https://youtu.be/tVW1ZDmZnwM"
|
|
* A general guideline can be found at URL:"https://github.com/CompPhysics/MachineLearning/blob/master/doc/Projects/EvaluationGrading/EvaluationForm.md".
|
|
!eblock
|
|
|
|
|
|
|
|
!split
|
|
===== Automatic differentiation =====
|
|
|
|
"Automatic differentiation (AD)":"https://en.wikipedia.org/wiki/Automatic_differentiation",
|
|
also called algorithmic
|
|
differentiation or computational differentiation,is a set of
|
|
techniques to numerically evaluate the derivative of a function
|
|
specified by a computer program. AD exploits the fact that every
|
|
computer program, no matter how complicated, executes a sequence of
|
|
elementary arithmetic operations (addition, subtraction,
|
|
multiplication, division, etc.) and elementary functions (exp, log,
|
|
sin, cos, etc.). By applying the chain rule repeatedly to these
|
|
operations, derivatives of arbitrary order can be computed
|
|
automatically, accurately to working precision, and using at most a
|
|
small constant factor more arithmetic operations than the original
|
|
program.
|
|
|
|
Automatic differentiation is neither:
|
|
|
|
* Symbolic differentiation, nor
|
|
* Numerical differentiation (the method of finite differences).
|
|
|
|
Symbolic differentiation can lead to inefficient code and faces the
|
|
difficulty of converting a computer program into a single expression,
|
|
while numerical differentiation can introduce round-off errors in the
|
|
discretization process and cancellation
|
|
|
|
|
|
|
|
Python has tools for so-called _automatic differentiation_.
|
|
Consider the following example
|
|
!bt
|
|
\[
|
|
f(x) = \sin\left(2\pi x + x^2\right)
|
|
\]
|
|
!et
|
|
which has the following derivative
|
|
!bt
|
|
\[
|
|
f'(x) = \cos\left(2\pi x + x^2\right)\left(2\pi + 2x\right)
|
|
\]
|
|
!et
|
|
Using _autograd_ we have
|
|
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
|
|
# To do elementwise differentiation:
|
|
from autograd import elementwise_grad as egrad
|
|
|
|
# To plot:
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
def f(x):
|
|
return np.sin(2*np.pi*x + x**2)
|
|
|
|
def f_grad_analytic(x):
|
|
return np.cos(2*np.pi*x + x**2)*(2*np.pi + 2*x)
|
|
|
|
# Do the comparison:
|
|
x = np.linspace(0,1,1000)
|
|
|
|
f_grad = egrad(f)
|
|
|
|
computed = f_grad(x)
|
|
analytic = f_grad_analytic(x)
|
|
|
|
plt.title('Derivative computed from Autograd compared with the analytical derivative')
|
|
plt.plot(x,computed,label='autograd')
|
|
plt.plot(x,analytic,label='analytic')
|
|
|
|
plt.xlabel('x')
|
|
plt.ylabel('y')
|
|
plt.legend()
|
|
|
|
plt.show()
|
|
|
|
print("The max absolute difference is: %g"%(np.max(np.abs(computed - analytic))))
|
|
!ec
|
|
|
|
!split
|
|
===== Using autograd =====
|
|
|
|
Here we
|
|
experiment with what kind of functions Autograd is capable
|
|
of finding the gradient of. The following Python functions are just
|
|
meant to illustrate what Autograd can do, but please feel free to
|
|
experiment with other, possibly more complicated, functions as well.
|
|
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
from autograd import grad
|
|
|
|
def f1(x):
|
|
return x**3 + 1
|
|
|
|
f1_grad = grad(f1)
|
|
|
|
# Remember to send in float as argument to the computed gradient from Autograd!
|
|
a = 1.0
|
|
|
|
# See the evaluated gradient at a using autograd:
|
|
print("The gradient of f1 evaluated at a = %g using autograd is: %g"%(a,f1_grad(a)))
|
|
|
|
# Compare with the analytical derivative, that is f1'(x) = 3*x**2
|
|
grad_analytical = 3*a**2
|
|
print("The gradient of f1 evaluated at a = %g by finding the analytic expression is: %g"%(a,grad_analytical))
|
|
!ec
|
|
|
|
|
|
!split
|
|
===== Autograd with more complicated functions =====
|
|
|
|
To differentiate with respect to two (or more) arguments of a Python
|
|
function, Autograd need to know at which variable the function if
|
|
being differentiated with respect to.
|
|
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
from autograd import grad
|
|
def f2(x1,x2):
|
|
return 3*x1**3 + x2*(x1 - 5) + 1
|
|
|
|
# By sending the argument 0, Autograd will compute the derivative w.r.t the first variable, in this case x1
|
|
f2_grad_x1 = grad(f2,0)
|
|
|
|
# ... and differentiate w.r.t x2 by sending 1 as an additional arugment to grad
|
|
f2_grad_x2 = grad(f2,1)
|
|
|
|
x1 = 1.0
|
|
x2 = 3.0
|
|
|
|
print("Evaluating at x1 = %g, x2 = %g"%(x1,x2))
|
|
print("-"*30)
|
|
|
|
# Compare with the analytical derivatives:
|
|
|
|
# Derivative of f2 w.r.t x1 is: 9*x1**2 + x2:
|
|
f2_grad_x1_analytical = 9*x1**2 + x2
|
|
|
|
# Derivative of f2 w.r.t x2 is: x1 - 5:
|
|
f2_grad_x2_analytical = x1 - 5
|
|
|
|
# See the evaluated derivations:
|
|
print("The derivative of f2 w.r.t x1: %g"%( f2_grad_x1(x1,x2) ))
|
|
print("The analytical derivative of f2 w.r.t x1: %g"%( f2_grad_x1(x1,x2) ))
|
|
|
|
print()
|
|
|
|
print("The derivative of f2 w.r.t x2: %g"%( f2_grad_x2(x1,x2) ))
|
|
print("The analytical derivative of f2 w.r.t x2: %g"%( f2_grad_x2(x1,x2) ))
|
|
!ec
|
|
|
|
Note that the grad function will not produce the true gradient of the function. The true gradient of a function with two or more variables will produce a vector, where each element is the function differentiated w.r.t a variable.
|
|
|
|
|
|
!split
|
|
===== More complicated functions using the elements of their arguments directly =====
|
|
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
from autograd import grad
|
|
def f3(x): # Assumes x is an array of length 5 or higher
|
|
return 2*x[0] + 3*x[1] + 5*x[2] + 7*x[3] + 11*x[4]**2
|
|
|
|
f3_grad = grad(f3)
|
|
|
|
x = np.linspace(0,4,5)
|
|
|
|
# Print the computed gradient:
|
|
print("The computed gradient of f3 is: ", f3_grad(x))
|
|
|
|
# The analytical gradient is: (2, 3, 5, 7, 22*x[4])
|
|
f3_grad_analytical = np.array([2, 3, 5, 7, 22*x[4]])
|
|
|
|
# Print the analytical gradient:
|
|
print("The analytical gradient of f3 is: ", f3_grad_analytical)
|
|
!ec
|
|
|
|
Note that in this case, when sending an array as input argument, the
|
|
output from Autograd is another array. This is the true gradient of
|
|
the function, as opposed to the function in the previous example. By
|
|
using arrays to represent the variables, the output from Autograd
|
|
might be easier to work with, as the output is closer to what one
|
|
could expect form a gradient-evaluting function.
|
|
|
|
!split
|
|
===== Functions using mathematical functions from Numpy =====
|
|
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
from autograd import grad
|
|
def f4(x):
|
|
return np.sqrt(1+x**2) + np.exp(x) + np.sin(2*np.pi*x)
|
|
|
|
f4_grad = grad(f4)
|
|
|
|
x = 2.7
|
|
|
|
# Print the computed derivative:
|
|
print("The computed derivative of f4 at x = %g is: %g"%(x,f4_grad(x)))
|
|
|
|
# The analytical derivative is: x/sqrt(1 + x**2) + exp(x) + cos(2*pi*x)*2*pi
|
|
f4_grad_analytical = x/np.sqrt(1 + x**2) + np.exp(x) + np.cos(2*np.pi*x)*2*np.pi
|
|
|
|
# Print the analytical gradient:
|
|
print("The analytical gradient of f4 at x = %g is: %g"%(x,f4_grad_analytical))
|
|
!ec
|
|
|
|
|
|
!split
|
|
===== More autograd =====
|
|
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
from autograd import grad
|
|
def f5(x):
|
|
if x >= 0:
|
|
return x**2
|
|
else:
|
|
return -3*x + 1
|
|
|
|
f5_grad = grad(f5)
|
|
|
|
x = 2.7
|
|
|
|
# Print the computed derivative:
|
|
print("The computed derivative of f5 at x = %g is: %g"%(x,f5_grad(x)))
|
|
!ec
|
|
|
|
|
|
!split
|
|
===== And with loops =====
|
|
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
from autograd import grad
|
|
def f6_for(x):
|
|
val = 0
|
|
for i in range(10):
|
|
val = val + x**i
|
|
return val
|
|
|
|
def f6_while(x):
|
|
val = 0
|
|
i = 0
|
|
while i < 10:
|
|
val = val + x**i
|
|
i = i + 1
|
|
return val
|
|
|
|
f6_for_grad = grad(f6_for)
|
|
f6_while_grad = grad(f6_while)
|
|
|
|
x = 0.5
|
|
|
|
# Print the computed derivaties of f6_for and f6_while
|
|
print("The computed derivative of f6_for at x = %g is: %g"%(x,f6_for_grad(x)))
|
|
print("The computed derivative of f6_while at x = %g is: %g"%(x,f6_while_grad(x)))
|
|
!ec
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
from autograd import grad
|
|
# Both of the functions are implementation of the sum: sum(x**i) for i = 0, ..., 9
|
|
# The analytical derivative is: sum(i*x**(i-1))
|
|
f6_grad_analytical = 0
|
|
for i in range(10):
|
|
f6_grad_analytical += i*x**(i-1)
|
|
|
|
print("The analytical derivative of f6 at x = %g is: %g"%(x,f6_grad_analytical))
|
|
!ec
|
|
|
|
!split
|
|
===== Using recursion =====
|
|
!bc pycod
|
|
import autograd.numpy as np
|
|
from autograd import grad
|
|
|
|
def f7(n): # Assume that n is an integer
|
|
if n == 1 or n == 0:
|
|
return 1
|
|
else:
|
|
return n*f7(n-1)
|
|
|
|
f7_grad = grad(f7)
|
|
|
|
n = 2.0
|
|
|
|
print("The computed derivative of f7 at n = %d is: %g"%(n,f7_grad(n)))
|
|
|
|
# The function f7 is an implementation of the factorial of n.
|
|
# By using the product rule, one can find that the derivative is:
|
|
|
|
f7_grad_analytical = 0
|
|
for i in range(int(n)-1):
|
|
tmp = 1
|
|
for k in range(int(n)-1):
|
|
if k != i:
|
|
tmp *= (n - k)
|
|
f7_grad_analytical += tmp
|
|
|
|
print("The analytical derivative of f7 at n = %d is: %g"%(n,f7_grad_analytical))
|
|
|
|
!ec
|
|
Note that if n is equal to zero or one, Autograd will give an error message. This message appears when the output is independent on input.
|
|
|
|
|
|
!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
|
|
===== 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
|
|
|
|
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
|
|
|
|
!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 = 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
|
|
|
|
|
|
!split
|
|
===== Same code but now with momentum gradient descent =====
|
|
!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)
|
|
|
|
|
|
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
|
|
|
|
|
|
!split
|
|
===== 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 = 1000
|
|
x = np.random.rand(n,1)
|
|
y = 2.0+3*x +4*x*x
|
|
|
|
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):
|
|
Giter = 0.0
|
|
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)
|
|
Giter += gradients*gradients
|
|
update = gradients*eta/(delta+np.sqrt(Giter))
|
|
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.
|
|
|
|
!split
|
|
===== RMSprop for adaptive learning rate with Stochastic Gradient Descent =====
|
|
!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 = 1000
|
|
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 = 0.0
|
|
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)
|
|
# Accumulated gradient
|
|
# Scaling with rho the new and the previous results
|
|
Giter = (rho*Giter+(1-rho)*gradients*gradients)
|
|
# Taking the diagonal only and inverting
|
|
update = gradients*eta/(delta+np.sqrt(Giter))
|
|
# Hadamard product
|
|
theta -= update
|
|
print("theta from own RMSprop")
|
|
print(theta)
|
|
!ec
|
|
|
|
!split
|
|
===== And finally "ADAM":"https://arxiv.org/pdf/1412.6980.pdf" =====
|
|
|
|
!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 = 1000
|
|
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 parameters beta1 and beta2, see https://arxiv.org/abs/1412.6980
|
|
beta1 = 0.9
|
|
beta2 = 0.999
|
|
# Including AdaGrad parameter to avoid possible division by zero
|
|
delta = 1e-7
|
|
iter = 0
|
|
for epoch in range(n_epochs):
|
|
first_moment = 0.0
|
|
second_moment = 0.0
|
|
iter += 1
|
|
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)
|
|
# Computing moments first
|
|
first_moment = beta1*first_moment + (1-beta1)*gradients
|
|
second_moment = beta2*second_moment+(1-beta2)*gradients*gradients
|
|
first_term = first_moment/(1.0-beta1**iter)
|
|
second_term = second_moment/(1.0-beta2**iter)
|
|
# Scaling with rho the new and the previous results
|
|
update = eta*first_term/(np.sqrt(second_term)+delta)
|
|
theta -= update
|
|
print("theta from own ADAM")
|
|
print(theta)
|
|
!ec
|
|
|
|
!split
|
|
===== 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
|
|
|
|
|
|
|
|
|
|
===== 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.
|
|
|
|
=== Getting started with Jax, note the way we import numpy ===
|
|
!bc pycod
|
|
import jax
|
|
import jax.numpy as jnp
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
|
|
from jax import grad as jax_grad
|
|
!ec
|
|
|
|
|
|
=== A warm-up example ===
|
|
|
|
!bc pycod
|
|
def function(x):
|
|
return x**2
|
|
|
|
def analytical_gradient(x):
|
|
return 2*x
|
|
|
|
def gradient_descent(starting_point, learning_rate, num_iterations, solver="analytical"):
|
|
x = starting_point
|
|
trajectory_x = [x]
|
|
trajectory_y = [function(x)]
|
|
|
|
if solver == "analytical":
|
|
grad = analytical_gradient
|
|
elif solver == "jax":
|
|
grad = jax_grad(function)
|
|
x = jnp.float64(x)
|
|
learning_rate = jnp.float64(learning_rate)
|
|
|
|
for _ in range(num_iterations):
|
|
|
|
x = x - learning_rate * grad(x)
|
|
trajectory_x.append(x)
|
|
trajectory_y.append(function(x))
|
|
|
|
return trajectory_x, trajectory_y
|
|
|
|
x = np.linspace(-5, 5, 100)
|
|
plt.plot(x, function(x), label="f(x)")
|
|
|
|
descent_x, descent_y = gradient_descent(5, 0.1, 10, solver="analytical")
|
|
jax_descend_x, jax_descend_y = gradient_descent(5, 0.1, 10, solver="jax")
|
|
|
|
plt.plot(descent_x, descent_y, label="Gradient descent", marker="o")
|
|
plt.plot(jax_descend_x, jax_descend_y, label="JAX", marker="x")
|
|
!ec
|
|
|
|
=== A more advanced example ===
|
|
|
|
!bc pycod
|
|
backend = np
|
|
|
|
def function(x):
|
|
return x*backend.sin(x**2 + 1)
|
|
|
|
def analytical_gradient(x):
|
|
return backend.sin(x**2 + 1) + 2*x**2*backend.cos(x**2 + 1)
|
|
|
|
|
|
x = np.linspace(-5, 5, 100)
|
|
plt.plot(x, function(x), label="f(x)")
|
|
|
|
descent_x, descent_y = gradient_descent(1, 0.01, 300, solver="analytical")
|
|
|
|
# Change the backend to JAX
|
|
backend = jnp
|
|
jax_descend_x, jax_descend_y = gradient_descent(1, 0.01, 300, solver="jax")
|
|
|
|
plt.scatter(descent_x, descent_y, label="Gradient descent", marker="v", s=10, color="red")
|
|
plt.scatter(jax_descend_x, jax_descend_y, label="JAX", marker="x", s=5, color="black")
|
|
!ec
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
!split
|
|
===== Introduction to Neural networks =====
|
|
|
|
Artificial neural networks are computational systems that can learn to
|
|
perform tasks by considering examples, generally without being
|
|
programmed with any task-specific rules. It is supposed to mimic a
|
|
biological system, wherein neurons interact by sending signals in the
|
|
form of mathematical functions between layers. All layers can contain
|
|
an arbitrary number of neurons, and each connection is represented by
|
|
a weight variable.
|
|
|
|
|
|
!split
|
|
===== Artificial neurons =====
|
|
|
|
The field of artificial neural networks has a long history of
|
|
development, and is closely connected with the advancement of computer
|
|
science and computers in general. A model of artificial neurons was
|
|
first developed by McCulloch and Pitts in 1943 to study signal
|
|
processing in the brain and has later been refined by others. The
|
|
general idea is to mimic neural networks in the human brain, which is
|
|
composed of billions of neurons that communicate with each other by
|
|
sending electrical signals. Each neuron accumulates its incoming
|
|
signals, which must exceed an activation threshold to yield an
|
|
output. If the threshold is not overcome, the neuron remains inactive,
|
|
i.e. has zero output.
|
|
|
|
This behaviour has inspired a simple mathematical model for an artificial neuron.
|
|
|
|
!bt
|
|
\begin{equation}
|
|
y = f\left(\sum_{i=1}^n w_ix_i\right) = f(u)
|
|
label{artificialNeuron}
|
|
\end{equation}
|
|
!et
|
|
Here, the output $y$ of the neuron is the value of its activation function, which have as input
|
|
a weighted sum of signals $x_i, \dots ,x_n$ received by $n$ other neurons.
|
|
|
|
Conceptually, it is helpful to divide neural networks into four
|
|
categories:
|
|
o general purpose neural networks for supervised learning,
|
|
o neural networks designed specifically for image processing, the most prominent example of this class being Convolutional Neural Networks (CNNs),
|
|
o neural networks for sequential data such as Recurrent Neural Networks (RNNs), and
|
|
o neural networks for unsupervised learning such as Deep Boltzmann Machines.
|
|
|
|
|
|
In natural science, DNNs and CNNs have already found numerous
|
|
applications. In statistical physics, they have been applied to detect
|
|
phase transitions in 2D Ising and Potts models, lattice gauge
|
|
theories, and different phases of polymers, or solving the
|
|
Navier-Stokes equation in weather forecasting. Deep learning has also
|
|
found interesting applications in quantum physics. Various quantum
|
|
phase transitions can be detected and studied using DNNs and CNNs,
|
|
topological phases, and even non-equilibrium many-body
|
|
localization. Representing quantum states as DNNs quantum state
|
|
tomography are among some of the impressive achievements to reveal the
|
|
potential of DNNs to facilitate the study of quantum systems.
|
|
|
|
In quantum information theory, it has been shown that one can perform
|
|
gate decompositions with the help of neural.
|
|
|
|
The applications are not limited to the natural sciences. There is a
|
|
plethora of applications in essentially all disciplines, from the
|
|
humanities to life science and medicine.
|
|
|
|
!split
|
|
===== Neural network types =====
|
|
|
|
An artificial neural network (ANN), is a computational model that
|
|
consists of layers of connected neurons, or nodes or units. We will
|
|
refer to these interchangeably as units or nodes, and sometimes as
|
|
neurons.
|
|
|
|
It is supposed to mimic a biological nervous system by letting each
|
|
neuron interact with other neurons by sending signals in the form of
|
|
mathematical functions between layers. A wide variety of different
|
|
ANNs have been developed, but most of them consist of an input layer,
|
|
an output layer and eventual layers in-between, called *hidden
|
|
layers*. All layers can contain an arbitrary number of nodes, and each
|
|
connection between two nodes is associated with a weight variable.
|
|
|
|
Neural networks (also called neural nets) are neural-inspired
|
|
nonlinear models for supervised learning. As we will see, neural nets
|
|
can be viewed as natural, more powerful extensions of supervised
|
|
learning methods such as linear and logistic regression and soft-max
|
|
methods we discussed earlier.
|
|
|
|
|
|
!split
|
|
===== Feed-forward neural networks =====
|
|
|
|
The feed-forward neural network (FFNN) was the first and simplest type
|
|
of ANNs that were devised. In this network, the information moves in
|
|
only one direction: forward through the layers.
|
|
|
|
Nodes are represented by circles, while the arrows display the
|
|
connections between the nodes, including the direction of information
|
|
flow. Additionally, each arrow corresponds to a weight variable
|
|
(figure to come). We observe that each node in a layer is connected
|
|
to *all* nodes in the subsequent layer, making this a so-called
|
|
*fully-connected* FFNN.
|
|
|
|
|
|
|
|
!split
|
|
===== Convolutional Neural Network =====
|
|
|
|
A different variant of FFNNs are *convolutional neural networks*
|
|
(CNNs), which have a connectivity pattern inspired by the animal
|
|
visual cortex. Individual neurons in the visual cortex only respond to
|
|
stimuli from small sub-regions of the visual field, called a receptive
|
|
field. This makes the neurons well-suited to exploit the strong
|
|
spatially local correlation present in natural images. The response of
|
|
each neuron can be approximated mathematically as a convolution
|
|
operation. (figure to come)
|
|
|
|
Convolutional neural networks emulate the behaviour of neurons in the
|
|
visual cortex by enforcing a *local* connectivity pattern between
|
|
nodes of adjacent layers: Each node in a convolutional layer is
|
|
connected only to a subset of the nodes in the previous layer, in
|
|
contrast to the fully-connected FFNN. Often, CNNs consist of several
|
|
convolutional layers that learn local features of the input, with a
|
|
fully-connected layer at the end, which gathers all the local data and
|
|
produces the outputs. They have wide applications in image and video
|
|
recognition.
|
|
|
|
!split
|
|
===== Recurrent neural networks =====
|
|
|
|
So far we have only mentioned ANNs where information flows in one
|
|
direction: forward. *Recurrent neural networks* on the other hand,
|
|
have connections between nodes that form directed *cycles*. This
|
|
creates a form of internal memory which are able to capture
|
|
information on what has been calculated before; the output is
|
|
dependent on the previous computations. Recurrent NNs make use of
|
|
sequential information by performing the same task for every element
|
|
in a sequence, where each element depends on previous elements. An
|
|
example of such information is sentences, making recurrent NNs
|
|
especially well-suited for handwriting and speech recognition.
|
|
|
|
!split
|
|
===== Other types of networks =====
|
|
|
|
There are many other kinds of ANNs that have been developed. One type
|
|
that is specifically designed for interpolation in multidimensional
|
|
space is the radial basis function (RBF) network. RBFs are typically
|
|
made up of three layers: an input layer, a hidden layer with
|
|
non-linear radial symmetric activation functions and a linear output
|
|
layer (''linear'' here means that each node in the output layer has a
|
|
linear activation function). The layers are normally fully-connected
|
|
and there are no cycles, thus RBFs can be viewed as a type of
|
|
fully-connected FFNN. They are however usually treated as a separate
|
|
type of NN due the unusual activation functions.
|
|
|
|
!split
|
|
===== Multilayer perceptrons =====
|
|
|
|
One uses often so-called fully-connected feed-forward neural networks
|
|
with three or more layers (an input layer, one or more hidden layers
|
|
and an output layer) consisting of neurons that have non-linear
|
|
activation functions.
|
|
|
|
Such networks are often called *multilayer perceptrons* (MLPs).
|
|
|
|
!split
|
|
===== Why multilayer perceptrons? =====
|
|
|
|
According to the *Universal approximation theorem*, a feed-forward
|
|
neural network with just a single hidden layer containing a finite
|
|
number of neurons can approximate a continuous multidimensional
|
|
function to arbitrary accuracy, assuming the activation function for
|
|
the hidden layer is a _non-constant, bounded and
|
|
monotonically-increasing continuous function_.
|
|
|
|
Note that the requirements on the activation function only applies to
|
|
the hidden layer, the output nodes are always assumed to be linear, so
|
|
as to not restrict the range of output values.
|
|
|
|
|
|
!split
|
|
===== Illustration of a single perceptron model and a multi-perceptron model =====
|
|
|
|
FIGURE: [figures/nns.png, width=600 frac=0.8] In a) we show a single perceptron model while in b) we dispay a network with two hidden layers, an input layer and an output layer.
|
|
|
|
|
|
!split
|
|
===== Examples of XOR, OR and AND gates =====
|
|
|
|
|
|
|
|
Let us first try to fit various gates using standard linear
|
|
regression. The gates we are thinking of are the classical XOR, OR and
|
|
AND gates, well-known elements in computer science. The tables here
|
|
show how we can set up the inputs $x_1$ and $x_2$ in order to yield a
|
|
specific target $y_i$.
|
|
|
|
|
|
|
|
|
|
!bc pycod
|
|
"""
|
|
Simple code that tests XOR, OR and AND gates with linear regression
|
|
"""
|
|
|
|
import numpy as np
|
|
# Design matrix
|
|
X = np.array([ [1, 0, 0], [1, 0, 1], [1, 1, 0],[1, 1, 1]],dtype=np.float64)
|
|
print(f"The X.TX matrix:{X.T @ X}")
|
|
Xinv = np.linalg.pinv(X.T @ X)
|
|
print(f"The invers of X.TX matrix:{Xinv}")
|
|
|
|
# The XOR gate
|
|
yXOR = np.array( [ 0, 1 ,1, 0])
|
|
ThetaXOR = Xinv @ X.T @ yXOR
|
|
print(f"The values of theta for the XOR gate:{ThetaXOR}")
|
|
print(f"The linear regression prediction for the XOR gate:{X @ ThetaXOR}")
|
|
|
|
|
|
# The OR gate
|
|
yOR = np.array( [ 0, 1 ,1, 1])
|
|
ThetaOR = Xinv @ X.T @ yOR
|
|
print(f"The values of theta for the OR gate:{ThetaOR}")
|
|
print(f"The linear regression prediction for the OR gate:{X @ ThetaOR}")
|
|
|
|
|
|
# The OR gate
|
|
yAND = np.array( [ 0, 0 ,0, 1])
|
|
ThetaAND = Xinv @ X.T @ yAND
|
|
print(f"The values of theta for the AND gate:{ThetaAND}")
|
|
print(f"The linear regression prediction for the AND gate:{X @ ThetaAND}")
|
|
!ec
|
|
|
|
What is happening here?
|
|
|
|
!split
|
|
===== Does Logistic Regression do a better Job? =====
|
|
|
|
!bc pycod
|
|
"""
|
|
Simple code that tests XOR and OR gates with linear regression
|
|
and logistic regression
|
|
"""
|
|
|
|
import matplotlib.pyplot as plt
|
|
from sklearn.linear_model import LogisticRegression
|
|
import numpy as np
|
|
|
|
# Design matrix
|
|
X = np.array([ [1, 0, 0], [1, 0, 1], [1, 1, 0],[1, 1, 1]],dtype=np.float64)
|
|
print(f"The X.TX matrix:{X.T @ X}")
|
|
Xinv = np.linalg.pinv(X.T @ X)
|
|
print(f"The invers of X.TX matrix:{Xinv}")
|
|
|
|
# The XOR gate
|
|
yXOR = np.array( [ 0, 1 ,1, 0])
|
|
ThetaXOR = Xinv @ X.T @ yXOR
|
|
print(f"The values of theta for the XOR gate:{ThetaXOR}")
|
|
print(f"The linear regression prediction for the XOR gate:{X @ ThetaXOR}")
|
|
|
|
|
|
# The OR gate
|
|
yOR = np.array( [ 0, 1 ,1, 1])
|
|
ThetaOR = Xinv @ X.T @ yOR
|
|
print(f"The values of theta for the OR gate:{ThetaOR}")
|
|
print(f"The linear regression prediction for the OR gate:{X @ ThetaOR}")
|
|
|
|
|
|
# The OR gate
|
|
yAND = np.array( [ 0, 0 ,0, 1])
|
|
ThetaAND = Xinv @ X.T @ yAND
|
|
print(f"The values of theta for the AND gate:{ThetaAND}")
|
|
print(f"The linear regression prediction for the AND gate:{X @ ThetaAND}")
|
|
|
|
# Now we change to logistic regression
|
|
|
|
|
|
# Logistic Regression
|
|
logreg = LogisticRegression()
|
|
logreg.fit(X, yOR)
|
|
print("Test set accuracy with Logistic Regression for OR gate: {:.2f}".format(logreg.score(X,yOR)))
|
|
|
|
logreg.fit(X, yXOR)
|
|
print("Test set accuracy with Logistic Regression for XOR gate: {:.2f}".format(logreg.score(X,yXOR)))
|
|
|
|
|
|
logreg.fit(X, yAND)
|
|
print("Test set accuracy with Logistic Regression for AND gate: {:.2f}".format(logreg.score(X,yAND)))
|
|
!ec
|
|
|
|
Not exactly impressive, but somewhat better.
|
|
|
|
!split
|
|
===== Adding Neural Networks =====
|
|
|
|
!bc pycod
|
|
|
|
# and now neural networks with Scikit-Learn and the XOR
|
|
|
|
from sklearn.neural_network import MLPClassifier
|
|
from sklearn.datasets import make_classification
|
|
X, yXOR = make_classification(n_samples=100, random_state=1)
|
|
FFNN = MLPClassifier(random_state=1, max_iter=300).fit(X, yXOR)
|
|
FFNN.predict_proba(X)
|
|
print(f"Test set accuracy with Feed Forward Neural Network for XOR gate:{FFNN.score(X, yXOR)}")
|
|
|
|
!ec
|
|
|
|
|
|
|
|
!split
|
|
===== Mathematical model =====
|
|
|
|
The output $y$ is produced via the activation function $f$
|
|
!bt
|
|
\[
|
|
y = f\left(\sum_{i=1}^n w_ix_i + b_i\right) = f(z),
|
|
\]
|
|
!et
|
|
This function receives $x_i$ as inputs.
|
|
Here the activation $z=(\sum_{i=1}^n w_ix_i+b_i)$.
|
|
In an FFNN of such neurons, the *inputs* $x_i$ are the *outputs* of
|
|
the neurons in the preceding layer. Furthermore, an MLP is
|
|
fully-connected, which means that each neuron receives a weighted sum
|
|
of the outputs of *all* neurons in the previous layer.
|
|
|
|
!split
|
|
===== Mathematical model =====
|
|
|
|
First, for each node $i$ in the first hidden layer, we calculate a weighted sum $z_i^1$ of the input coordinates $x_j$,
|
|
|
|
!bt
|
|
\begin{equation} z_i^1 = \sum_{j=1}^{M} w_{ij}^1 x_j + b_i^1
|
|
\end{equation}
|
|
!et
|
|
|
|
Here $b_i$ is the so-called bias which is normally needed in
|
|
case of zero activation weights or inputs. How to fix the biases and
|
|
the weights will be discussed below. The value of $z_i^1$ is the
|
|
argument to the activation function $f_i$ of each node $i$, The
|
|
variable $M$ stands for all possible inputs to a given node $i$ in the
|
|
first layer. We define the output $y_i^1$ of all neurons in layer 1 as
|
|
|
|
!bt
|
|
\begin{equation}
|
|
y_i^1 = f(z_i^1) = f\left(\sum_{j=1}^M w_{ij}^1 x_j + b_i^1\right)
|
|
label{outputLayer1}
|
|
\end{equation}
|
|
!et
|
|
|
|
where we assume that all nodes in the same layer have identical
|
|
activation functions, hence the notation $f$. In general, we could assume in the more general case that different layers have different activation functions.
|
|
In this case we would identify these functions with a superscript $l$ for the $l$-th layer,
|
|
|
|
!bt
|
|
\begin{equation}
|
|
y_i^l = f^l(u_i^l) = f^l\left(\sum_{j=1}^{N_{l-1}} w_{ij}^l y_j^{l-1} + b_i^l\right)
|
|
label{generalLayer}
|
|
\end{equation}
|
|
!et
|
|
|
|
where $N_l$ is the number of nodes in layer $l$. When the output of
|
|
all the nodes in the first hidden layer are computed, the values of
|
|
the subsequent layer can be calculated and so forth until the output
|
|
is obtained.
|
|
|
|
|
|
|
|
!split
|
|
===== Mathematical model =====
|
|
|
|
The output of neuron $i$ in layer 2 is thus,
|
|
|
|
!bt
|
|
\begin{align}
|
|
y_i^2 &= f^2\left(\sum_{j=1}^N w_{ij}^2 y_j^1 + b_i^2\right) \\
|
|
&= f^2\left[\sum_{j=1}^N w_{ij}^2f^1\left(\sum_{k=1}^M w_{jk}^1 x_k + b_j^1\right) + b_i^2\right]
|
|
label{outputLayer2}
|
|
\end{align}
|
|
!et
|
|
where we have substituted $y_k^1$ with the inputs $x_k$. Finally, the ANN output reads
|
|
|
|
!bt
|
|
\begin{align}
|
|
y_i^3 &= f^3\left(\sum_{j=1}^N w_{ij}^3 y_j^2 + b_i^3\right) \\
|
|
&= f_3\left[\sum_{j} w_{ij}^3 f^2\left(\sum_{k} w_{jk}^2 f^1\left(\sum_{m} w_{km}^1 x_m + b_k^1\right) + b_j^2\right)
|
|
+ b_1^3\right]
|
|
\end{align}
|
|
!et
|
|
|
|
!split
|
|
===== Mathematical model =====
|
|
|
|
We can generalize this expression to an MLP with $l$ hidden
|
|
layers. The complete functional form is,
|
|
|
|
!bt
|
|
\begin{align}
|
|
&y^{l+1}_i = f^{l+1}\left[\!\sum_{j=1}^{N_l} w_{ij}^3 f^l\left(\sum_{k=1}^{N_{l-1}}w_{jk}^{l-1}\left(\dots f^1\left(\sum_{n=1}^{N_0} w_{mn}^1 x_n+ b_m^1\right)\dots\right)+b_k^2\right)+b_1^3\right] &&
|
|
label{completeNN}
|
|
\end{align}
|
|
!et
|
|
|
|
which illustrates a basic property of MLPs: The only independent
|
|
variables are the input values $x_n$.
|
|
|
|
!split
|
|
===== Mathematical model =====
|
|
|
|
This confirms that an MLP, despite its quite convoluted mathematical
|
|
form, is nothing more than an analytic function, specifically a
|
|
mapping of real-valued vectors $\hat{x} \in \mathbb{R}^n \rightarrow
|
|
\hat{y} \in \mathbb{R}^m$.
|
|
|
|
Furthermore, the flexibility and universality of an MLP can be
|
|
illustrated by realizing that the expression is essentially a nested
|
|
sum of scaled activation functions of the form
|
|
|
|
!bt
|
|
\begin{equation}
|
|
f(x) = c_1 f(c_2 x + c_3) + c_4
|
|
\end{equation}
|
|
!et
|
|
|
|
where the parameters $c_i$ are weights and biases. By adjusting these
|
|
parameters, the activation functions can be shifted up and down or
|
|
left and right, change slope or be rescaled which is the key to the
|
|
flexibility of a neural network.
|
|
|
|
!split
|
|
=== Matrix-vector notation ===
|
|
|
|
We can introduce a more convenient notation for the activations in an A NN.
|
|
|
|
Additionally, we can represent the biases and activations
|
|
as layer-wise column vectors $\hat{b}_l$ and $\hat{y}_l$, so that the $i$-th element of each vector
|
|
is the bias $b_i^l$ and activation $y_i^l$ of node $i$ in layer $l$ respectively.
|
|
|
|
We have that $\mathrm{W}_l$ is an $N_{l-1} \times N_l$ matrix, while $\hat{b}_l$ and $\hat{y}_l$ are $N_l \times 1$ column vectors.
|
|
With this notation, the sum becomes a matrix-vector multiplication, and we can write
|
|
the equation for the activations of hidden layer 2 (assuming three nodes for simplicity) as
|
|
!bt
|
|
\begin{equation}
|
|
\hat{y}_2 = f_2(\mathrm{W}_2 \hat{y}_{1} + \hat{b}_{2}) =
|
|
f_2\left(\left[\begin{array}{ccc}
|
|
w^2_{11} &w^2_{12} &w^2_{13} \\
|
|
w^2_{21} &w^2_{22} &w^2_{23} \\
|
|
w^2_{31} &w^2_{32} &w^2_{33} \\
|
|
\end{array} \right] \cdot
|
|
\left[\begin{array}{c}
|
|
y^1_1 \\
|
|
y^1_2 \\
|
|
y^1_3 \\
|
|
\end{array}\right] +
|
|
\left[\begin{array}{c}
|
|
b^2_1 \\
|
|
b^2_2 \\
|
|
b^2_3 \\
|
|
\end{array}\right]\right).
|
|
\end{equation}
|
|
!et
|
|
|
|
!split
|
|
=== Matrix-vector notation and activation ===
|
|
|
|
The activation of node $i$ in layer 2 is
|
|
|
|
!bt
|
|
\begin{equation}
|
|
y^2_i = f_2\Bigr(w^2_{i1}y^1_1 + w^2_{i2}y^1_2 + w^2_{i3}y^1_3 + b^2_i\Bigr) =
|
|
f_2\left(\sum_{j=1}^3 w^2_{ij} y_j^1 + b^2_i\right).
|
|
\end{equation}
|
|
!et
|
|
|
|
This is not just a convenient and compact notation, but also a useful
|
|
and intuitive way to think about MLPs: The output is calculated by a
|
|
series of matrix-vector multiplications and vector additions that are
|
|
used as input to the activation functions. For each operation
|
|
$\mathrm{W}_l \hat{y}_{l-1}$ we move forward one layer.
|
|
|
|
|
|
!split
|
|
=== Activation functions ===
|
|
|
|
|
|
A property that characterizes a neural network, other than its
|
|
connectivity, is the choice of activation function(s). As described
|
|
in, the following restrictions are imposed on an activation function
|
|
for a FFNN to fulfill the universal approximation theorem
|
|
|
|
* Non-constant
|
|
|
|
* Bounded
|
|
|
|
* Monotonically-increasing
|
|
|
|
* Continuous
|
|
|
|
!split
|
|
=== Activation functions, Logistic and Hyperbolic ones ===
|
|
|
|
The second requirement excludes all linear functions. Furthermore, in
|
|
a MLP with only linear activation functions, each layer simply
|
|
performs a linear transformation of its inputs.
|
|
|
|
Regardless of the number of layers, the output of the NN will be
|
|
nothing but a linear function of the inputs. Thus we need to introduce
|
|
some kind of non-linearity to the NN to be able to fit non-linear
|
|
functions Typical examples are the logistic *Sigmoid*
|
|
|
|
!bt
|
|
\[
|
|
f(x) = \frac{1}{1 + e^{-x}},
|
|
\]
|
|
!et
|
|
and the *hyperbolic tangent* function
|
|
!bt
|
|
\[
|
|
f(x) = \tanh(x)
|
|
\]
|
|
!et
|
|
|
|
!split
|
|
=== Relevance ===
|
|
|
|
The *sigmoid* function are more biologically plausible because the
|
|
output of inactive neurons are zero. Such activation function are
|
|
called *one-sided*. However, it has been shown that the hyperbolic
|
|
tangent performs better than the sigmoid for training MLPs. has
|
|
become the most popular for *deep neural networks*
|
|
|
|
!bc pycod
|
|
"""The sigmoid function (or the logistic curve) is a
|
|
function that takes any real number, z, and outputs a number (0,1).
|
|
It is useful in neural networks for assigning weights on a relative scale.
|
|
The value z is the weighted sum of parameters involved in the learning algorithm."""
|
|
|
|
import numpy
|
|
import matplotlib.pyplot as plt
|
|
import math as mt
|
|
|
|
z = numpy.arange(-5, 5, .1)
|
|
sigma_fn = numpy.vectorize(lambda z: 1/(1+numpy.exp(-z)))
|
|
sigma = sigma_fn(z)
|
|
|
|
fig = plt.figure()
|
|
ax = fig.add_subplot(111)
|
|
ax.plot(z, sigma)
|
|
ax.set_ylim([-0.1, 1.1])
|
|
ax.set_xlim([-5,5])
|
|
ax.grid(True)
|
|
ax.set_xlabel('z')
|
|
ax.set_title('sigmoid function')
|
|
|
|
plt.show()
|
|
|
|
"""Step Function"""
|
|
z = numpy.arange(-5, 5, .02)
|
|
step_fn = numpy.vectorize(lambda z: 1.0 if z >= 0.0 else 0.0)
|
|
step = step_fn(z)
|
|
|
|
fig = plt.figure()
|
|
ax = fig.add_subplot(111)
|
|
ax.plot(z, step)
|
|
ax.set_ylim([-0.5, 1.5])
|
|
ax.set_xlim([-5,5])
|
|
ax.grid(True)
|
|
ax.set_xlabel('z')
|
|
ax.set_title('step function')
|
|
|
|
plt.show()
|
|
|
|
"""Sine Function"""
|
|
z = numpy.arange(-2*mt.pi, 2*mt.pi, 0.1)
|
|
t = numpy.sin(z)
|
|
|
|
fig = plt.figure()
|
|
ax = fig.add_subplot(111)
|
|
ax.plot(z, t)
|
|
ax.set_ylim([-1.0, 1.0])
|
|
ax.set_xlim([-2*mt.pi,2*mt.pi])
|
|
ax.grid(True)
|
|
ax.set_xlabel('z')
|
|
ax.set_title('sine function')
|
|
|
|
plt.show()
|
|
|
|
"""Plots a graph of the squashing function used by a rectified linear
|
|
unit"""
|
|
z = numpy.arange(-2, 2, .1)
|
|
zero = numpy.zeros(len(z))
|
|
y = numpy.max([zero, z], axis=0)
|
|
|
|
fig = plt.figure()
|
|
ax = fig.add_subplot(111)
|
|
ax.plot(z, y)
|
|
ax.set_ylim([-2.0, 2.0])
|
|
ax.set_xlim([-2.0, 2.0])
|
|
ax.grid(True)
|
|
ax.set_xlabel('z')
|
|
ax.set_title('Rectified linear unit')
|
|
|
|
plt.show()
|
|
!ec
|
|
|
|
|