0.3.5 ADAM

The ADAM method combines momentum with a step-size adapted to each coordinate of the iterates \(x_{n}.\) The latter is somewhat like the Newton-Raphson method, but ADAM method uses only the gradient of the objective function, not the Hessian. ADAM is an evolution of its predecessors AdaGrad and RMSProp, and is best understood its this historical context.

The AdaGrad (adaptive gradient) method keeps track of the magnitude of the gradients in each dimension through the sums \[v_{n,i}=\sum_{m=0}^{n}(\partial f_{i}(x_{m}))^{2},\quad i=1,\ldots,d,n\ge0.\tag{3.51}\] They are used to adapt the step-size using the update rule \[x_{n+1}=x_{n}-\eta\nabla f(x_{n})\odot v_{n}^{-1/2}\tag{3.52}\] where the inverse square root is taken component-wise, and \(\odot\) denotes component-wise multiplication. The effect of AdaGrad is to decrease the effective step-size in directions with large gradients, which should increase stability, while increasing it in directions with small gradients, which should speed up convergence. Note that (3.51) sums over the entire past, so in principle \(v_{n,i}\) can grow without bound with increasing \(n,\) though if the iterates get close to a local minimum the derivatives should become small, thus putting a bound on this growth. Furthermore, a few very large derivatives early on in the iteration has an effect which lasts for all \(n\ge0.\)

RMSProp (Root Mean Square Propagation) adds a “forgetting mechanism” to AdaGrad, which emphasizes recent derivatives in the computation of the effective step-size. For a parameter \(\gamma\in(0,1)\) it defines the vectors \(v_{n}\in\mathbb{R}^{d}\) by \[v_{n,i}=\gamma v_{n-1,i}+(1-\gamma)\partial_{i}f(x_{n})^{2},\quad i=1,\ldots,d,n\ge0,\tag{3.53}\] which means that (3.51) is replaced by \[v_{n,i}=(1-\gamma)\sum_{m=0}^{n}\gamma^{n-m}(\partial f_{i}(x_{n}))^{2},\quad i=1,\ldots,d,n\ge0.\] The update rule for \(x_{n}\) of RMSProp is the same as for AdaGrad, i.e. (3.52).

Finally ADAM combines the mechanism of RMSProp with momentum. The vector \(v_{n}\) is still defined as for RMSProp in (3.53). A normalized version \(\hat{v}_{n}\) is defined by \[\hat{v}_{n,i}=\frac{v_{n,i}}{1-\gamma^{n+1}}=\frac{\sum_{m=0}^{n}\gamma^{n-m}(\partial f_{i}(x_{n}))^{2}}{\sum_{m=0}^{n}\gamma^{n-m}},\quad i=1,\ldots,d,n\ge0\] The number \(\hat{v}_{n,i}\) can be seen as a temporally smoothed estimate of the square derivative \(\partial_{i}f(x_{n})^{2}.\) The momentum mechanism is implemented using the iteration \[\mu_{n,i}=\beta\mu_{n-1,i}+(1-\beta)\partial_{i}f(x_{n}),\quad i=1,\ldots,d,n\ge0,\] for which \[\mu_{n,i}=(1-\beta)\sum_{m=0}^{n}\beta^{n-m}\partial_{i}f(x_{n}),\quad i=1,\ldots,d,n\ge0.\] Its normalized version \[\hat{\mu}_{n,i}=\frac{\mu_{n,i}}{1-\beta^{n+1}}=\frac{\sum_{m=0}^{n}\beta^{n-m}\partial_{i}f(x_{n})}{\sum_{m=0}^{n}\beta^{n-m}},\quad i=1,\ldots,d,n\ge0,\] can be seen as a temporally smoothed estimate of the derivative \(\partial_{i}f(x_{n}).\) The update rule of ADAM is then \[x_{n+1}=x_{n}-\eta\odot\frac{\hat{\mu}_{n}}{\sqrt{\hat{v}_{n}}+\varepsilon},\quad n\ge0,\tag{3.54}\] where all operations are component-wise, and \(\varepsilon>0\) is used to prevent numerical issues that may arise from the division if \(\hat{v}_{n}\) has very small entries. This update uses the temporally smoothed estimates of \(\partial f_{i}(x_{n})\) to determine the direction of the update, which has a similar effect has momentum for GD, and uses essentially the same adaptive step-size mechanism based on an estimate of the a smoother estimate of \(\partial_{i}f(x_{n})^{2}\) as RMSProp does. Typical parameter values are \(\varepsilon=10^{-8},\beta=0.9,\gamma=0.999.\)

