Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Lecture 02, Notebook 02: Gradient Descent and Stochastic Gradient Descent

University of Lausanne

Course: Deep Learning for Solving and Estimating Dynamic Models in Economics and Finance
Script reference: §1.5 (Optimization: gradient descent, SGD, mini-batch, learning rate)
Notebook role: core
Author: Simon Scheidegger

RUN_MODE = "smoke"  # one of: "smoke", "teaching", "production"
SEED = 0
%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import random
from scipy import stats
from scipy.optimize import fmin

# Reproducibility: fix seeds for the SGD examples below.
np.random.seed(SEED)
random.seed(0)

Gradient Descent

  • Gradient descent, also known as steepest descent, is an optimization algorithm for finding the local minimum of a function.

  • To find a local minimum, the function “steps” in the direction of the negative of the gradient.

  • Gradient ascent is the same as gradient descent, except that it steps in the direction of the positive of the gradient and therefore finds local maxima instead of minima. The algorithm of gradient descent can be outlined as follows:

    1:   Choose initial guess x0x_0
    2:   for k = 0, 1, 2, ... do
    3:       sks_k = -f(xk)\nabla f(x_k)
    4:       choose αk\alpha_k to minimize f(xk+αksk)f(x_k+\alpha_k s_k)
    5:       xk+1=xk+αkskx_{k+1} = x_k + \alpha_k s_k
    6:   end for

As an analytical example, we try to find a local minimum for the function f(x)=x32x2+2f(x) = x^3-2x^2+2

f = lambda x: x**3-2*x**2+2
x = np.linspace(-1,2.5,1000)
plt.plot(x,f(x))
plt.xlim([-1,2.5])
plt.ylim([0,3])
plt.show()
<Figure size 432x288 with 1 Axes>

We can see from the graph above that our local minimum is gonna be around 1.4 (on the x-axis).

Let’s pretend that we don’t know that, so we set our starting point (arbitrarily, in this case) at x0=2x_0 = 2

x_old = 0
x_new = 2.0 # The algorithm starts at x=2
n_k = 0.1 # step size
precision = 0.0001

x_list, y_list = [x_new], [f(x_new)]

# returns the value of the derivative of our function
def f_prime(x):
    return 3*x**2-4*x
 
while abs(x_new - x_old) > precision:
    x_old = x_new
    s_k = -f_prime(x_old)
    x_new = x_old + n_k * s_k
    x_list.append(x_new)
    y_list.append(f(x_new))
print("Local minimum occurs at:", x_new)
print("Number of steps:", len(x_list))
Local minimum occurs at: 1.3334253508453249
Number of steps: 17

The figures below show step-by-step the route that was taken to find the local minimum.

plt.figure(figsize=[10,3])
plt.subplot(1,2,1)
plt.scatter(x_list,y_list,c="r")
plt.plot(x_list,y_list,c="r")
plt.plot(x,f(x), c="b")
plt.xlim([-1,2.5])
plt.ylim([0,3])
plt.title("Gradient descent with iteration step")
plt.subplot(1,2,2)
plt.scatter(x_list,y_list,c="r")
plt.plot(x_list,y_list,c="r")
plt.plot(x,f(x), c="b")
plt.xlim([1.2,2.1])
plt.ylim([0,3])
plt.title("Gradient descent -- zoomed-in")
plt.show()
<Figure size 720x216 with 2 Axes>

Note: that the step size (also called learning rate) in the implementation above is constant.

->Doing this makes it easier to implement the algorithm. 

However, it also presents some issues:

  • If the step size is too small, then convergence might be very slow.

  • If we make it too large, then the method may fail to converge at all.

    -> A solution to this is to use adaptive step sizes as the algorithm below does (using scipy’s fmin function to find optimal step sizes)

# we setup this function to pass into the fmin algorithm
def f2(n,x,s):
    x = x + n*s
    return f(x)

x_old = 0
x_new = 2 # The algorithm starts at x=2
precision = 0.0001

x_list, y_list = [x_new], [f(x_new)]

# returns the value of the derivative of our function
def f_prime(x):
    return 3*x**2-4*x

while abs(x_new - x_old) > precision:
    x_old = x_new
    s_k = -f_prime(x_old)
    
    # use scipy fmin function to find ideal step size.
    # float(...) keeps x_new a plain scalar across scipy versions
    # (newer scipy returns a 0-d array, which would make x_list inhomogeneous).
    n_k = float(fmin(f2,0.1,(x_old,s_k), full_output = False, disp = False)[0])

    x_new = x_old + n_k * s_k
    x_list.append(x_new)
    y_list.append(f(x_new))
    
