From 4550652273ce58c4c6d0755176de5e1cfefbc111 Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Tue, 27 Sep 2022 11:25:14 +0200 Subject: [PATCH] update week 39 --- doc/pub/week39/html/week39-bs.html | 193 ++- doc/pub/week39/html/week39-reveal.html | 472 +------- doc/pub/week39/html/week39-solarized.html | 482 +------- doc/pub/week39/html/week39.html | 482 +------- doc/pub/week39/ipynb/ipynb-week39-src.tar.gz | Bin 192 -> 192 bytes doc/pub/week39/ipynb/week39.ipynb | 1106 ++++++------------ doc/src/week39/week39.do.txt | 358 +----- 7 files changed, 635 insertions(+), 2458 deletions(-) diff --git a/doc/pub/week39/html/week39-bs.html b/doc/pub/week39/html/week39-bs.html index 83ee34d7e..fc49a1729 100644 --- a/doc/pub/week39/html/week39-bs.html +++ b/doc/pub/week39/html/week39-bs.html @@ -37,7 +37,6 @@ doconce format html week39.do.txt --html_style=bootstrap --pygments_html_style=d @@ -282,91 +265,83 @@ MathJax.Hub.Config({ Contents @@ -421,7 +396,7 @@ MathJax.Hub.Config({
  • 9
  • 10
  • ...
    • -
  • 87
    • +
  • 79
  • »
    • diff --git a/doc/pub/week39/html/week39-reveal.html b/doc/pub/week39/html/week39-reveal.html index 387d4b785..c2c9baf9a 100644 --- a/doc/pub/week39/html/week39-reveal.html +++ b/doc/pub/week39/html/week39-reveal.html @@ -198,7 +198,7 @@ MathJax.Hub.Config({

    Plan for week 39

    @@ -213,12 +213,6 @@ For a good discussion on gradient methods, we would like to recommend Goodfellow

    These sections summarize neatly what we have done till now and point to what is coming with respect to deep learning.

    -
    -

    Thursday September 29

    - -Overview Video, why do we care about gradient methods? -
    -

    Optimization, the central part of any Machine Learning algortithm

    @@ -1718,46 +1712,6 @@ minibatches. We denote these minibatches by \( B_k \) where

    -
    -

    Stochastic Gradient Descent

    - -

    Stochastic gradient descent (SGD) and variants thereof address some of -the shortcomings of the Gradient descent method discussed above. -

    - -

    The underlying idea of SGD comes from the observation that the cost -function, which we want to minimize, can almost always be written as a -sum over \( n \) data points \( \{\mathbf{x}_i\}_{i=1}^n \), -

    -

     
    -$$ -C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i, -\mathbf{\beta}). -$$ -

     
    -

    - -
    -

    Computation of gradients

    - -

    This in turn means that the gradient can be -computed as a sum over \( i \)-gradients -

    -

     
    -$$ -\nabla_\beta C(\mathbf{\beta}) = \sum_i^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}). -$$ -

     
    - -

    Stochasticity/randomness is introduced by only taking the -gradient on a subset of the data called minibatches. If there are \( n \) -data points and the size of each minibatch is \( M \), there will be \( n/M \) -minibatches. We denote these minibatches by \( B_k \) where -\( k=1,\cdots,n/M \). -

    -
    -

    SGD example

    As an example, suppose we have \( 10 \) data points \( (\mathbf{x}_1,\cdots, \mathbf{x}_{10}) \) @@ -1873,8 +1827,18 @@ 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. +reasonable time such that we do not move at all. Such approaches are +also called scaling. There are many such ways to scale the learning +rate +and discussions here. See +also +https://towardsdatascience.com/learning-rate-schedules-and-adaptive-learning-rate-methods-for-deep-learning-2c8f433990d1 +for a discussion of different scaling functions for the learning rate.

    +
    + +
    +

    Time decay rate

    As an example, let \( e = 0,1,2,3,\cdots \) denote the current epoch and let \( t_0, t_1 > 0 \) be two fixed numbers. Furthermore, let \( t = e \cdot m + i \) where \( m \) is the number of minibatches and \( i=0,\cdots,m-1 \). Then the function

     
    $$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ @@ -1934,99 +1898,6 @@ j = 0

    -
    -

    Program for stochastic gradient

    - - - -
    -
    -
    -
    -
    - 
    # 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)
-
-X = np.c_[np.ones((m,1)), 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_)
-
-
-theta = np.random.randn(2,1)
-eta = 0.1
-Niterations = 1000
-
-
-for iter in range(Niterations):
-    gradients = 2.0/m*X.T @ ((X @ theta)-y)
-    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)
-
-
-n_epochs = 50
-t0, t1 = 5, 50
-def learning_schedule(t):
-    return t0/(t+t1)
-
-theta = np.random.randn(2,1)
-
-# note: here the number of minibatches is equal to the number of points!!
-for epoch in range(n_epochs):
-    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)
-        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')
-plt.axis([0,2.0,0, 15.0])
-plt.xlabel(r'$x$')
-plt.ylabel(r'$y$')
-plt.title(r'Random numbers ')
-plt.show()
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    Challenge: try to write a similar code for a Logistic Regression case.

    -
    -

    Code with a Number of Minibatches which varies

    @@ -2078,6 +1949,7 @@ 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) @@ -2119,280 +1991,6 @@ plt.show()
    -
    -

    SGD example

    - -

    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 -\( B_1 = (\mathbf{x}_1,\mathbf{x}_2), \cdots, B_5 = -(\mathbf{x}_9,\mathbf{x}_{10}) \). Note that if you choose \( M=1 \) you -have only a single batch with all data points and on the other extreme, -you may choose \( M=n \) resulting in a minibatch for each datapoint, i.e -\( B_k = \mathbf{x}_k \). -

    - -

    The idea is now to approximate the gradient by replacing the sum over -all data points with a sum over the data points in one the minibatches -picked at random in each gradient descent step -

    -

     
    -$$ -\nabla_{\beta} -C(\mathbf{\beta}) = \sum_{i=1}^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}) \rightarrow \sum_{i \in B_k}^n \nabla_\beta -c_i(\mathbf{x}_i, \mathbf{\beta}). -$$ -

     
    -

    - -
    -

    The gradient step

    - -

    Thus a gradient descent step now looks like

    -

     
    -$$ -\beta_{j+1} = \beta_j - \gamma_j \sum_{i \in B_k}^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}) -$$ -

     
    - -

    where \( k \) is picked at random with equal -probability from \( [1,n/M] \). An iteration over the number of -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. -

    -
    - -
    -

    Simple example code

    - - - -
    -
    -
    -
    -
    - 
    import numpy as np 
-
-n = 100 #100 datapoints 
-M = 5   #size of each mini-batche
-m = int(n/M) #number of minibatches
-n_epochs = 10 #number of epochs
-
-j = 0
-for epoch in range(1,n_epochs+1):
-    for i in range(m):
-        k = np.random.randint(m) #Pick the k-th minibatch at random
-        #Compute the gradient using the data in minibatch Bk
-        #Compute new suggestion for 
-        j += 1
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    Taking the gradient only on a subset of the data has two important -benefits. First, it introduces randomness which decreases the chance -that our opmization scheme gets stuck in a local minima. Second, if -the size of the minibatches are small relative to the number of -datapoints (\( M < n \)), the computation of the gradient is much -cheaper since we sum over the datapoints in the \( k-th \) minibatch and not -all \( n \) datapoints. -

    -
    - -
    -

    When do we stop?

    - -

    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 -threshold and stop if true. However, the condition that the gradient -is zero is valid also for local minima, so this would only tell us -that we are close to a local/global minimum. However, we could also -evaluate the cost function at this point, store the result and -continue the search. If the test kicks in at a later stage we can -compare the values of the cost function and keep the \( \beta \) that -gave the lowest value. -

    -
    - -
    -

    Slightly different approach

    - -

    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. -

    - -

    As an example, let \( e = 0,1,2,3,\cdots \) denote the current epoch and let \( t_0, t_1 > 0 \) be two fixed numbers. Furthermore, let \( t = e \cdot m + i \) where \( m \) is the number of minibatches and \( i=0,\cdots,m-1 \). Then the function

     
    -$$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ -

     
    goes to zero as the number of epochs gets large. I.e. we start with a step length \( \gamma_j (0; t_0, t_1) = t_0/t_1 \) which decays in time \( t \).

    - -

    In this way we can fix the number of epochs, compute \( \beta \) and -evaluate the cost function at the end. Repeating the computation will -give a different result since the scheme is random by design. Then we -pick the final \( \beta \) that gives the lowest value of the cost -function. -

    - - - -
    -
    -
    -
    -
    - 
    import numpy as np 