When training artificial neural networks in modern machine learning one essentially always uses either GD with momentum or ADAM (or rather their stochastic variants - see the next section).

0.3.6 Stochastic Gradient Descent

Stochastic Gradient Descent (SGD) is a version of gradient descent for optimizing (the expectation of a) random objective function. One assumes that \(f_{0},f_{1},f_{2},\ldots\) is a sequence of random objective functions from \(\mathbb{R}^{d}\) to \(\mathbb{R}\) which are identically distributed, mutually independent (i.e. i.i.d.) and almost surely differentiable (i.e. \(\mathbb{P}(f_{m}\text{ is differentiable})=1\)). The goal is to minimize the expected objective function, i.e. to solve the optimization problem \[\text{minimize }\quad\mathbb{E}[f_{m}(x)]\quad\text{for}\quad x\in\mathbb{R}^{d},\tag{3.55}\] (the expectation is independent of \(m\) since the \(f_{m}\) are identically distributed). The update rule is as for gradient descent but using the random gradient \(\nabla f_{n}\) is the \(n\)-th step, i.e. \[x_{n+1}=x_{n}-\eta\nabla f_{n}(x_{n}),\quad n\ge0.\tag{3.56}\]

0.3.6.1 Motivation

The SGD algorithm is extremely important in modern machine learning (ML), where it is used to train e.g. artificial neural networks. Some context is needed to understand how one ends up with a sequence of random objective functions in this setting. We will study this in detail in the ML part of the course, but anticipate the most relevant aspects here. In ML (more specifically so called “supervised learning”) the learning step involves minimizing an (non-random) objective function \(\overline{f}(x),\) where \(x\in\mathbb{R}^{d}\) is a vector of parameters of the “model” (like the parameters \(a,b\) of the linear regression in Figure 2), which depends on a dataset \(z_{i},i=1,\ldots,p\) of vectors in \(\mathbb{R}^{\nu},\) for some \(\nu\ge1.\) More specifically, this objective function takes the form \[\overline{f}(x)=\frac{1}{p}\sum_{i=1}^{p}L(z_{i},x),\] where \(L:\mathbb{R}^{\nu}\times\mathbb{R}^{d}\to\mathbb{R}\) is a so called “loss function”, which is differentiable and such that \(L(z_{i},x)\) measures how well the model with parameter vector \(x\) fits the datapoint \(z_{i}\) (for linear regression this is the squared distance between the line with parameters \(x=(a,b)\) and the datapoint \(z_{i}=(\tilde{x}_{i},y_{i})\)). One can apply gradient descent to minimize this objective function, with updates \[x_{n+1}=x_{n}-\eta\frac{1}{p}\sum_{i=1}^{p}\nabla_{x}L(z_{i},x).\] For a large dataset, say \(p=10^{6},\) one step of GD requires the computation of the gradient \(\nabla_{x}L(z_{i},x)\) of the loss function \(p=10^{6}\) times. This is a lot of computation for just taking one step of gradient descent. To speed up the computation, it is common to instead pick just one random datapoint out of the dataset and take a step based only on the gradient of the loss function for that datapoint, i.e. \[x_{n+1}=x_{n}-\eta\nabla_{x}L(z_{i_{n}},x),\tag{3.57}\] where \(i_{n},n\ge0,\) is a i.i.d. sequence of independent integer random variables that are uniformly distributed on \(\{1,\ldots,p\}.\) Note that \[\mathbb{E}\left[L(z_{i_{m}},x)\right]=\frac{1}{p}\sum_{i=1}^{p}L(z_{i},x)=\overline{f}(x),\] and \[\mathbb{E}\left[\nabla_{x}L(z_{i_{m}},x)\right]=\frac{1}{p}\sum_{i=1}^{p}\nabla_{x}L(z_{i},x)=\nabla\overline{f}(x).\] Thus the direction of the update in (3.57) is \[\underset{\text{gradient of non-random obj. func.}}{\underbrace{\nabla\overline{f}(x)}}+\underset{\text{random noise with mean zero}}{\underbrace{E_{n}}}\] for \[E_{n}=\nabla_{x}L(z_{i_{m}},x)-\nabla\overline{f}(x).\] The sequence \(E_{n},n\ge0,\) is a i.i.d. sequence of mean-zero random noise that perturbs the the descent direction \(-\nabla\overline{f}(x)\) for the non-random objective function. Thus (3.57) amounts to a noisy version of GD. It is also an instance of the SGD algorithm as defined by the update rule (3.56), for \[f_{n}(x)=L(z_{i_{n}},x),\quad n\ge0.\]