print("Local minimum occurs at ", float(x_new))
print("Number of steps:", len(x_list))
Local minimum occurs at  1.3333333284505209
Number of steps: 4
  • With adaptive step sizes, the algorithm converges in just 4 iterations rather than 17.

  • Note: it also takes extra time to compute the appropriate step size at each iteration.

Below are some plots of the path taken below.

You can see that it converges very quickly to a point near the local minimum, so it’s hard to even discern the dots after the first two steps until we zoom in very close in the third frame below.

plt.figure(figsize=[15,3])
plt.subplot(1,3,1)
plt.scatter(x_list,y_list,c="r")
plt.plot(x_list,y_list,c="r")
plt.plot(x,f(x), c="b")
plt.xlim([-1,2.5])
plt.title("Gradient descent -- iteration steps")
plt.subplot(1,3,2)
plt.scatter(x_list,y_list,c="r")
plt.plot(x_list,y_list,c="r")
plt.plot(x,f(x), c="b")
plt.xlim([1.2,2.1])
plt.ylim([0,3])
plt.title("zoomed in -- iteration steps")
plt.subplot(1,3,3)
plt.scatter(x_list,y_list,c="r")
plt.plot(x_list,y_list,c="r")
plt.plot(x,f(x), c="b")
plt.xlim([1.3333,1.3335])
plt.ylim([0,3])
plt.title("zoomed in more -- iteration steps")
plt.show()
<Figure size 1080x216 with 3 Axes>

An alternative approach to update the step size is to choose a decrease constant dd that shrinks the step size over time: η(t+1)=η(t)/(1+t×d)\eta(t+1) = \eta(t) / (1+t \times d).

x_old = 0
x_new = 2  # The algorithm starts at x=2
n_k = 0.17 # step size
precision = 0.0001
t, d = 0, 1

x_list, y_list = [x_new], [f(x_new)]

# returns the value of the derivative of our function
def f_prime(x):
    return 3*x**2-4*x
 
while abs(x_new - x_old) > precision:
    x_old = x_new
    s_k = -f_prime(x_old)
    x_new = x_old + n_k * s_k
    x_list.append(x_new)
    y_list.append(f(x_new))
    n_k = n_k / (1 + t * d)
    t += 1

print("Local minimum occurs at:", x_new)
print("Number of steps:", len(x_list))
Local minimum occurs at: 1.3308506740900838
Number of steps: 6

A more complicated example

  • Consider a simple linear regression where we want to see how the temperature affects the noises made by crickets.

  • There is a data set of cricket chirp rates at various temperatures.

  • First we’ll load that data set in and plot it:

Data: cricket chirp rate vs temperature

SGD_data.txt is a 16-row, two-column CSV with no header. Column 1: chirps per second (striped ground cricket, Gryllus rubens).
Column 2: ambient temperature in degrees Fahrenheit.

This is a standard introductory dataset (originally from G. W. Pierce, The Songs of Insects, Harvard University Press, 1948) widely used to illustrate simple linear regression. It is included here as a minimal worked example for gradient descent: fitting heta0+heta1x heta_0 + heta_1 x to the chirp-temperature relationship.

#Load the dataset
data = np.loadtxt('SGD_data.txt', delimiter=',')
 
#Plot the data
plt.scatter(data[:, 0], data[:, 1], marker='o', c='b')
plt.title('cricket chirps versus temperature')
plt.xlabel('chirps/sec for striped ground crickets')
plt.ylabel('temperature in degrees Fahrenheit')
plt.xlim([13,21])
plt.ylim([65,95])
plt.show()
<Figure size 432x288 with 1 Axes>
  • Our goal is to find the parameters of a linear regression model---that is, the equation of the straight line

    hθ(x)=θ0+θ1xh_\theta(x) = \theta_0 + \theta_1 x

    that best fits our data points.

  • Thus, the function that we are trying to minimize in this case is:

    J(θ0,θ1)=1mi=1m(hθ(xi)yi)2J(\theta_0,\theta_1) = {1 \over m} \sum\limits_{i=1}^m (h_\theta(x_i)-y_i)^2

  • In this case, our gradient will be defined in two dimensions:

    θ0J(θ0,θ1)=1mi=1m(hθ(xi)yi)\frac{\partial}{\partial \theta_0} J(\theta_0,\theta_1) = \frac{1}{m} \sum\limits_{i=1}^m (h_\theta(x_i)-y_i)

    θ1J(θ0,θ1)=1mi=1m((hθ(xi)yi)xi)\frac{\partial}{\partial \theta_1} J(\theta_0,\theta_1) = \frac{1}{m} \sum\limits_{i=1}^m ((h_\theta(x_i)-y_i) \cdot x_i)