-
-def step_length(t,t0,t1):
-    return t0/(t+t1)
-
-n = 100 #100 datapoints 
-M = 5   #size of each minibatch
-m = int(n/M) #number of minibatches
-n_epochs = 500 #number of epochs
-t0 = 1.0
-t1 = 10
-
-gamma_j = t0/t1
-j = 0
-for epoch in range(1,n_epochs+1):
-    for i in range(m):
-        k = np.random.randint(m) #Pick the k-th minibatch at random
-        #Compute the gradient using the data in minibatch Bk
-        #Compute new suggestion for beta
-        t = epoch*m+i
-        gamma_j = step_length(t,t0,t1)
-        j += 1
-
-print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    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

    - - - -
    -
    -
    -
    -
    - 
    # 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
-
-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.inv(X.T @ 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
-
-
-for iter in range(Niterations):
-    gradients = 2.0/n*X.T @ ((X @ theta)-y)
-    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)
-
-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 = (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')
-plt.axis([0,2.0,0, 15.0])
-plt.xlabel(r'$x$')
-plt.ylabel(r'$y$')
-plt.title(r'Random numbers ')
-plt.show()
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -

    Replace or not

    @@ -3536,6 +3134,50 @@ weights = np.array([0.0, Introducing JAX + +

    Presently, instead of using autograd, we recommend using JAX

    + +

    JAX is Autograd and XLA (Accelerated Linear Algebra)), +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.

    + + + +
    +
    +
    +
    +
    + 
    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))
+
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + diff --git a/doc/pub/week39/html/week39-solarized.html b/doc/pub/week39/html/week39-solarized.html index d6a46d51d..6e235d7fd 100644 --- a/doc/pub/week39/html/week39-solarized.html +++ b/doc/pub/week39/html/week39-solarized.html @@ -64,7 +64,6 @@ div.toc p,a { @@ -319,7 +302,7 @@ MathJax.Hub.Config({

    Plan for week 39

    @@ -332,11 +315,6 @@ For a good discussion on gradient methods, we would like to recommend Goodfellow For more discussions of project 1, chapter 5 of Goodfellow et al is a good read, in particular sections 5.1-5.5 and 5.7-5.11.

    These sections summarize neatly what we have done till now and point to what is coming with respect to deep learning.

    -









    -

    Thursday September 29

    - -Overview Video, why do we care about gradient methods? -









    Optimization, the central part of any Machine Learning algortithm

    @@ -1663,41 +1641,6 @@ C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i, $$ -









    -

    Computation of gradients

    - -

    This in turn means that the gradient can be -computed as a sum over \( i \)-gradients -

    -$$ -\nabla_\beta C(\mathbf{\beta}) = \sum_i^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}). -$$ - -

    Stochasticity/randomness is introduced by only taking the -gradient on a subset of the data called minibatches. If there are \( n \) -data points and the size of each minibatch is \( M \), there will be \( n/M \) -minibatches. We denote these minibatches by \( B_k \) where -\( k=1,\cdots,n/M \). -

    - -









    -

    Stochastic Gradient Descent

    - -

    Stochastic gradient descent (SGD) and variants thereof address some of -the shortcomings of the Gradient descent method discussed above. -

    - -

    The underlying idea of SGD comes from the observation that the cost -function, which we want to minimize, can almost always be written as a -sum over \( n \) data points \( \{\mathbf{x}_i\}_{i=1}^n \), -

    -$$ -C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i, -\mathbf{\beta}). -$$ - -









    Computation of gradients

    @@ -1824,9 +1767,18 @@ 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. +reasonable time such that we do not move at all. Such approaches are +also called scaling. There are many such ways to scale the learning +rate +and discussions here. See +also +https://towardsdatascience.com/learning-rate-schedules-and-adaptive-learning-rate-methods-for-deep-learning-2c8f433990d1 +for a discussion of different scaling functions for the learning rate.

    +









    +

    Time decay rate

    +

    As an example, let \( e = 0,1,2,3,\cdots \) denote the current epoch and let \( t_0, t_1 > 0 \) be two fixed numbers. Furthermore, let \( t = e \cdot m + i \) where \( m \) is the number of minibatches and \( i=0,\cdots,m-1 \). Then the function $$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length \( \gamma_j (0; t_0, t_1) = t_0/t_1 \) which decays in time \( t \).

    In this way we can fix the number of epochs, compute \( \beta \) and @@ -1883,98 +1835,6 @@ j = 0 -









    -

    Program for stochastic gradient

    - - - -
    -
    -
    -
    -
    - 
    # 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)
-
-X = np.c_[np.ones((m,1)), 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_)
-
-
-theta = np.random.randn(2,1)
-eta = 0.1
-Niterations = 1000
-
-
-for iter in range(Niterations):
-    gradients = 2.0/m*X.T @ ((X @ theta)-y)
-    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)
-
-
-n_epochs = 50
-t0, t1 = 5, 50
-def learning_schedule(t):
-    return t0/(t+t1)
-
-theta = np.random.randn(2,1)
-
-# note: here the number of minibatches is equal to the number of points!!
-for epoch in range(n_epochs):
-    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)
-        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')
-plt.axis([0,2.0,0, 15.0])
-plt.xlabel(r'$x$')
-plt.ylabel(r'$y$')
-plt.title(r'Random numbers ')
-plt.show()
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    Challenge: try to write a similar code for a Logistic Regression case.

    -









    Code with a Number of Minibatches which varies

    @@ -2026,6 +1886,7 @@ 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) @@ -2067,270 +1928,6 @@ plt.show() -









    -

    SGD example

    - -

    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 -\( B_1 = (\mathbf{x}_1,\mathbf{x}_2), \cdots, B_5 = -(\mathbf{x}_9,\mathbf{x}_{10}) \). Note that if you choose \( M=1 \) you -have only a single batch with all data points and on the other extreme, -you may choose \( M=n \) resulting in a minibatch for each datapoint, i.e -\( B_k = \mathbf{x}_k \). -

    - -

    The idea is now to approximate the gradient by replacing the sum over -all data points with a sum over the data points in one the minibatches -picked at random in each gradient descent step -

    -$$ -\nabla_{\beta} -C(\mathbf{\beta}) = \sum_{i=1}^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}) \rightarrow \sum_{i \in B_k}^n \nabla_\beta -c_i(\mathbf{x}_i, \mathbf{\beta}). -$$ - - -









    -

    The gradient step

    - -

    Thus a gradient descent step now looks like

    -$$ -\beta_{j+1} = \beta_j - \gamma_j \sum_{i \in B_k}^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}) -$$ - -

    where \( k \) is picked at random with equal -probability from \( [1,n/M] \). An iteration over the number of -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. -

    - -









    -

    Simple example code

    - - - -
    -
    -
    -
    -
    - 
    import numpy as np 
-
-n = 100 #100 datapoints 
-M = 5   #size of each mini-batche
-m = int(n/M) #number of minibatches
-n_epochs = 10 #number of epochs
-
-j = 0
-for epoch in range(1,n_epochs+1):
-    for i in range(m):
-        k = np.random.randint(m) #Pick the k-th minibatch at random
-        #Compute the gradient using the data in minibatch Bk
-        #Compute new suggestion for 
-        j += 1
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    Taking the gradient only on a subset of the data has two important -benefits. First, it introduces randomness which decreases the chance -that our opmization scheme gets stuck in a local minima. Second, if -the size of the minibatches are small relative to the number of -datapoints (\( M < n \)), the computation of the gradient is much -cheaper since we sum over the datapoints in the \( k-th \) minibatch and not -all \( n \) datapoints. -

    - -









    -

    When do we stop?

    - -

    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 -threshold and stop if true. However, the condition that the gradient -is zero is valid also for local minima, so this would only tell us -that we are close to a local/global minimum. However, we could also -evaluate the cost function at this point, store the result and -continue the search. If the test kicks in at a later stage we can -compare the values of the cost function and keep the \( \beta \) that -gave the lowest value. -

    - -









    -

    Slightly different approach

    - -

    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. -

    - -

    As an example, let \( e = 0,1,2,3,\cdots \) denote the current epoch and let \( t_0, t_1 > 0 \) be two fixed numbers. Furthermore, let \( t = e \cdot m + i \) where \( m \) is the number of minibatches and \( i=0,\cdots,m-1 \). Then the function $$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length \( \gamma_j (0; t_0, t_1) = t_0/t_1 \) which decays in time \( t \).

    - -

    In this way we can fix the number of epochs, compute \( \beta \) and -evaluate the cost function at the end. Repeating the computation will -give a different result since the scheme is random by design. Then we -pick the final \( \beta \) that gives the lowest value of the cost -function. -

    - - - -
    -
    -
    -
    -
    - 
    import numpy as np 
