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
2022-10-05 07:55:07 +02:00
parent b00518e915
commit 031d352181
31 changed files with 3581 additions and 2153 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/LectureNotes/_build/jupyter_execute/chapteroptimization.py
+341 -32
View File
@@ -250,7 +250,19 @@
# $\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]$. 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.
# **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 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
# ensures convexity of a function $f$. We write $D_f$ to denote the
@@ -264,7 +276,7 @@
# 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.
#
@@ -1742,7 +1754,14 @@ a*= b
a /=b
# ## Using Autograd with OLS
# ## 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
@@ -1802,13 +1821,118 @@ plt.title(r'Random numbers ')
plt.show()
# ### 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**.
# ## Same code but now with momentum gradient descent
# In[27]:
# 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)
# 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.
# In[28]:
# 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)
# ## 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**.
# In[29]:
# Using Autograd to calculate gradients using SGD
# OLS example
from random import random, seed
@@ -1884,42 +2008,227 @@ print("theta from own sdg")
print(theta)
# ### And Logistic Regression
# Here we include momentum in the standard gradient descent approach.
# In[28]:
# In[30]:
# 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)
# ### Similar (second order function now) problem but now with AdaGrad
# In[31]:
# 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)
# Running this code we note an almost perfect agreement with the results from matrix inversion.
#
# Similarly, here is our implementation of RMSprop.
# In[32]:
# 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)
# ## 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.
# In[33]:
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))