The aforementioned motivation for SGD is computational efficiency. Another possible motivation for SGD is to improve stability or speed up convergence. If an objective function \(f:\mathbb{R}^{d}\to\mathbb{R}\) is very rough it may be unsuitable for GD, for instance if it has many non-global local minima that can trap GD. Adding a suitable random noise \(g(x)\) to \(f(x)\) at each point \(x\) and using SGD may allow the stochasticity of the iterates of SGD to avoid such “spurious” local minima. It may also encourage more exploration of the possible descent directions, and therefore the discovery of better directions that GD misses.

0.3.6.2 (Lack of) convergence and step-size schedules

An obvious difference between SGD and GD is that for SGD the iterates \(x_{n}\) are random variables. A very important difference between SGD and GD is that SGD with constant step-size (and non-trivial noise) never converges. Indeed, imagine starting GD at a global minimizer \(x_{*}\): then the iterates \(x_{n}\) all equal \(x_{*}\) and GD stays at \(x_{*}\) forever. But if SGD is started at \(x_{*}\) and \(\nabla f_{n}(x_{*})\) has non-trivial variance, then even though \(x_{0}=x_{*}\) the next iterate \(x_{1}\) will jump away from \(x_{*}.\) Instead of converging, SGD with constant step-size will in the best case settle into a region of a minimizer \(x_{*},\) where it will jump around randomly. This is illustrated in the following figure.

In such a scenario, the iterates \(x_{n}\) will converge in distribution as \(x_{n}\to\infty\) to some limiting random variable \(x_{\infty},\) whose variance corresponds to the typical distance between \(x_{n}\) and \(x_{*}\) (for \(n\) large). This variance will depend on the variance of the increments \(-\eta\nabla f_{n}(x)\) in the neighborhood of \(x_{*},\) which in turn depend on the variance of \(\nabla f_{n}(x)\) and on \(\eta>0.\) For smaller \(\eta>0\) the variance of the limiting distribution will be smaller, and the region where \(x_{n}\) settles down will be smaller, so \(x_{n}\) will typically be closer to \(x_{*}.\) But it will never converge to \(x_{*}.\)

If convergence is desired, one can achieve it with a step-size schedule, i.e. an \(n\)-dependent decreasing step-size. The rule (3.56) is then replaced with \[x_{n+1}=x_{n}-\eta_{n}\nabla f_{n}(x_{n}),\quad n\ge0,\] for a decreasing sequence \(\eta_{n},n\ge0.\) One can for instance take \(\eta_{n}\) proportional to \(n^{-1}.\) Below we will prove the convergence of SGD for such a schedule and strongly convex functions \(f_{n}\) with Lipschitz gradients.

In practice, a common approach is to start with some initial step-size that does not cause divergence, e.g. \(\eta_{0}=\eta=1,\) and then adaptively decrease \(\eta\) once a certain condition is met. For instance, if the objective function value \(f_{n}(x_{n})\) averaged over the last \(100\) steps did not decrease compared to average over the \(100\) steps before that, then \(\eta\) is decreased by half.

0.3.6.3 Momentum and ADAM

Momentum and ADAM generalize straight-forwardly from GD to SGD (in fact ADAM was introduces as a variant of SGD). By replacing the deterministic gradient \(\nabla f\) with the random gradient \(\nabla f_{n},\) the rule (3.35) for GD with momentum turns into the rule \[x_{n+1}=x_{n}-\eta\nabla f_{n}(x_{n})+\beta(x_{n}-x_{n-1}),\quad n\ge0,\] for SGD with momentum. Typically one uses \(\beta=0.9.\)