-
-def step_length(t,t0,t1):
-    return t0/(t+t1)
-
-n = 100 #100 datapoints 
-M = 5   #size of each minibatch
-m = int(n/M) #number of minibatches
-n_epochs = 500 #number of epochs
-t0 = 1.0
-t1 = 10
-
-gamma_j = t0/t1
-j = 0
-for epoch in range(1,n_epochs+1):
-    for i in range(m):
-        k = np.random.randint(m) #Pick the k-th minibatch at random
-        #Compute the gradient using the data in minibatch Bk
-        #Compute new suggestion for beta
-        t = epoch*m+i
-        gamma_j = step_length(t,t0,t1)
-        j += 1
-
-print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    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

    - - - -
    -
    -
    -
    -
    - 
    # 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
-
-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.inv(X.T @ 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
-
-
-for iter in range(Niterations):
-    gradients = 2.0/n*X.T @ ((X @ theta)-y)
-    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)
-
-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 = (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')
-plt.axis([0,2.0,0, 15.0])
-plt.xlabel(r'$x$')
-plt.ylabel(r'$y$')
-plt.title(r'Random numbers ')
-plt.show()
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -









    Replace or not

    @@ -3438,6 +3035,49 @@ weights = np.array([0.0, Introducing JAX + +

    Presently, instead of using autograd, we recommend using JAX

    + +

    JAX is Autograd and XLA (Accelerated Linear Algebra)), +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.

    + + + +
    +
    +
    +
    +
    + 
    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))
+
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    © 1999-2022, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license diff --git a/doc/pub/week39/html/week39.html b/doc/pub/week39/html/week39.html index d1efb903c..530f28186 100644 --- a/doc/pub/week39/html/week39.html +++ b/doc/pub/week39/html/week39.html @@ -141,7 +141,6 @@ div.toc p,a { @@ -396,7 +379,7 @@ MathJax.Hub.Config({

    Plan for week 39

      -
    • Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Examples on how to implement Logistic Regression
      • +
    • Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Examples on how to implement Logistic Regression and discussion of stochastic gradient descent
    • Friday: Stochastic Gradient descent with examples and automatic differentiation
    • Reading recommendations:
    @@ -409,11 +392,6 @@ For a good discussion on gradient methods, we would like to recommend Goodfellow For more discussions of project 1, chapter 5 of Goodfellow et al is a good read, in particular sections 5.1-5.5 and 5.7-5.11.

    These sections summarize neatly what we have done till now and point to what is coming with respect to deep learning.

    -









    -

    Thursday September 29

    - -Overview Video, why do we care about gradient methods? -









    Optimization, the central part of any Machine Learning algortithm

    @@ -1740,41 +1718,6 @@ C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i, $$ -









    -

    Computation of gradients

    - -

    This in turn means that the gradient can be -computed as a sum over \( i \)-gradients -

    -$$ -\nabla_\beta C(\mathbf{\beta}) = \sum_i^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}). -$$ - -

    Stochasticity/randomness is introduced by only taking the -gradient on a subset of the data called minibatches. If there are \( n \) -data points and the size of each minibatch is \( M \), there will be \( n/M \) -minibatches. We denote these minibatches by \( B_k \) where -\( k=1,\cdots,n/M \). -

    - -









    -

    Stochastic Gradient Descent

    - -

    Stochastic gradient descent (SGD) and variants thereof address some of -the shortcomings of the Gradient descent method discussed above. -

    - -

    The underlying idea of SGD comes from the observation that the cost -function, which we want to minimize, can almost always be written as a -sum over \( n \) data points \( \{\mathbf{x}_i\}_{i=1}^n \), -

    -$$ -C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i, -\mathbf{\beta}). -$$ - -









    Computation of gradients

    @@ -1901,9 +1844,18 @@ 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. +reasonable time such that we do not move at all. Such approaches are +also called scaling. There are many such ways to scale the learning +rate +and discussions here. See +also +https://towardsdatascience.com/learning-rate-schedules-and-adaptive-learning-rate-methods-for-deep-learning-2c8f433990d1 +for a discussion of different scaling functions for the learning rate.

    +









    +

    Time decay rate

    +

    As an example, let \( e = 0,1,2,3,\cdots \) denote the current epoch and let \( t_0, t_1 > 0 \) be two fixed numbers. Furthermore, let \( t = e \cdot m + i \) where \( m \) is the number of minibatches and \( i=0,\cdots,m-1 \). Then the function $$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length \( \gamma_j (0; t_0, t_1) = t_0/t_1 \) which decays in time \( t \).

    In this way we can fix the number of epochs, compute \( \beta \) and @@ -1960,98 +1912,6 @@ j = 0 -









    -

    Program for stochastic gradient

    - - - -
    -
    -
    -
    -
    - 
    # 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)
-
-X = np.c_[np.ones((m,1)), 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_)
-
-
-theta = np.random.randn(2,1)
-eta = 0.1
-Niterations = 1000
-
-
-for iter in range(Niterations):
-    gradients = 2.0/m*X.T @ ((X @ theta)-y)
-    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)
-
-
-n_epochs = 50
-t0, t1 = 5, 50
-def learning_schedule(t):
-    return t0/(t+t1)
-
-theta = np.random.randn(2,1)
-
-# note: here the number of minibatches is equal to the number of points!!
-for epoch in range(n_epochs):
-    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)
-        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')
-plt.axis([0,2.0,0, 15.0])
-plt.xlabel(r'$x$')
-plt.ylabel(r'$y$')
-plt.title(r'Random numbers ')
-plt.show()
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    Challenge: try to write a similar code for a Logistic Regression case.

    -









    Code with a Number of Minibatches which varies

    @@ -2103,6 +1963,7 @@ 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) @@ -2144,270 +2005,6 @@ plt.show() -









    -

    SGD example

    - -

    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 -\( B_1 = (\mathbf{x}_1,\mathbf{x}_2), \cdots, B_5 = -(\mathbf{x}_9,\mathbf{x}_{10}) \). Note that if you choose \( M=1 \) you -have only a single batch with all data points and on the other extreme, -you may choose \( M=n \) resulting in a minibatch for each datapoint, i.e -\( B_k = \mathbf{x}_k \). -

    - -

    The idea is now to approximate the gradient by replacing the sum over -all data points with a sum over the data points in one the minibatches -picked at random in each gradient descent step -

    -$$ -\nabla_{\beta} -C(\mathbf{\beta}) = \sum_{i=1}^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}) \rightarrow \sum_{i \in B_k}^n \nabla_\beta -c_i(\mathbf{x}_i, \mathbf{\beta}). -$$ - - -









    -

    The gradient step

    - -

    Thus a gradient descent step now looks like

    -$$ -\beta_{j+1} = \beta_j - \gamma_j \sum_{i \in B_k}^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}) -$$ - -

    where \( k \) is picked at random with equal -probability from \( [1,n/M] \). An iteration over the number of -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. -

    - -









    -

    Simple example code

    - - - -
    -
    -
    -
    -
    - 
    import numpy as np 
-
-n = 100 #100 datapoints 
-M = 5   #size of each mini-batche
-m = int(n/M) #number of minibatches
-n_epochs = 10 #number of epochs
-
-j = 0
-for epoch in range(1,n_epochs+1):
-    for i in range(m):
-        k = np.random.randint(m) #Pick the k-th minibatch at random
-        #Compute the gradient using the data in minibatch Bk
-        #Compute new suggestion for 
-        j += 1
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    Taking the gradient only on a subset of the data has two important -benefits. First, it introduces randomness which decreases the chance -that our opmization scheme gets stuck in a local minima. Second, if -the size of the minibatches are small relative to the number of -datapoints (\( M < n \)), the computation of the gradient is much -cheaper since we sum over the datapoints in the \( k-th \) minibatch and not -all \( n \) datapoints. -

    - -









    -

    When do we stop?

    - -

    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 -threshold and stop if true. However, the condition that the gradient -is zero is valid also for local minima, so this would only tell us -that we are close to a local/global minimum. However, we could also -evaluate the cost function at this point, store the result and -continue the search. If the test kicks in at a later stage we can -compare the values of the cost function and keep the \( \beta \) that -gave the lowest value. -

    - -









    -

    Slightly different approach

    - -

    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. -

    - -

    As an example, let \( e = 0,1,2,3,\cdots \) denote the current epoch and let \( t_0, t_1 > 0 \) be two fixed numbers. Furthermore, let \( t = e \cdot m + i \) where \( m \) is the number of minibatches and \( i=0,\cdots,m-1 \). Then the function $$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length \( \gamma_j (0; t_0, t_1) = t_0/t_1 \) which decays in time \( t \).

    - -

    In this way we can fix the number of epochs, compute \( \beta \) and -evaluate the cost function at the end. Repeating the computation will -give a different result since the scheme is random by design. Then we -pick the final \( \beta \) that gives the lowest value of the cost -function. -

    - - - -
    -
    -
    -
    -
    - 
    import numpy as np 
