Stochastic gradient descent

As mentioned earlier, classical stochastic gradient descent is generally sensitive to the step size η. Fast convergence requires large step sizes but this may induce numerical instability. The problem can be largely solved[13] by considering implicit updates whereby the stochastic gradient is evaluated at the next iterate rather than the current one:

${\displaystyle w^{new}:=w^{old}-\eta \nabla Q_{i}(w^{new}).}$

This equation is implicit since ${\displaystyle w^{new}}$ appears on both sides of the equation. It is a stochastic form of the proximal gradient method since the update can also be written as:

${\displaystyle w^{new}:=\arg \min _{w}\{Q_{i}(w)+{\frac {1}{2\eta }}||w-w^{old}||^{2}\}.}$

As an example, consider least squares with features ${\displaystyle x_{1},\ldots ,x_{n}\in \mathbb {R} ^{p}}$ and observations ${\displaystyle y_{1},\ldots ,y_{n}\in \mathbb {R} }$. We wish to solve:

${\displaystyle \min _{w}\sum _{j=1}^{n}(y_{j}-x_{j}'w)^{2},}$

where ${\displaystyle x_{j}'w=x_{j1}w_{1}+x_{j,2}w_{2}+...+x_{j,p}w_{p}}$ indicates the inner product. Note that ${\displaystyle x}$ could have "1" as the first element to include an intercept. Classical stochastic gradient descent proceeds as follows:

${\displaystyle w^{new}=w^{old}+\eta (y_{i}-x_{i}'w^{old})x_{i}}$

where ${\displaystyle i}$ is uniformly sampled between 1 and ${\displaystyle n}$. Although theoretical convergence of this procedure happens under relatively mild assumptions, in practice the procedure can be quite unstable. In particular, when ${\displaystyle \eta }$ is misspecified so that ${\displaystyle I-\eta x_{i}x_{i}'}$ has large absolute eigenvalues with high probability, the procedure may diverge numerically within a few iterations. In contrast, implicit stochastic gradient descent (shortened as ISGD) can be solved in closed-form as:

${\displaystyle w^{new}=w^{old}+{\frac {\eta }{1+\eta ||x_{i}||^{2}}}(y_{i}-x_{i}'w^{old})x_{i}.}$

This procedure will remain numerically stable virtually for all ${\displaystyle \eta }$ as the step size is now normalized. Such comparison between classical and implicit stochastic gradient descent in the least squares problem is very similar to the comparison between least mean squares (LMS) and normalized least mean squares filter (NLMS).

Even though a closed-form solution for ISGD is only possible in least squares, the procedure can be efficiently implemented in a wide range of models. Specifically, suppose that ${\displaystyle Q_{i}(w)}$ depends on ${\displaystyle w}$ only through a linear combination with features ${\displaystyle x_{i}}$, so that we can write ${\displaystyle \nabla _{w}Q_{i}(w)=-q(x_{i}'w)x_{i}}$, where ${\displaystyle q()\in \mathbb {R} }$ may depend on ${\displaystyle x_{i},y_{i}}$ as well but not on ${\displaystyle w}$ except through ${\displaystyle x_{i}'w}$. Least squares obeys this rule, and so does logistic regression, and most generalized linear models. For instance, in least squares, ${\displaystyle q(x_{i}'w)=y_{i}-x_{i}'w}$, and in logistic regression ${\displaystyle q(x_{i}'w)=y_{i}-S(x_{i}'w)}$, where ${\displaystyle S(u)=e^{u}/(1+e^{u})}$ is the logistic function. In Poisson regression, ${\displaystyle q(x_{i}'w)=y_{i}-e^{x_{i}'w}}$, and so on.

In such settings, ISGD is simply implemented as follows. Let ${\displaystyle f(\xi )=\eta q(x_{i}'w^{old}+\xi ||x_{i}||^{2})}$, where ${\displaystyle \xi }$ is scalar. Then, ISGD is equivalent to:

${\displaystyle w^{new}=w^{old}+\xi ^{\ast }x_{i},~{\text{where}}~\xi ^{\ast }=f(\xi ^{\ast }).}$

The scaling factor ${\displaystyle \xi ^{\ast }\in \mathbb {R} }$ can be found through the bisection method since in most regular models, such as the aforementioned generalized linear models, function ${\displaystyle q()}$ is decreasing, and thus the search bounds for ${\displaystyle \xi ^{\ast }}$ are ${\displaystyle [\min(0,f(0)),\max(0,f(0))]}$.

Further proposals include the momentum method, which appeared in Rumelhart, Hinton and Williams' paper on backpropagation learning.[14] Stochastic gradient descent with momentum remembers the update Δ w at each iteration, and determines the next update as a linear combination of the gradient and the previous update:[15][16]

${\displaystyle \Delta w:=\alpha \Delta w-\eta \nabla Q_{i}(w)}$

${\displaystyle w:=w+\Delta w}$

that leads to:

${\displaystyle w:=w-\eta \nabla Q_{i}(w)+\alpha \Delta w}$

where the parameter ${\displaystyle w}$ which minimizes ${\displaystyle Q(w)}$ is to be estimated, ${\displaystyle \eta }$ is a step size (sometimes called the learning rate in machine learning) and ${\displaystyle \alpha }$ is an exponential decay factor between 0 and 1 that determines the relative contribution of the current gradient and earlier gradients to the weight change.

The name momentum stems from an analogy to momentum in physics: the weight vector ${\displaystyle w}$, thought of as a particle traveling through parameter space,[14] incurs acceleration from the gradient of the loss ("force"). Unlike in classical stochastic gradient descent, it tends to keep traveling in the same direction, preventing oscillations. Momentum has been used successfully by computer scientists in the training of artificial neural networks for several decades.[17]

Averaged stochastic gradient descent, invented independently by Ruppert and Polyak in the late 1980s, is ordinary stochastic gradient descent that records an average of its parameter vector over time. That is, the update is the same as for ordinary stochastic gradient descent, but the algorithm also keeps track of[18]

${\displaystyle {\bar {w}}={\frac {1}{t}}\sum _{i=0}^{t-1}w_{i}}$.

When optimization is done, this averaged parameter vector takes the place of w.

AdaGrad (for adaptive gradient algorithm) is a modified stochastic gradient descent algorithm with per-parameter learning rate, first published in 2011.[19] Informally, this increases the learning rate for sparser parameters and decreases the learning rate for ones that are less sparse. This strategy often improves convergence performance over standard stochastic gradient descent in settings where data is sparse and sparse parameters are more informative. Examples of such applications include natural language processing and image recognition.[19] It still has a base learning rate η, but this is multiplied with the elements of a vector {G_j,j} which is the diagonal of the outer product matrix

${\displaystyle G=\sum _{\tau =1}^{t}g_{\tau }g_{\tau }^{\mathsf {T}}}$

where ${\displaystyle g_{\tau }=\nabla Q_{i}(w)}$, the gradient, at iteration τ. The diagonal is given by

${\displaystyle G_{j,j}=\sum _{\tau =1}^{t}g_{\tau ,j}^{2}}$.

This vector is updated after every iteration. The formula for an update is now

${\displaystyle w:=w-\eta \,\mathrm {diag} (G)^{-{\frac {1}{2}}}\circ g}$[lower-alpha 1]

or, written as per-parameter updates,

${\displaystyle w_{j}:=w_{j}-{\frac {\eta }{\sqrt {G_{j,j}}}}g_{j}.}$

Each {G_(i,i)} gives rise to a scaling factor for the learning rate that applies to a single parameter w_i. Since the denominator in this factor, ${\displaystyle {\sqrt {G_{i}}}={\sqrt {\sum _{\tau =1}^{t}g_{\tau }^{2}}}}$ is the ℓ₂ norm of previous derivatives, extreme parameter updates get dampened, while parameters that get few or small updates receive higher learning rates.[17]

While designed for convex problems, AdaGrad has been successfully applied to non-convex optimization.[20]

RMSProp (for Root Mean Square Propagation) is also a method in which the learning rate is adapted for each of the parameters. The idea is to divide the learning rate for a weight by a running average of the magnitudes of recent gradients for that weight.[21] So, first the running average is calculated in terms of means square,

${\displaystyle v(w,t):=\gamma v(w,t-1)+(1-\gamma )(\nabla Q_{i}(w))^{2}}$

where ${\displaystyle \gamma }$ is the forgetting factor.

And the parameters are updated as,

${\displaystyle w:=w-{\frac {\eta }{\sqrt {v(w,t)}}}\nabla Q_{i}(w)}$

RMSProp has shown good adaptation of learning rate in different applications. RMSProp can be seen as a generalization of Rprop and is capable to work with mini-batches as well opposed to only full-batches.[22]

Adam[23] (short for Adaptive Moment Estimation) is an update to the RMSProp optimizer. In this optimization algorithm, running averages of both the gradients and the second moments of the gradients are used. Given parameters ${\displaystyle w^{(t)}}$ and a loss function ${\displaystyle L^{(t)}}$, where ${\displaystyle t}$ indexes the current training iteration (indexed at ${\displaystyle 0}$), Adam's parameter update is given by:

${\displaystyle m_{w}^{(t+1)}\leftarrow \beta _{1}m_{w}^{(t)}+(1-\beta _{1})\nabla _{w}L^{(t)}}$

${\displaystyle v_{w}^{(t+1)}\leftarrow \beta _{2}v_{w}^{(t)}+(1-\beta _{2})(\nabla _{w}L^{(t)})^{2}}$

${\displaystyle {\hat {m}}_{w}={\frac {m_{w}^{(t+1)}}{1-\beta _{1}^{t+1}}}}$

${\displaystyle {\hat {v}}_{w}={\frac {v_{w}^{(t+1)}}{1-\beta _{2}^{t+1}}}}$

${\displaystyle w^{(t+1)}\leftarrow w^{(t)}-\eta {\frac {{\hat {m}}_{w}}{{\sqrt {{\hat {v}}_{w}}}+\epsilon }}}$

where ${\displaystyle \epsilon }$ is a small scalar (e.g. ${\displaystyle 10^{-8}}$) used to prevent division by 0, and ${\displaystyle \beta _{1}}$ (e.g. 0.9) and ${\displaystyle \beta _{2}}$ (e.g. 0.999) are the forgetting factors for gradients and second moments of gradients, respectively. Squaring and square-rooting is done elementwise.

Backtracking line search is another variant of gradient descent. All of the below are sourced from the mentioned link. It is based on a condition known as the Armijo–Goldstein condition. Both methods allow learning rates to change at each iteration; however, the manner of the change is different. Backtracking line search uses function evaluations to check Armijo's condition, and in principle the loop in the algorithm for determining the learning rates can be long and unknown in advance. Adaptive SGD does not need a loop in determining learning rates. On the other hand, adaptive SGD does not guarantee the "descent property" – which Backtracking line search enjoys – which is that ${\displaystyle f(x_{n+1})\leq f(x_{n})}$ for all n. If the gradient of the cost function is globally Lipschitz continuous, with Lipschitz constant L, and learning rate is chosen of the order 1/L, then the standard version of SGD is a special case of backtracking line search.

A stochastic analogue of the standard (deterministic) Newton–Raphson algorithm (a "second-order" method) provides an asymptotically optimal or near-optimal form of iterative optimization in the setting of stochastic approximation. A method that uses direct measurements of the Hessian matrices of the summands in the empirical risk function was developed by Byrd, Hansen, Nocedal, and Singer.[24] However, directly determining the required Hessian matrices for optimization may not be possible in practice. Practical and theoretically sound methods for second-order versions of SGD that do not require direct Hessian information are given by Spall and others.[25][26][27] (A less efficient method based on finite differences, instead of simultaneous perturbations, is given by Ruppert.[28]) These methods not requiring direct Hessian information are based on either values of the summands in the above empirical risk function or values of the gradients of the summands (i.e., the SGD inputs). In particular, second-order optimality is asymptotically achievable without direct calculation of the Hessian matrices of the summands in the empirical risk function.