The ADAM rules for GD from (3.53)-(3.54) turn into \[v_{n,i}=\gamma v_{n-1,i}+(1-\gamma)\partial_{i}f_{n}(x_{n})^{2},\quad i=1,\ldots,d,n\ge0,\tag{3.64}\] \[\hat{v}_{n,i}=\frac{v_{n,i}}{1-\gamma^{n+1}},\quad i=1,\ldots,d,n\ge0,\] \[\mu_{n,i}=\beta\mu_{n-1,i}+(1-\beta)\partial_{i}f_{n}(x_{n}),\quad i=1,\ldots,d,n\ge0,\] \[\hat{\mu}_{n,i}=\frac{\mu_{n,i}}{1-\beta^{n+1}},\quad i=1,\ldots,d,n\ge0,\] \[x_{n+1}=x_{n}-\eta\odot\frac{\hat{\mu}_{n}}{\sqrt{\hat{v}_{n}}+\varepsilon},\quad n\ge0,\tag{3.65}\] with typical parameter values \(\varepsilon=10^{-8},\beta=0.9,\gamma=0.999.\)

These can of course also be combined with a step-size schedule by replacing \(\eta\) by an \(n\)-dependent \(\eta_{n}.\)

The rest of this section is concerned with the theoretical analysis of SGD.

0.3.6.4 Theoretical analysis

We start by investigating SGD it in the simplest case possible, namely for a univariate (\(d=1\)) quadratic objective function with Gaussian noise.

Next we consider SGD and the same objective function, but with a decreasing step-size schedule.

Next we consider a higher dimensional quadratic objective function with isotropic Gaussian noise. Isotropic means invariant to direction. For a \(\mathbb{R}^{d}\)-valued Gaussian random vector \(\varepsilon=(\varepsilon_{1},\ldots,\varepsilon_{d})\) it means that the entries \(\varepsilon_{i}\) are i.i.d. Gaussians with mean zero. Then the inner products \(\varepsilon\cdot e_{i}=\varepsilon_{i}\) with any of the standard basis vectors \(e_{1},\ldots,e_{d}\) of \(\mathbb{R}^{d}\) are Gaussians with the same variance and mean zero. More generally, for any fixed unit vector \(v\in\mathbb{R}^{d},\) the inner product \(\varepsilon\cdot v\) is Gaussian with that variance and mean zero. This follows from the fact that any linear combination of jointly Gaussian random variables is again Gaussian, and the following easy computation of the mean and variance of \(\varepsilon\cdot v\): \[\mathbb{E}[\varepsilon\cdot v]=\mathbb{E}\left[\sum_{i=1}^{d}\varepsilon_{i}v_{i}\right]=\sum_{i=1}^{d}\mathbb{E}[\varepsilon_{i}]v_{i}=\sum_{i=1}^{d}0\cdot v_{i}=0,\] (mean) and using (3.67) \[\mathbb{E}[(\varepsilon\cdot v)^{2}]=\mathbb{E}\left[\left(\sum_{i=1}^{d}\varepsilon_{i}v_{i}\right)^{2}\right]=\mathbb{E}[\varepsilon_{1}]^{2}\sum_{i,j=1}^{d}v_{i}^{2}=\mathbb{E}[\varepsilon_{1}]^{2},\] (variance). Thus the distribution of the inner product \(\varepsilon\cdot v\) is independent of the direction of \(v,\) which is an expression of the isotropy of the distribution of \(\varepsilon.\) Moreover, if \(w_{1},\ldots,w_{d}\) is any orthonormal basis of \(\mathbb{R}^{d},\) then the inner products \(\varepsilon\cdot w_{1},\ldots,\varepsilon\cdot w_{n}\) are independent. One can verify this by showing that the covariances \(\mathbb{E}[(\varepsilon\cdot w_{i})(\varepsilon\cdot w_{j})]\) are zero for \(i\ne j,\) using a computation similar to (3.67). This means that the the joint probability distribution of the coordinates of \(\varepsilon\) is the same regardless of the orthonormal basis used to compute them. This is an even stronger expression of the isotropy of \(\varepsilon.\)

Appendix

0.3.7 Proof of Lagrange Multiplier Theorem

0.3.8 Additional code for gradient descent Example 3.2

For a complicated objective function, it is easy to make a mistake when computing the derivatives, or when implementing the computation in code. This can lead to strange behavior, such as the objective function increasing as GD progresses. To find such errors it is convenient to compare the output of ones gradient computation routines with numerically computed approximate gradients. The code below gives an example of this.

import numpy as np

def g(a,b,x):
    return b*np.exp(-0.5*(x-a)**2)

# Gradient of g wrt to a
def grad_g_a(a,b,x):
    return (x-a)*g(a,b,x)