-
-def step_length(t,t0,t1):
-    return t0/(t+t1)
-
-n = 100 #100 datapoints 
-M = 5   #size of each minibatch
-m = int(n/M) #number of minibatches
-n_epochs = 500 #number of epochs
-t0 = 1.0
-t1 = 10
-
-gamma_j = t0/t1
-j = 0
-for epoch in range(1,n_epochs+1):
-    for i in range(m):
-        k = np.random.randint(m) #Pick the k-th minibatch at random
-        #Compute the gradient using the data in minibatch Bk
-        #Compute new suggestion for beta
-        t = epoch*m+i
-        gamma_j = step_length(t,t0,t1)
-        j += 1
-
-print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j))
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -

    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

    - - - -
    -
    -
    -
    -
    - 
    # 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
-
-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.inv(X.T @ 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
-
-
-for iter in range(Niterations):
-    gradients = 2.0/n*X.T @ ((X @ theta)-y)
-    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)
-
-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 = (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')
-plt.axis([0,2.0,0, 15.0])
-plt.xlabel(r'$x$')
-plt.ylabel(r'$y$')
-plt.title(r'Random numbers ')
-plt.show()
-
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    -
    - -









    Replace or not

    @@ -3515,6 +3112,49 @@ weights = np. +









    +

    Introducing JAX

    + +

    Presently, instead of using autograd, we recommend using JAX

    + +

    JAX is Autograd and XLA (Accelerated Linear Algebra)), +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.

    + + + +
    +
    +
    +
    +
    + 
    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))