Below, we set up our function for h, J and the gradient:

h = lambda theta_0,theta_1,x: theta_0 + theta_1*x

def J(x,y,m,theta_0,theta_1):
    returnValue = 0
    for i in range(m):
        returnValue += (h(theta_0,theta_1,x[i])-y[i])**2
    returnValue = returnValue/(2*m)
    return returnValue

def grad_J(x,y,m,theta_0,theta_1):
    returnValue = np.array([0.,0.])
    for i in range(m):
        returnValue[0] += (h(theta_0,theta_1,x[i])-y[i])
        returnValue[1] += (h(theta_0,theta_1,x[i])-y[i])*x[i]
    returnValue = returnValue/(m)
    return returnValue

Now, we load our training data into the x and y variables:

x = data[:, 0]
y = data[:, 1]
m = len(x)

Now, we run the gradient descent algorithm (without adaptive step sizes in this example):

theta_old = np.array([0.,0.])
theta_new = np.array([1.,1.]) # The algorithm starts at [1,1]
n_k = 0.001 # step size
precision = 0.001
num_steps = 0
s_k = float("inf")

while np.linalg.norm(s_k) > precision:
    num_steps += 1
    theta_old = theta_new
    s_k = -grad_J(x,y,m,theta_old[0],theta_old[1])
    theta_new = theta_old + n_k * s_k

print("Local minimum occurs where:")
print("theta_0 =", theta_new[0])
print("theta_1 =", theta_new[1])
print("This took",num_steps,"steps to converge")
Local minimum occurs where:
theta_0 = 25.128552558595363
theta_1 = 3.297264756251897
This took 565859 steps to converge

For comparison, let’s get the actual values for θ0\theta_0 and θ1\theta_1:

actualvalues = sp.stats.linregress(x,y)
print("Actual values for theta are:")
print("theta_0 =", actualvalues.intercept)
print("theta_1 =", actualvalues.slope)
Actual values for theta are:
theta_0 = 25.232304983426026
theta_1 = 3.2910945679475647
  • We can see that the values are relatively close to the actual values (even though our method was pretty slow). If you look at the source code of linregress, it uses the convariance matrix of x and y to compute fastly. Below, you can see a plot of the line drawn with our theta values against the data:

xx = np.linspace(0,21,1000)
plt.scatter(data[:, 0], data[:, 1], marker='o', c='b')
plt.plot(xx,h(theta_new[0],theta_new[1],xx))
plt.xlim([13,21])
plt.ylim([65,95])
plt.title('cricket chirps versus temperature')
plt.xlabel('chirps/sec for striped ground crickets')
plt.ylabel('temperature in degrees Fahrenheit')
plt.show()
<Figure size 432x288 with 1 Axes>

Notice that in the method above we need to calculate the gradient in every step of our algorithm.

  • In the example with the crickets, this is not a big deal since there are only 15 data points.

  • But imagine that we had 1 million data points.

  • If this were the case, it would certainly make the method above far less efficient.

  • In machine learning, the algorithm above is often called batch gradient descent to contrast it withmini-batch gradient descent (which we will not go into here) and stochastic gradient descent.