# Gradient of g wrt to b
def grad_g_b(a,b,x):
    return np.exp(-0.5*(x-a)**2)

def h(avec,bvec,x):
    ret = 0
    for i in range(len(avec)):
        ret += g(avec[i],bvec[i],x)
    return ret

# Gradient of h
def grad_h(avec,bvec,x):
    reta = np.zeros(len(avec))
    retb = np.zeros(len(bvec))
    for i in range(len(avec)):
        reta[i] = grad_g_a(avec[i],bvec[i],x)
        retb[i] = grad_g_b(avec[i],bvec[i],x)
    return (reta,retb)

def L(avec,bvec,xvec,yvec):
    ret = 0
    for k in range(len(xvec)):
        ret += (h(avec,bvec,xvec[k])-yvec[k])**2/2.0
    return ret

# Gradient of L
def grad_L(avec,bvec,xvec,yvec):
    reta = np.zeros(len(avec))
    retb = np.zeros(len(bvec))
    for k in range(len(xvec)):
        (ga,gb) = grad_h(avec,bvec,xvec[k])
        reta += (h(avec,bvec,xvec[k])-yvec[k])*ga
        retb += (h(avec,bvec,xvec[k])-yvec[k])*gb
    return (reta,retb)


xvec = np.array( [-1,0,1] )
yvec = xvec


# Code to debug gradient calculation:
# compare to numerically computed derivatives
avec = np.random.uniform( -1, 1, num_bumps )
bvec = np.random.uniform( -1, 1, num_bumps )
evec = np.array( [0,1,0,0])
eps = 0.00001
# numerically computed derivative wrt to a_1
print( (h(avec+eps*evec,bvec,xvec[0])-h(avec,bvec,xvec[0]))/eps )
# numerically computed derivative wrt to b_2
print( (h(avec,bvec+eps*evec,xvec[0])-h(avec,bvec,xvec[0]))/eps ) 
# exact gradient
print( grad_h(avec,bvec,xvec[0]) )

The following code can be used to create plots based on the output of gradient descent.

import matplotlib.pyplot as plt

# Data points
xvec = np.array( [-1,0,1] )
yvec = xvec

# Initial parameters
avec0 = np.random.uniform(-1, 1, num_bumps)
bvec0 = np.zeros(num_bumps)

# Run GD
(histL,hista,histb) = GD_minimize(avec0,bvec0, \
                                 xvec,yvec,num_steps)



# Plot results

plt.plot(histL)
plt.title("Objective value $L(a_n,b_n)$ for each step of GD")
plt.xlabel("Step")
plt.ylabel("$L$")
plt.show()

plt.semilogy(histL)
plt.title("Objective value $L(a_n,b_n)$ for each step of GD")
plt.xlabel("Step")
plt.ylabel("$L$ (log scale)")
plt.show()


xs = np.linspace( -2, 2, 100 )
for n in [0,1,4,10,len(hista)-1]:
    plt.plot( xs,
             [h(hista[n],histb[n],x) for x in xs],
             label=f"$n={n}$" )
plt.scatter(xvec, yvec,color='black')
plt.title( "Function $x \\to h(a_n,b_n,x)$" + \
           "encoded by values $(a_n,b_n)$ at $n$-th step of GD" )
plt.show()

for i in range(4):    
    colors = ['blue', 'red', 'green', 'purple']
    plt.scatter([hista[0][i]],
                [histb[0][i]],
                marker='x', color=colors[i])
    plt.scatter([hista[-1][i]],
                [histb[-1][i]],
                marker='o', color=colors[i])
    plt.plot([avec[i] for avec in hista],
             [bvec[i] for bvec in histb],
             color=colors[i], label=f"$i={i+1}$")
last_n = len(hista)-1
plt.title(f"Trajectories of $(a_{{n,i}},b_{{n,i}})$" + \
          f"in plane throughout optimization\n" + \
          f"X: starting point $(a_{{0,i}},b_{{0,i}}),$" + \
          f"Dot: ending point " + \
          f"$(a_{{{last_n},i}},b_{{{last_n},i}})$")
plt.show()

0.3.9 SGD for non-convex function

The best we can do is the following proposition, which guarantees that the gradient \(\nabla f(x_{n})\) decays under appropriate conditions. For technical reasons it actually bounds the minimum of \(|\nabla f(x_{m})|^{2}\) over \(m=0,\ldots,n.\)

Home

Contents

Study Weeks

✕