+
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    © 1999-2022, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license diff --git a/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz b/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz index 403824ba80d3e57fc7b3dccb7aa8641a92b844d8..a9eadf0378b50d6f8ce7015f3a1ff2a7a3ef1c27 100644 GIT binary patch literal 192 zcmV;x06+g9iwFR2!7^h21MSbv3c@f92k@Qu6nTQtZtcXQ;0_)H5nrHVnX9^XwjH{+ zcORf9#mf+(zssMH5R!eiT5q$+-CZynLP*LO47o`8m?WO+5v3d`ODN-vaGnC82~Xkx z$b2WgwAKmJpHf#RR2J2{xqhrHKkS)afoJ}SLnSS2cAcxV0;L`1TA$&Dcutm)Y&w-g uq0tU4FnDdHK@h3~Q54cit;8j4j6NDs+bI0?GoI&p-q#*&2NNv-2mk;y=~grV literal 192 zcmV;x06+g9iwFQxrZQsy1MSaC3c@fD2H>uHia9~av*J>)3m1Zj7f5Mpqc*8YiuU&Q z0lHG$6cO@meuf!_nSHieZ?nYTeKZ?{P|6q#xk&k#h)neeV-A>6%s8Q#5I~r6f+9fk zo%GT=FKmBGU7e(MQooz)$I9}?E8Gu)&+(`BNXPUlb< u>4p|qd2N-IAan\n", + "for a discussion of different scaling functions for the learning rate." + ] + }, + { + "cell_type": "markdown", + "id": "665fc5ff", + "metadata": { + "editable": true + }, + "source": [ + "## Time decay rate\n", "\n", "As an example, let $e = 0,1,2,3,\\cdots$ denote the current epoch and let $t_0, t_1 > 0$ be two fixed numbers. Furthermore, let $t = e \\cdot m + i$ where $m$ is the number of minibatches and $i=0,\\cdots,m-1$. Then the function $$\\gamma_j(t; t_0, t_1) = \\frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length $\\gamma_j (0; t_0, t_1) = t_0/t_1$ which decays in *time* $t$.\n", "\n", @@ -2694,7 +2628,7 @@ { "cell_type": "code", "execution_count": 11, - "id": "faa43d57", + "id": "1514a558", "metadata": { "collapsed": false, "editable": true @@ -2729,104 +2663,7 @@ }, { "cell_type": "markdown", - "id": "026cdc3d", - "metadata": { - "editable": true - }, - "source": [ - "## Program for stochastic gradient" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "id": "36d689f3", - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], - "source": [ - "# Importing various packages\n", - "from math import exp, sqrt\n", - "from random import random, seed\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "from sklearn.linear_model import SGDRegressor\n", - "\n", - "m = 100\n", - "x = 2*np.random.rand(m,1)\n", - "y = 4+3*x+np.random.randn(m,1)\n", - "\n", - "X = np.c_[np.ones((m,1)), x]\n", - "theta_linreg = np.linalg.inv(X.T @ X) @ (X.T @ y)\n", - "print(\"Own inversion\")\n", - "print(theta_linreg)\n", - "sgdreg = SGDRegressor(max_iter = 50, penalty=None, eta0=0.1)\n", - "sgdreg.fit(x,y.ravel())\n", - "print(\"sgdreg from scikit\")\n", - "print(sgdreg.intercept_, sgdreg.coef_)\n", - "\n", - "\n", - "theta = np.random.randn(2,1)\n", - "eta = 0.1\n", - "Niterations = 1000\n", - "\n", - "\n", - "for iter in range(Niterations):\n", - " gradients = 2.0/m*X.T @ ((X @ theta)-y)\n", - " theta -= eta*gradients\n", - "print(\"theta from own gd\")\n", - "print(theta)\n", - "\n", - "xnew = np.array([[0],[2]])\n", - "Xnew = np.c_[np.ones((2,1)), xnew]\n", - "ypredict = Xnew.dot(theta)\n", - "ypredict2 = Xnew.dot(theta_linreg)\n", - "\n", - "\n", - "n_epochs = 50\n", - "t0, t1 = 5, 50\n", - "def learning_schedule(t):\n", - " return t0/(t+t1)\n", - "\n", - "theta = np.random.randn(2,1)\n", - "\n", - "# note: here the number of minibatches is equal to the number of points!!\n", - "for epoch in range(n_epochs):\n", - " for i in range(m):\n", - " random_index = np.random.randint(m)\n", - " xi = X[random_index:random_index+1]\n", - " yi = y[random_index:random_index+1]\n", - " gradients = 2 * xi.T @ ((xi @ theta)-yi)\n", - " eta = learning_schedule(epoch*m+i)\n", - " theta = theta - eta*gradients\n", - "print(\"theta from own sdg\")\n", - "print(theta)\n", - "\n", - "plt.plot(xnew, ypredict, \"r-\")\n", - "plt.plot(xnew, ypredict2, \"b-\")\n", - "plt.plot(x, y ,'ro')\n", - "plt.axis([0,2.0,0, 15.0])\n", - "plt.xlabel(r'$x$')\n", - "plt.ylabel(r'$y$')\n", - "plt.title(r'Random numbers ')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "e352b052", - "metadata": { - "editable": true - }, - "source": [ - "**Challenge**: try to write a similar code for a Logistic Regression case." - ] - }, - { - "cell_type": "markdown", - "id": "6df2d658", + "id": "3c37d935", "metadata": { "editable": true }, @@ -2838,8 +2675,8 @@ }, { "cell_type": "code", - "execution_count": 13, - "id": "ef94f697", + "execution_count": 12, + "id": "358337d1", "metadata": { "collapsed": false, "editable": true @@ -2886,6 +2723,7 @@ "M = 5 #size of each minibatch\n", "m = int(n/M) #number of minibatches\n", "t0, t1 = 5, 50\n", + "\n", "def learning_schedule(t):\n", " return t0/(t+t1)\n", "\n", @@ -2915,316 +2753,7 @@ }, { "cell_type": "markdown", - "id": "1a493e7e", - "metadata": { - "editable": true - }, - "source": [ - "## SGD example\n", - "\n", - "As an example, suppose we have $10$ data points $(\\mathbf{x}_1,\\cdots, \\mathbf{x}_{10})$ \n", - "and we choose to have $M=5$ minibathces,\n", - "then each minibatch contains two data points. In particular we have\n", - "$B_1 = (\\mathbf{x}_1,\\mathbf{x}_2), \\cdots, B_5 =\n", - "(\\mathbf{x}_9,\\mathbf{x}_{10})$. Note that if you choose $M=1$ you\n", - "have only a single batch with all data points and on the other extreme,\n", - "you may choose $M=n$ resulting in a minibatch for each datapoint, i.e\n", - "$B_k = \\mathbf{x}_k$.\n", - "\n", - "The idea is now to approximate the gradient by replacing the sum over\n", - "all data points with a sum over the data points in one the minibatches\n", - "picked at random in each gradient descent step" - ] - }, - { - "cell_type": "markdown", - "id": "df7909bc", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\nabla_{\\beta}\n", - "C(\\mathbf{\\beta}) = \\sum_{i=1}^n \\nabla_\\beta c_i(\\mathbf{x}_i,\n", - "\\mathbf{\\beta}) \\rightarrow \\sum_{i \\in B_k}^n \\nabla_\\beta\n", - "c_i(\\mathbf{x}_i, \\mathbf{\\beta}).\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "afdf2e10", - "metadata": { - "editable": true - }, - "source": [ - "## The gradient step\n", - "\n", - "Thus a gradient descent step now looks like" - ] - }, - { - "cell_type": "markdown", - "id": "c2e37d69", - "metadata": { - "editable": true - }, - "source": [ - "$$\n", - "\\beta_{j+1} = \\beta_j - \\gamma_j \\sum_{i \\in B_k}^n \\nabla_\\beta c_i(\\mathbf{x}_i,\n", - "\\mathbf{\\beta})\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "6a125356", - "metadata": { - "editable": true - }, - "source": [ - "where $k$ is picked at random with equal\n", - "probability from $[1,n/M]$. An iteration over the number of\n", - "minibathces (n/M) is commonly referred to as an epoch. Thus it is\n", - "typical to choose a number of epochs and for each epoch iterate over\n", - "the number of minibatches, as exemplified in the code below." - ] - }, - { - "cell_type": "markdown", - "id": "abaf8e3d", - "metadata": { - "editable": true - }, - "source": [ - "## Simple example code" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "id": "4628460a", - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], - "source": [ - "import numpy as np \n", - "\n", - "n = 100 #100 datapoints \n", - "M = 5 #size of each mini-batche\n", - "m = int(n/M) #number of minibatches\n", - "n_epochs = 10 #number of epochs\n", - "\n", - "j = 0\n", - "for epoch in range(1,n_epochs+1):\n", - " for i in range(m):\n", - " k = np.random.randint(m) #Pick the k-th minibatch at random\n", - " #Compute the gradient using the data in minibatch Bk\n", - " #Compute new suggestion for \n", - " j += 1" - ] - }, - { - "cell_type": "markdown", - "id": "ec7e3114", - "metadata": { - "editable": true - }, - "source": [ - "Taking the gradient only on a subset of the data has two important\n", - "benefits. First, it introduces randomness which decreases the chance\n", - "that our opmization scheme gets stuck in a local minima. Second, if\n", - "the size of the minibatches are small relative to the number of\n", - "datapoints ($M < n$), the computation of the gradient is much\n", - "cheaper since we sum over the datapoints in the $k-th$ minibatch and not\n", - "all $n$ datapoints." - ] - }, - { - "cell_type": "markdown", - "id": "8edd25f6", - "metadata": { - "editable": true - }, - "source": [ - "## When do we stop?\n", - "\n", - "A natural question is when do we stop the search for a new minimum?\n", - "One possibility is to compute the full gradient after a given number\n", - "of epochs and check if the norm of the gradient is smaller than some\n", - "threshold and stop if true. However, the condition that the gradient\n", - "is zero is valid also for local minima, so this would only tell us\n", - "that we are close to a local/global minimum. However, we could also\n", - "evaluate the cost function at this point, store the result and\n", - "continue the search. If the test kicks in at a later stage we can\n", - "compare the values of the cost function and keep the $\\beta$ that\n", - "gave the lowest value." - ] - }, - { - "cell_type": "markdown", - "id": "4a05e627", - "metadata": { - "editable": true - }, - "source": [ - "## Slightly different approach\n", - "\n", - "Another approach is to let the step length $\\gamma_j$ depend on the\n", - "number of epochs in such a way that it becomes very small after a\n", - "reasonable time such that we do not move at all.\n", - "\n", - "As an example, let $e = 0,1,2,3,\\cdots$ denote the current epoch and let $t_0, t_1 > 0$ be two fixed numbers. Furthermore, let $t = e \\cdot m + i$ where $m$ is the number of minibatches and $i=0,\\cdots,m-1$. Then the function $$\\gamma_j(t; t_0, t_1) = \\frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length $\\gamma_j (0; t_0, t_1) = t_0/t_1$ which decays in *time* $t$.\n", - "\n", - "In this way we can fix the number of epochs, compute $\\beta$ and\n", - "evaluate the cost function at the end. Repeating the computation will\n", - "give a different result since the scheme is random by design. Then we\n", - "pick the final $\\beta$ that gives the lowest value of the cost\n", - "function." - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "id": "c016a06b", - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], - "source": [ - "import numpy as np \n", - "\n", - "def step_length(t,t0,t1):\n", - " return t0/(t+t1)\n", - "\n", - "n = 100 #100 datapoints \n", - "M = 5 #size of each minibatch\n", - "m = int(n/M) #number of minibatches\n", - "n_epochs = 500 #number of epochs\n", - "t0 = 1.0\n", - "t1 = 10\n", - "\n", - "gamma_j = t0/t1\n", - "j = 0\n", - "for epoch in range(1,n_epochs+1):\n", - " for i in range(m):\n", - " k = np.random.randint(m) #Pick the k-th minibatch at random\n", - " #Compute the gradient using the data in minibatch Bk\n", - " #Compute new suggestion for beta\n", - " t = epoch*m+i\n", - " gamma_j = step_length(t,t0,t1)\n", - " j += 1\n", - "\n", - "print(\"gamma_j after %d epochs: %g\" % (n_epochs,gamma_j))" - ] - }, - { - "cell_type": "markdown", - "id": "547b28ad", - "metadata": { - "editable": true - }, - "source": [ - "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$." - ] - }, - { - "cell_type": "markdown", - "id": "ca68cd6e", - "metadata": { - "editable": true - }, - "source": [ - "## Program for stochastic gradient" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "id": "1890a370", - "metadata": { - "collapsed": false, - "editable": true - }, - "outputs": [], - "source": [ - "# Importing various packages\n", - "# Importing various packages\n", - "from math import exp, sqrt\n", - "from random import random, seed\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "\n", - "n = 100\n", - "x = 2*np.random.rand(n,1)\n", - "y = 4+3*x+np.random.randn(n,1)\n", - "\n", - "X = np.c_[np.ones((n,1)), x]\n", - "XT_X = X.T @ X\n", - "theta_linreg = np.linalg.inv(X.T @ X) @ (X.T @ y)\n", - "print(\"Own inversion\")\n", - "print(theta_linreg)\n", - "# Hessian matrix\n", - "H = (2.0/n)* XT_X\n", - "EigValues, EigVectors = np.linalg.eig(H)\n", - "print(f\"Eigenvalues of Hessian Matrix:{EigValues}\")\n", - "\n", - "theta = np.random.randn(2,1)\n", - "eta = 1.0/np.max(EigValues)\n", - "Niterations = 1000\n", - "\n", - "\n", - "for iter in range(Niterations):\n", - " gradients = 2.0/n*X.T @ ((X @ theta)-y)\n", - " theta -= eta*gradients\n", - "print(\"theta from own gd\")\n", - "print(theta)\n", - "\n", - "xnew = np.array([[0],[2]])\n", - "Xnew = np.c_[np.ones((2,1)), xnew]\n", - "ypredict = Xnew.dot(theta)\n", - "ypredict2 = Xnew.dot(theta_linreg)\n", - "\n", - "n_epochs = 50\n", - "M = 5 #size of each minibatch\n", - "m = int(n/M) #number of minibatches\n", - "t0, t1 = 5, 50\n", - "def learning_schedule(t):\n", - " return t0/(t+t1)\n", - "\n", - "theta = np.random.randn(2,1)\n", - "\n", - "for epoch in range(n_epochs):\n", - "# Can you figure out a better way of setting up the contributions to each batch?\n", - " for i in range(m):\n", - " random_index = M*np.random.randint(m)\n", - " xi = X[random_index:random_index+M]\n", - " yi = y[random_index:random_index+M]\n", - " gradients = (2.0/M)* xi.T @ ((xi @ theta)-yi)\n", - " eta = learning_schedule(epoch*m+i)\n", - " theta = theta - eta*gradients\n", - "print(\"theta from own sdg\")\n", - "print(theta)\n", - "\n", - "\n", - "\n", - "\n", - "plt.plot(xnew, ypredict, \"r-\")\n", - "plt.plot(xnew, ypredict2, \"b-\")\n", - "plt.plot(x, y ,'ro')\n", - "plt.axis([0,2.0,0, 15.0])\n", - "plt.xlabel(r'$x$')\n", - "plt.ylabel(r'$y$')\n", - "plt.title(r'Random numbers ')\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "49371ea1", + "id": "e9f0cd69", "metadata": { "editable": true }, @@ -3239,7 +2768,7 @@ }, { "cell_type": "markdown", - "id": "6f44130c", + "id": "c653bcf7", "metadata": { "editable": true }, @@ -3254,7 +2783,7 @@ }, { "cell_type": "markdown", - "id": "3dfe1344", + "id": "9104cfd9", "metadata": { "editable": true }, @@ -3266,7 +2795,7 @@ }, { "cell_type": "markdown", - "id": "7438ccdb", + "id": "40d4c6f1", "metadata": { "editable": true }, @@ -3284,7 +2813,7 @@ }, { "cell_type": "markdown", - "id": "f6c3290d", + "id": "a7314538", "metadata": { "editable": true }, @@ -3303,7 +2832,7 @@ }, { "cell_type": "markdown", - "id": "3e3bf9cd", + "id": "f1fe0cd2", "metadata": { "editable": true }, @@ -3315,7 +2844,7 @@ }, { "cell_type": "markdown", - "id": "30fb1856", + "id": "0474ed6d", "metadata": { "editable": true }, @@ -3325,7 +2854,7 @@ }, { "cell_type": "markdown", - "id": "59802bd4", + "id": "f5f149c7", "metadata": { "editable": true }, @@ -3341,7 +2870,7 @@ }, { "cell_type": "markdown", - "id": "4a93e7c9", + "id": "a5d75f0f", "metadata": { "editable": true }, @@ -3353,7 +2882,7 @@ }, { "cell_type": "markdown", - "id": "2bd5f401", + "id": "9cc78aa1", "metadata": { "editable": true }, @@ -3363,7 +2892,7 @@ }, { "cell_type": "markdown", - "id": "8a70082f", + "id": "4d366fb4", "metadata": { "editable": true }, @@ -3375,7 +2904,7 @@ }, { "cell_type": "markdown", - "id": "18808d48", + "id": "a8c4e049", "metadata": { "editable": true }, @@ -3385,7 +2914,7 @@ }, { "cell_type": "markdown", - "id": "26406ac8", + "id": "39623005", "metadata": { "editable": true }, @@ -3397,7 +2926,7 @@ }, { "cell_type": "markdown", - "id": "c462f769", + "id": "b24d23a3", "metadata": { "editable": true }, @@ -3413,7 +2942,7 @@ }, { "cell_type": "markdown", - "id": "6bf678e7", + "id": "a2a5babc", "metadata": { "editable": true }, @@ -3425,7 +2954,7 @@ }, { "cell_type": "markdown", - "id": "4cdbd09d", + "id": "d85f23d3", "metadata": { "editable": true }, @@ -3458,7 +2987,7 @@ }, { "cell_type": "markdown", - "id": "5d252515", + "id": "5dfc6524", "metadata": { "editable": true }, @@ -3470,7 +2999,7 @@ }, { "cell_type": "markdown", - "id": "596d58a2", + "id": "ec075e9c", "metadata": { "editable": true }, @@ -3488,7 +3017,7 @@ }, { "cell_type": "markdown", - "id": "50e01c2e", + "id": "4895b015", "metadata": { "editable": true }, @@ -3498,7 +3027,7 @@ }, { "cell_type": "markdown", - "id": "386f9167", + "id": "0f1c8078", "metadata": { "editable": true }, @@ -3529,7 +3058,7 @@ }, { "cell_type": "markdown", - "id": "709b2114", + "id": "6db13950", "metadata": { "editable": true }, @@ -3544,7 +3073,7 @@ }, { "cell_type": "markdown", - "id": "116ee962", + "id": "40c09f63", "metadata": { "editable": true }, @@ -3562,7 +3091,7 @@ }, { "cell_type": "markdown", - "id": "fe810dcf", + "id": "2e48f563", "metadata": { "editable": true }, @@ -3574,7 +3103,7 @@ }, { "cell_type": "markdown", - "id": "cd81ad67", + "id": "2cfde768", "metadata": { "editable": true }, @@ -3586,7 +3115,7 @@ }, { "cell_type": "markdown", - "id": "06498abc", + "id": "6ae6eb5a", "metadata": { "editable": true }, @@ -3604,7 +3133,7 @@ }, { "cell_type": "markdown", - "id": "8a3c0446", + "id": "3d0ded53", "metadata": { "editable": true }, @@ -3627,7 +3156,7 @@ }, { "cell_type": "markdown", - "id": "82aeaa28", + "id": "9961ce37", "metadata": { "editable": true }, @@ -3645,7 +3174,7 @@ }, { "cell_type": "markdown", - "id": "8c71ca82", + "id": "93752a86", "metadata": { "editable": true }, @@ -3657,7 +3186,7 @@ }, { "cell_type": "markdown", - "id": "f0b27db1", + "id": "7c6a2f1c", "metadata": { "editable": true }, @@ -3669,7 +3198,7 @@ }, { "cell_type": "markdown", - "id": "5e3da9e9", + "id": "1c3f2b1d", "metadata": { "editable": true }, @@ -3681,7 +3210,7 @@ }, { "cell_type": "markdown", - "id": "6c5c06b4", + "id": "94535f21", "metadata": { "editable": true }, @@ -3693,7 +3222,7 @@ }, { "cell_type": "markdown", - "id": "3b247467", + "id": "472665ef", "metadata": { "editable": true }, @@ -3705,7 +3234,7 @@ }, { "cell_type": "markdown", - "id": "22079281", + "id": "daf0557f", "metadata": { "editable": true }, @@ -3722,7 +3251,7 @@ }, { "cell_type": "markdown", - "id": "a24eece5", + "id": "e0c15ef4", "metadata": { "editable": true }, @@ -3741,7 +3270,7 @@ }, { "cell_type": "markdown", - "id": "634c6f26", + "id": "3ce44bb3", "metadata": { "editable": true }, @@ -3753,7 +3282,7 @@ }, { "cell_type": "markdown", - "id": "71ef4fbf", + "id": "c3afdb63", "metadata": { "editable": true }, @@ -3773,7 +3302,7 @@ }, { "cell_type": "markdown", - "id": "40720b38", + "id": "52c05554", "metadata": { "editable": true }, @@ -3811,7 +3340,7 @@ }, { "cell_type": "markdown", - "id": "5c07021c", + "id": "3c43ffc0", "metadata": { "editable": true }, @@ -3823,7 +3352,7 @@ }, { "cell_type": "markdown", - "id": "ae7a6c62", + "id": "224a2860", "metadata": { "editable": true }, @@ -3833,7 +3362,7 @@ }, { "cell_type": "markdown", - "id": "5b9164e2", + "id": "bd6bca05", "metadata": { "editable": true }, @@ -3845,7 +3374,7 @@ }, { "cell_type": "markdown", - "id": "d56054c0", + "id": "ca1225a2", "metadata": { "editable": true }, @@ -3855,8 +3384,8 @@ }, { "cell_type": "code", - "execution_count": 17, - "id": "15d3d863", + "execution_count": 13, + "id": "1160dbda", "metadata": { "collapsed": false, "editable": true @@ -3901,7 +3430,7 @@ }, { "cell_type": "markdown", - "id": "5ab6d83d", + "id": "61789cf2", "metadata": { "editable": true }, @@ -3917,8 +3446,8 @@ }, { "cell_type": "code", - "execution_count": 18, - "id": "751d1395", + "execution_count": 14, + "id": "ff54127a", "metadata": { "collapsed": false, "editable": true @@ -3946,7 +3475,7 @@ }, { "cell_type": "markdown", - "id": "838d2638", + "id": "9c0b5769", "metadata": { "editable": true }, @@ -3960,8 +3489,8 @@ }, { "cell_type": "code", - "execution_count": 19, - "id": "2241609f", + "execution_count": 15, + "id": "15addd82", "metadata": { "collapsed": false, "editable": true @@ -4005,7 +3534,7 @@ }, { "cell_type": "markdown", - "id": "b8a70461", + "id": "553b9c37", "metadata": { "editable": true }, @@ -4015,7 +3544,7 @@ }, { "cell_type": "markdown", - "id": "a5bc4872", + "id": "59a5e9dd", "metadata": { "editable": true }, @@ -4025,8 +3554,8 @@ }, { "cell_type": "code", - "execution_count": 20, - "id": "33f0541e", + "execution_count": 16, + "id": "7bc0712b", "metadata": { "collapsed": false, "editable": true @@ -4054,7 +3583,7 @@ }, { "cell_type": "markdown", - "id": "50e2ba5d", + "id": "2d35dfcf", "metadata": { "editable": true }, @@ -4069,7 +3598,7 @@ }, { "cell_type": "markdown", - "id": "fe1a3058", + "id": "e9196a8d", "metadata": { "editable": true }, @@ -4079,8 +3608,8 @@ }, { "cell_type": "code", - "execution_count": 21, - "id": "fda91b69", + "execution_count": 17, + "id": "1f8e358e", "metadata": { "collapsed": false, "editable": true @@ -4108,7 +3637,7 @@ }, { "cell_type": "markdown", - "id": "1a48d07c", + "id": "d7ffc201", "metadata": { "editable": true }, @@ -4118,8 +3647,8 @@ }, { "cell_type": "code", - "execution_count": 22, - "id": "6c8baf8d", + "execution_count": 18, + "id": "04eaacee", "metadata": { "collapsed": false, "editable": true @@ -4144,7 +3673,7 @@ }, { "cell_type": "markdown", - "id": "3ac7c6c0", + "id": "9103d45c", "metadata": { "editable": true }, @@ -4154,8 +3683,8 @@ }, { "cell_type": "code", - "execution_count": 23, - "id": "ac60cd6c", + "execution_count": 19, + "id": "d458260f", "metadata": { "collapsed": false, "editable": true @@ -4190,8 +3719,8 @@ }, { "cell_type": "code", - "execution_count": 24, - "id": "59735281", + "execution_count": 20, + "id": "63951ce2", "metadata": { "collapsed": false, "editable": true @@ -4211,7 +3740,7 @@ }, { "cell_type": "markdown", - "id": "b247b6ad", + "id": "cfa9377a", "metadata": { "editable": true }, @@ -4221,8 +3750,8 @@ }, { "cell_type": "code", - "execution_count": 25, - "id": "da05505a", + "execution_count": 21, + "id": "df5d6a99", "metadata": { "collapsed": false, "editable": true @@ -4260,7 +3789,7 @@ }, { "cell_type": "markdown", - "id": "feb7935f", + "id": "d2ebfe71", "metadata": { "editable": true }, @@ -4270,7 +3799,7 @@ }, { "cell_type": "markdown", - "id": "714f2b8d", + "id": "8241e706", "metadata": { "editable": true }, @@ -4283,8 +3812,8 @@ }, { "cell_type": "code", - "execution_count": 26, - "id": "0ac7f4e6", + "execution_count": 22, + "id": "2e7bcd08", "metadata": { "collapsed": false, "editable": true @@ -4306,7 +3835,7 @@ }, { "cell_type": "markdown", - "id": "2db47b66", + "id": "72674415", "metadata": { "editable": true }, @@ -4316,7 +3845,7 @@ }, { "cell_type": "markdown", - "id": "0c29b75a", + "id": "ce70a961", "metadata": { "editable": true }, @@ -4326,8 +3855,8 @@ }, { "cell_type": "code", - "execution_count": 27, - "id": "6cfc0ef8", + "execution_count": 23, + "id": "1096a064", "metadata": { "collapsed": false, "editable": true @@ -4349,7 +3878,7 @@ }, { "cell_type": "markdown", - "id": "ee6bf2f1", + "id": "6287deb6", "metadata": { "editable": true }, @@ -4361,8 +3890,8 @@ }, { "cell_type": "code", - "execution_count": 28, - "id": "37c14206", + "execution_count": 24, + "id": "ff505ddd", "metadata": { "collapsed": false, "editable": true @@ -4387,7 +3916,7 @@ }, { "cell_type": "markdown", - "id": "645d9ee0", + "id": "6d424a52", "metadata": { "editable": true }, @@ -4398,8 +3927,8 @@ }, { "cell_type": "code", - "execution_count": 29, - "id": "27ffa613", + "execution_count": 25, + "id": "e2bfa0ec", "metadata": { "collapsed": false, "editable": true @@ -4414,7 +3943,7 @@ }, { "cell_type": "markdown", - "id": "32628932", + "id": "4ca92eaf", "metadata": { "editable": true }, @@ -4428,8 +3957,8 @@ }, { "cell_type": "code", - "execution_count": 30, - "id": "47a1833d", + "execution_count": 26, + "id": "a87b6ce0", "metadata": { "collapsed": false, "editable": true @@ -4489,7 +4018,7 @@ }, { "cell_type": "markdown", - "id": "d4dde158", + "id": "022f10f3", "metadata": { "editable": true }, @@ -4500,8 +4029,8 @@ }, { "cell_type": "code", - "execution_count": 31, - "id": "b86c1477", + "execution_count": 27, + "id": "3cd687a8", "metadata": { "collapsed": false, "editable": true @@ -4585,7 +4114,7 @@ }, { "cell_type": "markdown", - "id": "27576d79", + "id": "e843956a", "metadata": { "editable": true }, @@ -4595,8 +4124,8 @@ }, { "cell_type": "code", - "execution_count": 32, - "id": "c6d0db97", + "execution_count": 28, + "id": "3e715ead", "metadata": { "collapsed": false, "editable": true @@ -4637,6 +4166,45 @@ "\n", "print(\"Trained loss:\", training_loss(weights))" ] + }, + { + "cell_type": "markdown", + "id": "117c6315", + "metadata": { + "editable": true + }, + "source": [ + "## Introducing [JAX](https://jax.readthedocs.io/en/latest/)\n", + "\n", + "Presently, instead of using **autograd**, we recommend using [JAX](https://jax.readthedocs.io/en/latest/)\n", + "\n", + "**JAX** is Autograd and [XLA (Accelerated Linear Algebra))](https://www.tensorflow.org/xla),\n", + "brought together for high-performance numerical computing and machine learning research.\n", + "It provides composable transformations of Python+NumPy programs: differentiate, vectorize, parallelize, Just-In-Time compile to GPU/TPU, and more.\n", + "\n", + "Here's a simple example on how you can use **JAX** to compute the derivate of the logistic function." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "59eedf62", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "import jax.numpy as jnp\n", + "from jax import grad, jit, vmap\n", + "\n", + "def sum_logistic(x):\n", + " return jnp.sum(1.0 / (1.0 + jnp.exp(-x)))\n", + "\n", + "x_small = jnp.arange(3.)\n", + "derivative_fn = grad(sum_logistic)\n", + "print(derivative_fn(x_small))" + ] } ], "metadata": {}, diff --git a/doc/src/week39/week39.do.txt b/doc/src/week39/week39.do.txt index c41d0a431..6a9ed5815 100644 --- a/doc/src/week39/week39.do.txt +++ b/doc/src/week39/week39.do.txt @@ -5,7 +5,7 @@ DATE: today !split ===== Plan for week 39 ===== -* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Examples on how to implement Logistic Regression +* Thursday: Repetition of Logistic regression equations and classification problems and discussion of Gradient methods. Examples on how to implement Logistic Regression and discussion of stochastic gradient descent * Friday: Stochastic Gradient descent with examples and automatic differentiation * Reading recommendations: @@ -18,11 +18,6 @@ For Stochastic Gradient Descent, we recommend chapter 4 of Geron's text. _For more discussions of project 1, chapter 5 of Goodfellow et al is a good read, in particular sections 5.1-5.5 and 5.7-5.11_. These sections summarize neatly what we have done till now and point to what is coming with respect to deep learning. -!split -===== Thursday September 29 ===== - -"Overview Video, why do we care about gradient methods?":"https://www.uio.no/studier/emner/matnat/fys/FYS-STK3155/h20/forelesningsvideoer/OverarchingAimsWeek39.mp4?vrtx=view-as-webpage" - !split @@ -1152,39 +1147,6 @@ minibatches. We denote these minibatches by $B_k$ where $k=1,\cdots,n/M$. -!split -===== Stochastic Gradient Descent ===== - -Stochastic gradient descent (SGD) and variants thereof address some of -the shortcomings of the Gradient descent method discussed above. - -The underlying idea of SGD comes from the observation that the cost -function, which we want to minimize, can almost always be written as a -sum over $n$ data points $\{\mathbf{x}_i\}_{i=1}^n$, -!bt -\[ -C(\mathbf{\beta}) = \sum_{i=1}^n c_i(\mathbf{x}_i, -\mathbf{\beta}). -\] -!et - -!split -===== Computation of gradients ===== - -This in turn means that the gradient can be -computed as a sum over $i$-gradients -!bt -\[ -\nabla_\beta C(\mathbf{\beta}) = \sum_i^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}). -\] -!et - -Stochasticity/randomness is introduced by only taking the -gradient on a subset of the data called minibatches. If there are $n$ -data points and the size of each minibatch is $M$, there will be $n/M$ -minibatches. We denote these minibatches by $B_k$ where -$k=1,\cdots,n/M$. !split ===== SGD example ===== @@ -1273,7 +1235,16 @@ 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. +reasonable time such that we do not move at all. Such approaches are +also called scaling. There are many such ways to "scale the learning +rate":"https://towardsdatascience.com/gradient-descent-the-learning-rate-and-the-importance-of-feature-scaling-6c0b416596e1" +and "discussions here":"https://www.jmlr.org/papers/volume23/20-1258/20-1258.pdf". See +also +URL:"https://towardsdatascience.com/learning-rate-schedules-and-adaptive-learning-rate-methods-for-deep-learning-2c8f433990d1" +for a discussion of different scaling functions for the learning rate. + +!split +===== Time decay rate ===== As an example, let $e = 0,1,2,3,\cdots$ denote the current epoch and let $t_0, t_1 > 0$ be two fixed numbers. Furthermore, let $t = e \cdot m + i$ where $m$ is the number of minibatches and $i=0,\cdots,m-1$. Then the function $$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length $\gamma_j (0; t_0, t_1) = t_0/t_1$ which decays in *time* $t$. @@ -1314,81 +1285,6 @@ print("gamma_j after %d epochs: %g" % (n_epochs,gamma_j)) -!split -===== Program for stochastic gradient ===== - -!bc pycod -# 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) - -X = np.c_[np.ones((m,1)), 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_) - - -theta = np.random.randn(2,1) -eta = 0.1 -Niterations = 1000 - - -for iter in range(Niterations): - gradients = 2.0/m*X.T @ ((X @ theta)-y) - 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) - - -n_epochs = 50 -t0, t1 = 5, 50 -def learning_schedule(t): - return t0/(t+t1) - -theta = np.random.randn(2,1) - -# note: here the number of minibatches is equal to the number of points!! -for epoch in range(n_epochs): - 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) - 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') -plt.axis([0,2.0,0, 15.0]) -plt.xlabel(r'$x$') -plt.ylabel(r'$y$') -plt.title(r'Random numbers ') -plt.show() - -!ec - -_Challenge_: try to write a similar code for a Logistic Regression case. - - !split ===== Code with a Number of Minibatches which varies ===== @@ -1435,6 +1331,7 @@ 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) @@ -1465,214 +1362,6 @@ plt.show() - - - -!split -===== SGD example ===== - -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 -$B_1 = (\mathbf{x}_1,\mathbf{x}_2), \cdots, B_5 = -(\mathbf{x}_9,\mathbf{x}_{10})$. Note that if you choose $M=1$ you -have only a single batch with all data points and on the other extreme, -you may choose $M=n$ resulting in a minibatch for each datapoint, i.e -$B_k = \mathbf{x}_k$. - -The idea is now to approximate the gradient by replacing the sum over -all data points with a sum over the data points in one the minibatches -picked at random in each gradient descent step -!bt -\[ -\nabla_{\beta} -C(\mathbf{\beta}) = \sum_{i=1}^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}) \rightarrow \sum_{i \in B_k}^n \nabla_\beta -c_i(\mathbf{x}_i, \mathbf{\beta}). -\] -!et - -!split -===== The gradient step ===== - -Thus a gradient descent step now looks like -!bt -\[ -\beta_{j+1} = \beta_j - \gamma_j \sum_{i \in B_k}^n \nabla_\beta c_i(\mathbf{x}_i, -\mathbf{\beta}) -\] -!et - -where $k$ is picked at random with equal -probability from $[1,n/M]$. An iteration over the number of -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. - -!split -===== Simple example code ===== - -!bc pycod -import numpy as np - -n = 100 #100 datapoints -M = 5 #size of each mini-batche -m = int(n/M) #number of minibatches -n_epochs = 10 #number of epochs - -j = 0 -for epoch in range(1,n_epochs+1): - for i in range(m): - k = np.random.randint(m) #Pick the k-th minibatch at random - #Compute the gradient using the data in minibatch Bk - #Compute new suggestion for - j += 1 -!ec - -Taking the gradient only on a subset of the data has two important -benefits. First, it introduces randomness which decreases the chance -that our opmization scheme gets stuck in a local minima. Second, if -the size of the minibatches are small relative to the number of -datapoints ($M < n$), the computation of the gradient is much -cheaper since we sum over the datapoints in the $k-th$ minibatch and not -all $n$ datapoints. - -!split -===== When do we stop? ===== - -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 -threshold and stop if true. However, the condition that the gradient -is zero is valid also for local minima, so this would only tell us -that we are close to a local/global minimum. However, we could also -evaluate the cost function at this point, store the result and -continue the search. If the test kicks in at a later stage we can -compare the values of the cost function and keep the $\beta$ that -gave the lowest value. - -!split -===== Slightly different approach ===== - -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. - -As an example, let $e = 0,1,2,3,\cdots$ denote the current epoch and let $t_0, t_1 > 0$ be two fixed numbers. Furthermore, let $t = e \cdot m + i$ where $m$ is the number of minibatches and $i=0,\cdots,m-1$. Then the function $$\gamma_j(t; t_0, t_1) = \frac{t_0}{t+t_1} $$ goes to zero as the number of epochs gets large. I.e. we start with a step length $\gamma_j (0; t_0, t_1) = t_0/t_1$ which decays in *time* $t$. - -In this way we can fix the number of epochs, compute $\beta$ and -evaluate the cost function at the end. Repeating the computation will -give a different result since the scheme is random by design. Then we -pick the final $\beta$ that gives the lowest value of the cost -function. - -!bc pycod -import numpy as np - -def step_length(t,t0,t1): - return t0/(t+t1) - -n = 100 #100 datapoints -M = 5 #size of each minibatch -m = int(n/M) #number of minibatches -n_epochs = 500 #number of epochs -t0 = 1.0 -t1 = 10 - -gamma_j = t0/t1 -j = 0 -for epoch in range(1,n_epochs+1): - for i in range(m): - k = np.random.randint(m) #Pick the k-th minibatch at random - #Compute the gradient using the data in minibatch Bk - #Compute new suggestion for beta - t = epoch*m+i - gamma_j = step_length(t,t0,t1) - j += 1 - -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$. - - -!split -===== 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 - -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.inv(X.T @ 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 - - -for iter in range(Niterations): - gradients = 2.0/n*X.T @ ((X @ theta)-y) - 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) - -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 = (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') -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 ===== Replace or not ===== @@ -2468,3 +2157,26 @@ print("Trained loss:", training_loss(weights)) !ec +!split +===== Introducing "JAX":"https://jax.readthedocs.io/en/latest/" ===== + +Presently, instead of using _autograd_, we recommend using "JAX":"https://jax.readthedocs.io/en/latest/" + +_JAX_ is Autograd and "XLA (Accelerated Linear Algebra))":"https://www.tensorflow.org/xla", +brought together for high-performance numerical computing and machine learning research. +It provides composable transformations of Python+NumPy programs: differentiate, vectorize, parallelize, Just-In-Time compile to GPU/TPU, and more. + +Here's a simple example on how you can use _JAX_ to compute the derivate of the logistic function. + +!bc pycod +import jax.numpy as jnp +from jax import grad, jit, vmap + +def sum_logistic(x): + return jnp.sum(1.0 / (1.0 + jnp.exp(-x))) + +x_small = jnp.arange(3.) +derivative_fn = grad(sum_logistic) +print(derivative_fn(x_small)) + +!ec