Stochastic gradient descent

  • In batch gradient descent, we must look at every example in the entire training set on every step (in cases where a training set is used for gradient descent).

  • This can be quite slow if the training set is sufficiently large.

  • In stochastic gradient descent, we update our values after looking at each item in the training set, so that we can start making progress right away.

  • Recall the linear regression example above. In that example, we calculated the gradient for each of the two theta values as follows:

    θ0J(θ0,θ1)=1mi=1m(hθ(xi)yi)\frac{\partial}{\partial \theta_0} J(\theta_0,\theta_1) = \frac{1}{m} \sum\limits_{i=1}^m (h_\theta(x_i)-y_i)

    θ1J(θ0,θ1)=1mi=1m((hθ(xi)yi)xi)\frac{\partial}{\partial \theta_1} J(\theta_0,\theta_1) = \frac{1}{m} \sum\limits_{i=1}^m ((h_\theta(x_i)-y_i) \cdot x_i)

    Where hθ(x)=θ0+θ1xh_\theta(x) = \theta_0 + \theta_1 x and the per-example cost is 12(hθ(xi)yi)2\tfrac{1}{2}(h_\theta(x_i)-y_i)^2 (the 12\tfrac{1}{2} absorbs the chain-rule factor of 2, which is why no 2 appears in the gradient or in the SGD update below).

  • Then we followed this algorithm (where α\alpha was a non-adapting stepsize):

    1:   Choose initial guess x0x_0
    2:   for k = 0, 1, 2, ... do
    3:       sks_k = -f(xk)\nabla f(x_k)
    4:       xk+1=xk+αskx_{k+1} = x_k + \alpha s_k
    5:   end for

  • When the sample data had 15 data points as in the example above, calculating the gradient was not very costly.

  • But for very large data sets, this would not be the case. So instead, we consider a stochastic gradient descent algorithm for simple linear regression such as the following, where m is the size of the data set:

    1:   Randomly shuffle the data set
    2:   for k = 0, 1, 2, ... do
    3:       for i = 1 to m do
    4:            [θ0θ1]=[θ0θ1]α[(hθ(xi)yi)xi(hθ(xi)yi)]\begin{bmatrix} \theta_{0} \\ \theta_1 \\ \end{bmatrix}=\begin{bmatrix} \theta_0 \\ \theta_1 \\ \end{bmatrix}-\alpha\begin{bmatrix} (h_\theta(x_i)-y_i) \\ x_i(h_\theta(x_i)-y_i) \\ \end{bmatrix}
    5:       end for
    6:   end for

  • Typically, with stochastic gradient descent, you will run through the entire data set 1 to 10 times (see value for k in line 2 of the pseudocode above), depending on how fast the data is converging and how large the data set is.

  • With batch gradient descent, we must go through the entire data set before we make any progress. With this algorithm though, we can make progress right away and continue to make progress as we go through the data set. Therefore, stochastic gradient descent is often preferred when dealing with large data sets.

  • Unlike gradient descent, stochastic gradient descent will tend to oscillate near a minimum value rather than continuously getting closer. It may never actually converge to the minimum though.

  • One way around this is to slowly decrease the step size α\alpha as the algorithm runs. However, this is less common than using a fixed α\alpha.

An analytical example

Here we demonstrate the use of stochastic gradient descent for linear regression. In the example below, we’ll create a set of 500,000 points around the line y=4x+10+ϵy = 4x+10+\epsilon, for values of x between 0 and 100:

f = lambda x: x*4+10+np.random.randn(len(x))*10

x = np.random.random(500000)*100
y = f(x) 
m = len(y)
  • We first randomly shuffle around our dataset.

  • Note: in this example, this step isn’t strictly necessary since the data is already in a random order.

  • However, that obviously may not always be the case:

from random import shuffle

x_shuf = []
y_shuf = []
index_shuf =  list(range(len(x)))
shuffle(index_shuf)
for i in index_shuf:
    x_shuf.append(x[i])
    y_shuf.append(y[i])
  • Now we setup our h function and our cost function, which we will use to check how the value is improving.

h = lambda theta_0,theta_1,x: theta_0 + theta_1*x
cost = lambda theta_0,theta_1, x_i, y_i: 0.5*(h(theta_0,theta_1,x_i)-y_i)**2
  • Next, we run our stochastic gradient descent algorithm. To see it’s progress, we’ll take a cost measurement at every step.

  • Every 10,000 steps, we’ll get an average cost from the last 10,000 steps and then append that to our cost_list variable.

  • We will run through the entire list 10 times here:

theta_old = np.array([0.,0.])
theta_new = np.array([1.,1.]) # The algorithm starts at [1,1]
n_k = 0.000005 # step size

iter_num = 0
s_k = np.array([float("inf"),float("inf")])
sum_cost = 0
cost_list = []

for j in range(10):
    for i in range(m):
        iter_num += 1
        theta_old = theta_new
        s_k[0] = (h(theta_old[0],theta_old[1],x[i])-y[i])
        s_k[1] = (h(theta_old[0],theta_old[1],x[i])-y[i])*x[i]
        s_k = (-1)*s_k
        theta_new = theta_old + n_k * s_k
        sum_cost += cost(theta_old[0],theta_old[1],x[i],y[i])
        if (i+1) % 10000 == 0:
            cost_list.append(sum_cost/10000.0)
            sum_cost = 0   
            
print("Local minimum occurs where:")
print("theta_0 =", theta_new[0] )
print("theta_1 =", theta_new[1])
Local minimum occurs where:
theta_0 = 9.935073783110743
theta_1 = 3.9996685608275655

As one can see, the values for θ0\theta_0 and θ1\theta_1 are close to their true values of 10 and 4.

Next, we plot our cost versus the number of iterations. As you can see, the cost goes down quickly at first, but starts to level off as we go through more iterations:

iterations = np.arange(len(cost_list))*10000
plt.plot(iterations,cost_list)
plt.xlabel("iterations")
plt.ylabel("average cost")
plt.show()
<Figure size 432x288 with 1 Axes>