Fractional Underdamped Langevin Dynamics: Retargeting SGD with Momentum under Heavy-Tailed Gradient Noise

1 Introduction

Gradient-based optimization algorithms have been the de facto choice in deep learning for solving the optimization problems of the form:

x^{⋆} = a r g min x \in R^{d} {f (x) ≜ \frac{1}{n} n \sum i = 1 f^{(i)} (x)},

(1)

where $f : R^{d} \to R$

denotes the non-convex loss function,

f^{(i)}

denotes the loss contributed by an individual data point

i \in {1, \dots, n}

x \in R^{d}

denotes the collection of all the parameters of the neural network. Among others, stochastic gradient descent with momentum (SGDm) is one of the most popular algorithms for solving such optimization tasks (see e.g., Sutskever et al. (2013); Smith et al. (2018)), and is based on the following iterative scheme:

{~ v}^{k + 1} = ~ γ {~ v}^{k} - ~ η \nabla {~ f}_{k + 1} (x^{k}), x^{k + 1} = x^{k} + {~ v}^{k + 1},

(2)

where $k$ denotes the iteration number, $~ η$ is the step-size, $~ γ$ is the friction, and $~ v$ denotes the velocity (also referred to as momentum). Here, $\nabla {~ f}_{k}$ denotes the stochastic gradients defined as follows:

\nabla {~ f}_{k} (x) ≜ \frac{1}{b} \sum i \in Ω_{k} f^{(i)} (x),

(3)

where $Ω_{k} \subset {1, \dots, n}$ denotes a random subset drawn from the set of data points with $| Ω_{k} | = b ≪ n$ for all $k$ .

When the gradients are computed on all the data points (i.e., $\nabla {~ f}_{k} = \nabla f$ ), SGDm becomes deterministic and can be viewed as a discretization of the following continuous-time system (Gao et al., 2018a; Maddison et al., 2018):

d v_{t} = - (γ v_{t} + \nabla f (x_{t})) d t, d x_{t} = v_{t} d t,

(4)

where $v_{t}$ is still called the velocity. The connection between this system and (2) becomes clearer, if we discretize this system by using the Euler scheme with step-size $η$ :

		$v^{k + 1} = v^{k} - η (γ v^{k + 1} + \nabla f (x^{k})),$
		$x^{k + 1} = x^{k} + η v^{k + 1},$		(5)

and make the change of variables ${~ v}^{k} ≜ η v^{k}$ , $~ γ ≜ (1 - η γ)$ , and $~ η ≜ η^{2}$ . However, due to the presence of the stochastic gradient noise $U_{k} (x) ≜ \nabla {~ f}_{k} (x) - \nabla f (x)$ , the sequence ${x_{k}, v_{k}}_{k \in N_{+}}$ will be a stochastic process and the deterministic system (4) would not be an appropriate proxy.

Understanding the statistical properties of ${x_{k}, v_{k}}_{k \in N_{+}}$ would be of crucial importance as it might reveal the peculiar properties that lie behind the performance of SGDm for learning with neural networks. A popular approach for understanding the dynamics of stochastic optimization algorithms in deep learning is to impose some structure on the noise $U_{k}$ and relate the process (2) to a stochastic differential equation (SDE) (Mandt et al., 2016; Jastrzebski et al., 2017; Hu et al., 2017; Chaudhari and Soatto, 2018; Zhu et al., 2019; Simsekli et al., 2019)

. For instance, by assuming that the second-order moments of the stochastic gradient noise are bounded (i.e.,

E ∥ U_{k} (x) ∥^{2} < \infty

for all admissible

k

x

), one might argue that

U_{k}

can be approximated by a Gaussian random vector due to the central limit theorem (CLT)

(Fischer, 2010). Under this assumption, we might view (2) as a discretization of the following SDE, which is also known as the underdamped or kinetic Langevin dynamics:

		$d v_{t} = - (γ v_{t} + \nabla f (x_{t})) d t + \sqrt{2 γ / β} d B_{t}$
		$d x_{t} = v_{t} d t,$		(6)

where $B_{t}$ denotes the $d$ -dimensional Brownian motion and $β > 0$ is called the inverse temperature variable, measuring the noise intensity along with $γ$ . It is easy to check that, under very mild assumptions, the solution process ${x_{t}, v_{t}}_{t \geq 0}$ admits an invariant distribution whose density is proportional to $exp (- β (f (x) + ∥ v ∥^{2} / 2))$ , where the function $∥ v ∥^{2} / 2$ is often called the Gaussian kinetic energy (see e.g. (Betancourt et al., 2017)) and the distribution itself is called the Boltzmann-Gibbs measure (Pavliotis, 2014; Gao et al., 2018a; Hérau and Nier, 2004; Dalalyan and Riou-Durand, 2018). We then observe that the marginal distribution $x$ in the stationarity has a density proportional to $exp (- β f (x))$ , which indicated that any local minimum of $f$ appears as a local maximum of this density. This is a desirable property since it implies that, when the gradient noise $U_{k}$ has light tails, the process will spend more time near the local minima of $f$ . Furthermore, it has been shown that as $β$ goes to infinity, the marginal distribution of $x$ concentrates around the global optimum $x^{⋆}$ . This observation has yielded interesting results for understanding the dynamics of SGDm in the contexts of both sampling and optimization with convex and non-convex potentials $f$ (Gao et al., 2018a, b; Zou et al., 2018; Lu et al., 2016).

While the Gaussianity assumption can be accurate in certain settings such as small networks or ResNets (Martin and Mahoney, 2019; Panigrahi et al., 2019), recently it has been empirically demonstrated that in several deep learning setups, the stochastic gradient noise can exhibit a heavy-tailed behavior (Şimşekli et al., 2019; Zhang et al., 2019b). While the Gaussianity assumption would not be appropriate in this case since the conventional CLT would not hold anymore, nevertheless we can invoke the generalized CLT, which states that the asymptotic distribution of $U_{k}$ will be a symmetric $α$ -stable distribution ( $S α S$

); a class of distributions that are commonly used in the statistical physics literature as an approximation to heavy-tailed random variables

(Sliusarenko et al., 2013; Dubkov et al., 2008). As we will define in more detail in the next section, in the core of

S α S

, lies the parameter

α \in (0, 2]

, which determines the heaviness of the tail of the distribution. The tails get heavier as

α

gets smaller, the case

α = 2

reduces to the Gaussian random variables. This is illustrated in Figure 1.

Figure 1: $S α S$ densities and $L_{t}^{α}$ .

With the assumption of $U_{k}$ being $S α S$ distributed, the choice of Brownian motion will be no longer appropriate and should be replaced with an $α$ -stable Lévy motion, which motivates the following Lévy-driven SDE:


		$d x_{t} = v_{t} d t,$		(7)

where $v_{t -}$ denotes the left limit of $v_{t}$ and $L_{t}^{α}$ denotes the $α$ -stable Lévy process with independent components, which coincides with $\sqrt{2} B_{t}$ when $α = 2$ . Unfortunately, when $α < 2$ , as opposed to its Brownian counterpart, the invariant measures of such SDEs do not admit an analytical form in general; yet, one can still show that the invariant measure cannot be in the form of the Boltzmann-Gibbs measure (Eliazar and Klafter, 2003).

A more striking property of (7) was very recently revealed in a statistical physics study (Capała and Dybiec, 2019), where the authors numerically illustrated that, even when $f$ has a single minimum, the invariant measure of (7) can exhibit multiple maxima, none of which coincides with the minimum of $f$ . A similar property has been formally proven in the overdamped dynamics with Cauchy noise (i.e., $α = 1$ and $γ \to \infty$ ) by Sliusarenko et al. (2013). Since the process (7

) would spend more time around the modes of its invariant measure (i.e., the high probability region), in an optimization context (i.e., for larger

β

) the sample paths would concentrate around these modes, which might be arbitrarily distant from the optima of

f

. In other words, the heavy-tails of the gradient noise could result in an undesirable bias, which would be still present even when the step-size is taken to be arbitrarily small. As we will detail in Section 3, informally, this phenomenon stems from the fact that the heavy-tailed noise leads to aggressive updates on

v

, which are then directly transmitted to

x

due to the dynamics. Unless ‘tamed’, these updates create an hurling effect on

x

and drift it away from the modes of the “potential”

f

that is sought to be minimized.

Contributions: In this study, we develop a fractional underdamped Langevin dynamics whose invariant distribution is guaranteed to be in the form of the Boltzmann-Gibbs measure, hence its optima exactly match the optima of $f$ . We first prove a general theorem which holds for any kinetic energy function, which is not necessarily the Gaussian kinetic energy. However, it turns out that some components of the dynamics might not admit an analytical form for an arbitrary choice of the kinetic energy. Then we identify two choices of kinetic energies, where all the terms in dynamics can be written in an analytical form or accurately computable. We also analyze the Euler discretization of (14) and identify sufficient conditions for ensuring weak convergence of the ergodic averages computed over the iterates.

We observe that the discretization of the proposed dynamics has interesting algorithmic similarities with natural gradient descent (Amari, 1998) and gradient clipping (Pascanu et al., 2013), which we believe bring further theoretical understanding for their role in deep learning. Finally, we support our theory with experiments conducted on both synthetic settings and neural networks.

2 Technical Background & Related Work

The stable distributions are heavy-tailed distributions that appear as the limiting distribution of the generalized CLT for a sum of i.i.d. random variables with infinite variance

(Lévy, 1937). In this paper, we are interested in centered symmetric

α

-stable distribution. A scalar random variable

X

follows a symmetric

α

-stable distribution denoted as

X \sim S α S (σ)

if its characteristic function takes the form:

E [e^{i ω X}] = exp (- σ | ω |^{α})

ω \in R

, where

α \in (0, 2]

and

σ > 0

. Here,

α \in (0, 2]

is known as the tail-index, which determines the tail thickness of the distribution.

S α S

becomes heavier-tailed as

α

gets smaller.

σ > 0

is known as the scale parameter that measures the spread of

X

around

0

. The probability density function of a symmetric

α

-stable distribution,

α \in (0, 2]

, does not yield closed-form expression in general except for a few special cases. When

α = 1

and

α = 2

S α S

reduces to the Cauchy and the Gaussian distributions, respectively. When

0 < α < 2

α

-stable distributions have heavy-tails so that their moments are finite only up to the order

α

in the sense that

E [| X |^{p}] < \infty

if and only if

p < α

, which implies infinite variance.

Lévy motions are stochastic processes with independent and stationary increments. Their successive displacements are random and independent, and statistically identical over different time intervals of the same length, and can be viewed as the continuous-time analogue of random walks. The best known and most important examples are the Poisson process, Brownian motion, the Cauchy process and more generally stable processes. Lévy motions are prototypes of Markov processes and of semimartingales, and concern many aspects of probability theory. We refer to

(Bertoin, 1996) for a survey on the theory of Lévy motions.

In general, Lévy motions are heavy-tailed, which make it appropriate to model natural phenomena with possibly large variations, that offer occurs in statistical physics (Eliazar and Klafter, 2003), signal processing (Kuruoglu, 1999), and finance (Mandelbrot, 2013).

We define $L_{t}^{α}$ , a $d$ -dimensional symmetric $α$ -stable Lévy motion with independent components as follows. Each component of $L_{t}^{α}$ is an independent scalar $α$ -stable Lévy process, which is defined as follows: (cf. Figure 1)

$L_{0}^{α} = 0$ almost surely.
For any $t_{0} < t_{1} < \dots < t_{N}$ , the increments $L_{t_{n}}^{α} - L_{t_{n - 1}}^{α}$ are independent, $n = 1, 2, \dots, N$ .
The difference $L_{t}^{α} - L_{s}^{α}$ and $L_{t - s}^{α}$ have the same distribution: $S α S ((t - s)^{1 / α})$ for $s < t$ .
$L_{t}^{α}$ has stochastically continuous sample paths, i.e. for any $δ > 0$ and $s \geq 0$ , $P (| L_{t}^{α} - L_{s}^{α} | > δ) \to 0$ as $t \to s$ .

When $α = 2$ , we obtain a scaled Brownian motion $\sqrt{2} B_{t}$ as a special case so that the difference $L_{t}^{α} - L_{s}^{α}$ follows a Gaussian distribution $N (0, 2 (t - s))$ and $L_{t}^{α}$ is almost surely continuous. When $0 < α < 2$ , due to the stochastic continuity property, symmetric $α$ -stable Lévy motions can have have a countable number of discontinuities, which are often known as jumps. The sample paths are continuous from the right and they have left limits, a property known as càdlàg (Duan, 2015).

Recently, Şimşekli (2017) extended the overdamped Langevin dynamics to an SDE driven by $L_{t}^{α}$ , given as:¹¹1In Şimşekli (2017), (8) does not contain an inverse temperature $β$ , which was later on introduced in Nguyen et al. (2019).

d x_{t} = b (x_{t -}, α) d t + β^{- 1 / α} d L_{t}^{α},

(8)

where the drift $b (x, α) = ((b (x, α))_{i}, 1 \leq i \leq d)$ is defined as follows:

(b (x, α))_{i} = - D_{x_{i}}^{α - 2} (ϕ (x) \partial_{x_{i}} f (x)) / ϕ (x) .

(9)

Here, $ϕ (x) = exp (- f (x))$ and $D$ denotes the fractional Riesz derivative (Riesz, 1949):

D^{γ} u (x) := F^{- 1} {| ω |^{γ} (F (u)) (ω)} (x),

(10)

where $F$

denotes the Fourier transform. Briefly,

D^{γ}

extends usual differentiation to fractional orders and when

γ = 2

it coincides (up to a sign difference) with the usual second-order derivative

- d^{2} f (x) / d x^{2}

The important property of the process (8) is that it admits an invariant distribution whose density is proportional to $exp (- β f (x))$ (Nguyen et al., 2019). It is easy to show that, when $α = 2$ , the drift reduces to $b (x, 2) = - \nabla f (x)$ , hence we recover the classical overdamped dynamics:

d x_{t} = - \nabla f (x_{t}) d t + \sqrt{2 / β} d B_{t} .

(11)

Since the fractional Riesz derivative is costly to compute, Şimşekli (2017) proposed an approximation of $b (x, α)$ based on the alternative definition of $D$ given in (Ortigueira, 2006), such that:

b (x, α) \approx - c_{α} \nabla f (x),

(12)

where $c_{α} := Γ (α - 1) / Γ (α / 2)^{2}$ . This approximation essentially results in replacing $B_{t}$ with $L_{t}^{α}$ in (11) in a rather straightforward manner. While avoiding the computational issues originated from the Riesz derivatives, as shown in (Nguyen et al., 2019), this approximation can induce an arbitrary bias in a non-convex optimization context. Besides, the stationary distribution of this approximated dynamics was analytically derived in (Sliusarenko et al., 2013) under the choice of $α = 1$ and $f (x) = x^{4} / 4 - a x^{2} / 2$ for $x \in R^{1}$ and $a > 0$ . These results show that, in the presence of heavy-tailed perturbations, the drift should be modified, otherwise an inaccurate approximation of the Riesz derivatives can result in an explicit bias, which moves the modes of the distribution away from the modes of $f$ .

From a pure Monte Carlo perspective, Ye and Zhu (2018) extended the fractional overdamped dynamics (8) to higher-order dynamics and proposed the so-called fractional Hamiltonian dynamics (FHD), given as follows:

$d x_{t} =$	$D^{α - 2} {ϕ (z_{t}) v_{t}} / ϕ (z_{t}) d t,$
$d v_{t} =$	$- D^{α - 2} {ϕ (z_{t}) \nabla f (x_{t})} / ϕ (z_{t}) d t$
	$- γ D^{α - 2} {ϕ (z_{t}) v_{t}} / ϕ (z_{t}) d t + γ^{1 / α} d L_{t}^{α},$	(13)

where $z_{t} = (x_{t}, v_{t})$ , and $ϕ (z) = e^{- f (x) - \frac{1}{2} ∥ v ∥^{2}}$ . They showed that the invariant measure of the process has a density proportional to $ϕ (z)$ , i.e., the Boltzmann-Gibbs measure. Similar to the overdamped case (8), the Riesz derivatives do not admit an analytical form in general. Hence they approximated them by using the same approximation given in (12), which yields the SDE given in (7) (up to a scaling factor). This observation also confirms that the heavy-tailed noise requires an adjustment in the dynamics, otherwise the induced bias might drive the dynamics away from the minima of $f$ (Capała and Dybiec, 2019).

3 Fractional Underdamped Langevin Dynamics

In this section, we develop the fractional underdamped Langevin dynamics (FULD), which is expressed by the following SDE:

		$d v_{t} = - (γ c (v_{t -}, α) + \nabla f (x_{t})) d t + (γ / β)^{1 / α} d L_{t}^{α},$
		$d x_{t} = \nabla g (v_{t}) d t,$		(14)

where $c : R^{d} \times (0, 2] \mapsto R^{d}$ is the drift function for the velocity and $g : R^{d} \mapsto R$ denotes a general notion of kinetic energy. In the next theorem, which is the main theoretical result of this paper, we will identify the relation between these two functions such that the solution process will keep the generalized Boltzmann-Gibbs measure, invariant. All the proofs are given in the appendix. Let $c (v, α) = ((c (v, α))_{i}, 1 \leq i \leq d)$ has the following form:

(c (v, α))_{i} := \frac{D_{v_{i}}^{α - 2} (ψ (v) \partial_{v_{i}} g (v))}{ψ (v)}, ψ (v) := e^{- g (v)} .

(15)

The the measure $π (d x, d v) \propto e^{- β (f (x) + g (v))} d x d v$ on $R^{d} \times R^{d}$ is an invariant probability measure for the Markov process $(x_{t}, v_{t})$ . One of the main features of FULD is that the fractional Riesz derivatives only appears in the drift $c$ , which only depends on $v$ . This is highly in contrast with FHD (13), where the Riesz derivatives are taken over both $x$ and $v$ , which is the source of intractablility. Moreover, FULD enjoys the freedom to choose different kinetic energy functions $g (v)$ . In the sequel, we will investigate two options for $g$ , such that the drift $c$ can be analytically obtained.

3.1 Gaussian kinetic energy

In classical overdamped Langevin dynamics and Hamiltonian dynamics, the default choice of kinetic energy is the Gaussian kinetic energy, which corresponds to taking $g (v) = \frac{1}{2} ∥ v ∥^{2}$ (Neal, 2010; Livingstone et al., 2019; Dalalyan and Riou-Durand, 2018). With this choice, the fractional dynamics becomes:


	$d x_{t} = v_{t} d t .$		(16)

In the next result, we will show that in this case, the drift $c$ admits an analytical solution. Let $g (v) = \frac{1}{2} ∥ v ∥^{2}$ . Then, for any $1 \leq i \leq d$ ,

(c (v, α))_{i} = \frac{2^{\frac{α}{2}} v_{i}}{\sqrt{π}} Γ (\frac{α + 1}{2}) {_{1} F}_{1} (\frac{2 - α}{2}; \frac{3}{2}; \frac{v_{i}^{2}}{2}),

(17)

where $Γ$ is the gamma function and $_{1} F_{1}$ is the Kummer confluent hypergeometric function. In particular, when $α = 2$ , we have $(c (v, α))_{i} = v_{i}$ . We observe that the fractional dynamics (16) strictly extends the underdamped Langevin dynamics (6) as $c (v, 2) = v$ .

Let us now investigate the form of the new drift $c$ and its implications. In Figure 2, we illustrate $c$ for the $d = 1$ dimensional case (note that for $d > 1$ , each component of $c$ still behaves like Figure 2). We observe that due to the hypergeometric function $_{1} F_{1}$ , the drift grows exponentially fast with $| v |$ whenever $α < 2$ . Semantically, this means that, in order to be able to compensate the large jumps incurred by $L_{t}^{α}$ , the drift has to react very strongly and hence prevent $v$ to take large values. To illustrate this behavior, we provide more visual illustrations in the appendix.

Figure 2: Illustration of one dimensional a) drift function $c$ for the Gaussian kinetic energy, b) $\nabla G_{α}$ for the $S α S$ kinetic energy.

Even though this aggressive behavior of $c$ can be beneficial for the continuous-time system, it is unfortunately clear that its Euler-Maruyama discretization will not yield a practical algorithm due to the same behavior. Indeed, we would need the function $c$ to be Lipschitz continuous in order to guarantee the algorithmic stability of its discretization (Kloeden and Platen, 2013); however, if we consider the integral form of $_{1} F_{1}$ (cf. (Abramowitz and Stegun, 1972)), we observe that the function

(c (v, α))_{i} = \frac{2^{\frac{α}{2}} v_{i}}{\sqrt{π}} \cdot \frac{Γ (\frac{3}{2})}{Γ (\frac{2 - α}{2})} \int_{0}^{1} e^{\frac{v_{i}^{2}}{2} t} t^{- \frac{α}{2}} (1 - t)^{\frac{α - 1}{2}} d t

is clearly not Lipschitz continuous in $v_{i}$ . Therefore, we conclude that FULD with the Gaussian kinetic energy is mostly of theoretical interest.

3.2 Alpha-stable kinetic energy

The dynamics with the Gaussian kinetic energy requires a very strong drift $c$ mainly because we force the dynamics to make sure that the invariant distribution of $v$ to be a Gaussian. Since the Gaussian distribution has light-tails, it cannot tolerate samples with large magnitudes, hence requires a large dissipation to make sure $v$ does not take large values.

In order to avoid such an explosive drift that potentially degrades practicality, next we explore heavy-tailed kinetic energies, which would allow the components of $v$ to take large values, while still making sure that the drift $c$ in (15) admits an analytical form.

In our next result, we show that, when we choose an $S α S$ kinetic energy, such that the tail-index $α$ of this kinetic energy matches the one of the driving process $L_{t}^{α}$ , the drift $c$ simplifies and becomes the identity function. Let $e^{- g_{α} (v)}$ be the probability density function of $S α S (\frac{1}{α})$ . Choose $ψ (v) = e^{- G_{α} (v)}$ in (15), where $G_{α} (v) = \sum_{i = 1}^{d} g_{α} (v_{i})$ for any $v = (v_{1}, \dots, v_{d})$ . Then,

(c (v, α))_{i} = v_{i}, 1 \leq i \leq d .

(18)

This result hints that, perhaps $G_{α} (v)$ is the natural choice of kinetic energy for the systems driven by $L_{t}^{α}$ .

It now follows from Theorem 3.2 that the FULD with $α$ -stable kinetic energy reduces to the following SDE:

	$d v_{t} = - γ v_{t -} d t - \nabla f (x_{t}) d t + γ^{1 / α} d L_{t}^{α},$
			(19)

It can be easily verified that $\nabla G_{α} (v_{t}) = v_{t}$ for $α = 2$ , as $g_{2} (v) = \frac{1}{2} log 2 π + \frac{1}{2} v^{2}$ , hence, the SDE (19) also reduces to the classical underdamped Langevin dynamics (6).

While this choice of $g$ results in an analytically available $c$ , unfortunately the function $\nabla G_{α}$ itself admits a closed-form analytical formula only when $α = 1$ or $α = 2$ , due to the properties of the $S α S$ densities. Nevertheless, as $\nabla G_{α}$ is based on one-dimensional $S α S$ densities, it can be very accurately computed by using the recent methods developed in (Ament and O’Neil, 2018).

We visually inspect the behavior of $\nabla G_{α}$ in Figure 2 for dimension one. We observe that, as soon as $α < 2$ , $\nabla G_{α}$ takes a very smooth form. Besides, for small $| v |$ the function behaves like a linear function and when $| v |$ goes to infinity, it vanishes. This behavior can be interpreted as follows: since $v$ can take larger values due to the heavy tails of the kinetic energy, in order to be able target the correct distribution, the dynamics compensates the potential bursts in $v$ by passing it through the asymptotically vanishing $\nabla G_{α}$ .

3.3 Euler discretization and weak convergence analysis

As visually hinted in Figure 2, the function $\nabla G_{α}$ has strong regularity, which makes (19) to be potentially beneficial for practical implementations. Indeed, it is easy to verify that $\nabla G_{α}$ is Lipschitz continuous for $α = 1$ and $2$ , and in our next result, we show that this observation is true for any admissible $α$ , which is a desired property when discretizing continuous-time dynamics. For $0 < α \leq 2$ , the map $v \mapsto g_{α} (v)$ is Lipschitz continuous, hence $v \mapsto \nabla G_{α} (v)$ is also Lipschitz continuous. Accordingly we consider the following Euler-Maruyama discretization for (19):

	$v^{k + 1}$	$= {~ γ}_{k} v^{k} - η_{k} \nabla f (x^{k}) + (η_{k} γ / β)^{1 / α} s^{k + 1},$
	$x^{k + 1}$	$= x^{k} + η_{k} \nabla G_{α} (v^{k + 1}),$		(20)

where ${~ γ}_{k} = 1 - γ η_{k}$ , $s^{k}$ is a random vector whose components are independently $S α S (1)$ distributed, and $(η_{k})_{k \in N_{+}}$ is a sequence of step-sizes.

In this section, we analyze the weak convergence of the ergodic averages computed by using (20). Given a test function $h$ , consider its expectation with respect to the target measure $π$ , i.e. $π (h) := E_{X \sim π} [h (X)] = \int h (x) π (d x)$ with $π (d x) \propto exp (- β f (x)) d x$ . We will discuss next how this expectation can be approximated through the sample averages

{¯ π}_{K} (h) := (1 / S_{K}) K \sum k = 1 η_{k} h (x^{k}),

(21)

where $S_{K} := \sum_{k = 1}^{K} η_{k}$ is the cumulative sum of the step-size sequence. We note that this notion of convergence is stronger than the convergence of (20) near a local minimum since it requires the convergence of the measure itself, and our analysis can be extended to a global optimization context by using the techniques presented in (Raginsky et al., 2017).

We now present the assumptions that imply our results.

A 1

The step-size sequence ${η_{k}}$ is non-increasing and satisfies ${lim}_{k \to \infty} η_{k} = 0$ and ${lim}_{K \to \infty} S_{K} = \infty$ .

A 2

Let $V : R^{2 d} \to R_{+}$ be a twice continuously differentiable function, satisfying ${lim}_{∥ z ∥ \to \infty} V (z) = \infty$ , $∥ \nabla V ∥ \leq C \sqrt{V}$ for some $C > 0$ and has a bounded Hessian $\nabla^{2} V$ . Given $p \in (0, \frac{1}{2}]$ , there exists $a \in (1 - \frac{p}{2}, 1]$ , $β_{1} \in R$ , $β_{2} > 0$ such that $∥ b ∥^{2} \leq C V^{a}$ and $⟨ \nabla V, b ⟩ \leq β_{1} - β_{2} V^{a}$ where $b (v, x) = (- γ v - \nabla f (x), \nabla G_{α} (v))$ is the drift of the $(v_{t}, x_{t})$ process defined in (19).

These are common assumptions ensuring that the SDE is simulated with infinite time-horizon and the process is not explosive (Panloup, 2008; Şimşekli, 2017). We can now establish the weak convergence of (21) and present it as a corollary to Theorem 3, Proposition 3.3, and (Panloup, 2008) (Theorem 2). Assume that the gradient $\nabla f$ is Lipschitz continuous and has linear growth i.e., there exits $C > 0$ such that $∥ \nabla f (x) ∥ \leq C (1 + ∥ x ∥)$ for all $x$ . Furthermore, assume that Assumptions 1 and 2 hold for some $p \in (0, 1 / 2]$ . If the test function $h = o (V^{\frac{p}{2} + a - 1})$ then

{¯ π}_{K} (h) \to π (h) almost surely as K \to \infty .

3.4 Connections to existing approaches

We now point out interesting algorithmic connections between (20) and two methods that are commonly used in practice. We first roll back our initial hypothesis that the gradient noise is $S α S$ distributed, i.e., $\nabla {~ f}_{k} (x) = \nabla f (x) + (η_{k} γ / β)^{1 / α} s^{k}$ , and modify (20) as follows:

	$v^{k + 1}$	$= {~ γ}_{k} v^{k} - η_{k} \nabla {~ f}_{k + 1} (x^{k}),$
	$x^{k + 1}$	$= x^{k} + η_{k} \nabla G_{α} (v^{k + 1}) .$		(22)

As a special case when ${~ γ}_{k} = 0$ , we obtain a stochastic gradient descent-type recursion:

x^{k + 1} = x^{k} - η_{k}^{2} \nabla G_{α} (\nabla {~ f}_{k + 1} (x^{k})) .

(23)

Let us now consider gradient-clipping

, a heuristic approach for eliminating the problem of ‘exploding gradients’, which often appear in training neural networks

(Pascanu et al., 2013; Zhang et al., 2019a). Very recently, Zhang et al. (2019b) empirically illustrated that such explosions stem from heavy-tailed gradients and formally proved that gradient clipping indeed improves convergence rates under heavy-tailed perturbations. We notice that, the behavior of (22) is reminiscent of gradient clipping: due to the vanishing behavior of

\nabla G_{α}

for

α < 2

, as the components of

v^{k}

gets larger in magnitude, the update applied on

x^{k}

gets smaller. The behavior becomes more prominent in (23). On the other hand, (22) is more aggressive in the sense that the updates can get arbitrarily small as the value of

α

decreases as opposed to being ‘clipped’ with a threshold.

The second connection is with the natural gradient descent algorithm, where the stochastic gradients are pre-conditioned with the inverse Fisher information matrix (FIM) (Amari, 1998). Here FIM is defined as $E [\nabla f (x) \nabla f (x)^{⊤}]$ , where the expectation is taken over the data. Notice that when $α = 1$

(i.e., Cauchy distribution), we have the following form:

\nabla G_{1} (v) = (\frac{2 v_{1}}{v_{1}^{2} + 1}, \dots, \frac{2 v_{d}}{} v_{d}^{2} + 1)

. Therefore, we observe that, in (23),

\nabla G_{1} (\nabla {~ f}_{k} (x))

can be equivalently written as

M_{k} (x)^{- 1} \nabla {~ f}_{k} (x)

, where

M_{k} (x)

is a diagonal matrix with entries

m_{i i} = ((\nabla {~ f}_{k} (x))_{i}^{2} + 1) / 2

. Therefore, we can see

M_{k}

as an estimator of the diagonal part of FIM, as they will be in the same order when

| (\nabla {~ f}_{k} (x))_{i} |

is large. Besides, (22) then appears as its momentum extension. However,

M_{k}

will be biased mainly due to the fact that FIM is the average of the squared gradients, whereas

M_{k}

is based on the square of the average gradients. This connection is rather surprising, since a seemingly unrelated, differential geometric approach turns out to have strong algorithmic similarities with a method that naturally arises when the gradient noise is Cauchy distributed.

Figure 3: Estimated invariant measures for the quartic potential: top $α = 1$ , bottom $α = 1.9$ .

4 Numerical Study

In this section, we will illustrate our theory on several experiments which are conducted in both synthetic and real-data settings. We note that, as expected, FULD with Gaussian kinetic energy did not yield a numerically stable discretization due to the explosive behavior of $c$ . Hence, in this section, we only focus on FULD with $S α S$ kinetic energy in order to avoid obscuring the results and from now on we will simply refer to FULD with $S α S$ kinetic energy as FULD.

4.1 Synthetic setting

We first consider a one-dimensional synthetic setting, similar to the one considered in (Capała and Dybiec, 2019). We consider a quartic potential function with a quadratic component, $f (x) = x^{4} / 4 - x^{2} / 2$ . We then simulate the ‘uncorrected dynamics’ (UD) given in (7) and FULD (19) by using the Euler-Maruyama discretization to compare their behavior for different $α$ . For $α \notin {1, 2}$ , we used the software given in (Ament and O’Neil, 2018) for computing $\nabla G_{α}$ .

Figure 3 illustrates the distribution of the samples generated by simulating the two dynamics. In this setup, we set $β = 1$ , $η = 0.01$ , $γ = 10$ with number of iterations $K = 50000$ . We observe that, for $α = 1.9$ , FULD very accurately captures the form of the distribution, whereas UD exhibits a visible bias and the shape of its resulting distribution is slightly distorted. Nevertheless, since the perturbations are close to a Gaussian in this case (i.e., $α$ is close to $2$ ), the difference is not substantial and can be tolerable in an optimization context. However, this behavior becomes much more emphasized when we use a heavier-tailed driving process: when $α = 1$ , we observe that the target distribution of UD becomes distant from the Gibbs measure $exp (- f (x))$ , and more importantly its modes no longer match the minima of $f$ ; agreeing with the observations presented in (Capała and Dybiec, 2019). On the other hand, thanks to the correction brought by $\nabla G_{α}$ , FULD still captures the target distribution very accurately, even when the driving force is Cauchy.

On the other hand, in our experiments we observed that, for small values of $α$ , UD can quickly become numerically unstable and even diverge for slightly larger step-sizes, whereas this problem never occurred for FULD. This outcome also stems from the fact that UD does not have any mechanism to compensate the potential large updates originating from the heavy-tailed perturbations. To illustrate this observation more clearly, in Figure 4 we illustrate the iterates $(x^{k})_{k = 1}^{K}$ which were used for producing Figure 3. We observe that, while the iterates of UD are well-behaved for $α = 1.9$ , the magnitude range of the iterates gets quite large when $α$ is set to $1$ . On the other hand, for both values of $α$ , FULD iterates are always kept in a reasonable range, thanks to the clipping-like effect of $\nabla G_{α}$ .

Figure 4: Illustration of the iterates for the quartic potential: top $α = 1$ , bottom $α = 1.9$ .

4.2 Neural networks

In our next set of experiments, we evaluate our theory on neural networks. In particular, we apply the iterative scheme given in (22) as an optimization algorithm for training neural networks, and compare its behavior with classical SGDm defined in (2). In this setting, we do not add any explicit noise, all the stochasticity comes from the potentially heavy-tailed stochastic gradients (3).

Figure 5: Neural network results on MNIST.

We consider a fully-connected network for a classification task on the MNIST and CIFAR10 datasets, with different depths (i.e. number of layers) and widths (i.e. number of neurons per layer). For each depth-width pair, we train two neural networks by using SGDm (

2) and our modified version (22), and compare their final train/test accuracies and loss values. We use the conventional train-test split of the datasets: for MNIST we have

60

K training and

10

K test samples, and for CIFAR10 these numbers are

50

K and

10

K, respectively. We use the cross entropy loss (also referred to as the ‘negative-log-likelihood’).

We note that the modified scheme (22) reduces to (2) when $α = 2$ , since $\nabla G_{2} (v) = v$ . Hence in this section, we will refer to SGDm as the special case of (22) with $α = 2$ . On the other hand, in these experiments, computing $\nabla G_{α}$ becomes impractical for $α \notin {1, 2}$ , since the algorithms given in (Ament and O’Neil, 2018) become prohibitively slow with the increased dimension $d$ . Hence, we will focus on the analytically available Cauchy case, i.e., $α = 1$ , which can be efficiently implemented (cf. its definition in Section 3.4). We expect that, if the stochastic gradient noise can be well-approximated by using a Cauchy distribution, then the modified dynamics should exhibit an improved performance since it would eliminate the potential bias brought by the heavy-tailed noise.

In these experiments, we set $η = 0.1$ , $γ = 0.1$ , and run the algorithms for $K = 10000$ iterations ²²2Since the scale of the gradient noise is proportional to $(γ / β)^{\frac{1}{α}}$ (see (20)), in this setup, a fixed $γ$ implicitly determines $β$ . . We measure the accuracy and the loss at every 100th iteration and we report the average of the last two measurements. Figure 5 shows the results obtained on the MNIST dataset. We observe that, in most of the cases, setting $α = 1$ yields a better performance in terms both training and testing accuracies/losses. This difference becomes more visible when the width is set to $256$ : the accuracy difference between the algorithms reaches to $0.5 %$ . We obtain a similar result on the CIFAR10 dataset, as illustrated in Figure 6. In most of the cases $α = 1$ performs better, with the maximum accuracy difference being $\approx 1.0 %$ , implying the gradient noise can be approximated by a Cauchy random variable.

We observed a similar behavior when the width was set to $64$ . However, when we set the width to $32$ we did not perceive a significant difference in terms of the performance of the algorithms. We suspect that this behavior is due to either the gradient noise is not well-approximated either by a Gaussian or Cauchy. On the other hand, when the width was set to $512$ , $α = 2$ resulted in a slightly better performance, which would be an indication that the Gaussian approximation is closer. The corresponding figures are provided in the appendix.

Figure 6: Neural network results on CIFAR10.

5 Conclusion and Future Directions

We considered the continuous-time variant of SGDm, known as the underdamped Langevin dynamics (ULD), and developed theory for the case where the gradient noise can be well-approximated by a heavy-tailed $α$ -stable random vector. As opposed to naïvely replacing the driving stochastic force in ULD, which correspondonds to running SGDm with heavy-tailed gradient noise, the dynamics that we developed exactly target the Boltzmann-Gibbs distribution, and hence do not introduce an implicit bias. We further established the weak convergence of the Euler-Maruyama discretization and illustrated interesting connections between the discretized algorithm and existing approaches commonly used in practice. We supported our theory with experiments on a synthetic setting and fully connected neural networks.

Our framework opens up interesting future directions. Our current modeling strategy requires a state-independent, isotropic noise assumption, which would not accurately reflect the reality. While anisotropic noise can be incorporated to our framework by using the approach of Ye and Zhu (2018), state-dependent noise introduces challenging technical difficulties. Similarly, it has been illustrated that the tail-index $α$ can depend on the state and different components of the noise can have a different $α$ (Şimşekli et al., 2019). Incorporating such state dependencies would be an important direction of future research.

Acknowledgments

We thank Jingzhao Zhang for fruitful discussions. The contribution of Umut Şimşekli to this work is partly supported by the French National Research Agency (ANR) as a part of the FBIMATRIX (ANR-16-CE23-0014) project, and by the industrial chair Data science & Artificial Intelligence from Télécom Paris. Lingjiong Zhu is grateful to the support from Simons Foundation Collaboration Grant. Mert Gürbüzbalaban acknowledges support from the grants NSF DMS-1723085 and NSF CCF-1814888.

References

Abramowitz and Stegun (1972) M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables. Dover, New York, 1972.
Amari (1998) Shun-Ichi Amari. Natural gradient works efficiently in learning. Neural Computation, 10(2):251–276, 1998.
Ament and O’Neil (2018) Sebastian Ament and Michael O’Neil. Accurate and efficient numerical calculation of stable densities via optimized quadrature and asymptotics. Statistics and Computing, 28:171–185, 2018.
Bertoin (1996) Jean Bertoin. Lévy Processes. Cambridge University Press, 1996.
Betancourt et al. (2017) Michael Betancourt, Simon Byrne, Sam Livingstone, Mark Girolami, et al. The geometric foundations of Hamiltonian Monte Carlo. Bernoulli, 23(4A):2257–2298, 2017.
Brosse et al. (2019) Nicolas Brosse, Alain Durmus, Éric Moulines, and Sotirios Sabanis. The tamed unadjusted Langevin algorithm. Stochastic Processes and their Applications, 129(10):3638–3663, 2019.
Capała and Dybiec (2019) Karol Capała and Bartłomiej Dybiec. Stationary states for underdamped anharmonic oscillators driven by Cauchy noise. arXiv preprint arXiv:1905.12078, 2019.
Chaudhari and Soatto (2018) P. Chaudhari and S. Soatto. Stochastic gradient descent performs variational inference, converges to limit cycles for deep networks. In International Conference on Learning Representations, 2018.
Şimşekli (2017) Umut Şimşekli.
Fractional Langevin Monte Carlo: Exploring Lévy driven stochastic differential equations for Markov Chain Monte Carlo.
In International Conference on Machine Learning
, pages 3200–3209, 2017.
Dalalyan and Riou-Durand (2018) Arnak S Dalalyan and Lionel Riou-Durand. On sampling from a log-concave density using kinetic Langevin diffusions. arXiv preprint arXiv:1807.09382, 2018.
Duan (2015) J. Duan. An Introduction to Stochastic Dynamics. Cambridge University Press, New York, 2015.
Dubkov et al. (2008) Alexander A Dubkov, Bernardo Spagnolo, and Vladimir V Uchaikin. Lévy flight superdiffusion: An introduction. International Journal of Bifurcation and Chaos, 18(09):2649–2672, 2008.
Eliazar and Klafter (2003) Iddo Eliazar and Joseph Klafter. Lévy-driven Langevin systems: Targeted stochasticity. Journal of Statistical Physics, 111(3-4):739–768, 2003.
Fischer (2010) Hans Fischer. A History of the Central Limit Theorem: From Classical to Modern Probability Theory. Springer Science & Business Media, 2010.
Gao et al. (2018a) Xuefeng Gao, Mert Gürbüzbalaban, and Lingjiong Zhu. Global convergence of Stochastic Gradient Hamiltonian Monte Carlo for non-convex stochastic optimization: Non-asymptotic performance bounds and momentum-based acceleration. arXiv:1809.04618, 2018a.
Gao et al. (2018b) Xuefeng Gao, Mert Gürbüzbalaban, and Lingjiong Zhu. Breaking reversibility accelerates Langevin dynamics for global non-convex optimization. arXiv:1812.07725, 2018b.
Hérau and Nier (2004) Frédéric Hérau and Francis Nier. Isotropic hypoellipticity and trend to equilibrium for the Fokker-Planck equation with a high-degree potential. Archive for Rational Mechanics and Analysis, 171(2):151–218, 2004.
Hu et al. (2017) W. Hu, C. J. Li, L. Li, and J.-G. Liu. On the diffusion approximation of nonconvex stochastic gradient descent. arXiv preprint arXiv:1705.07562, 2017.
Jastrzebski et al. (2017) S. Jastrzebski, Z. Kenton, D. Arpit, N. Ballas, A. Fischer, Y. Bengio, and A. Storkey. Three factors influencing minima in SGD. arXiv preprint arXiv:1711.04623, 2017.
Kloeden and Platen (2013) Peter E Kloeden and Eckhard Platen. Numerical Solution of Stochastic Differential Equations, volume 23. Springer Science & Business Media, 2013.
Kuruoglu (1999) E. E. Kuruoglu. Signal Processing in $α$ -Stable Noise Environments: A Least $ℓ_{p}$ -Norm Approach. PhD Thesis, University of Cambridge, 1999.
Lévy (1937) P. Lévy. Théorie de l’addition des variables aléatoires. Gauthiers-Villars, Paris, 1937.
Livingstone et al. (2019) Samuel Livingstone, Michael F Faulkner, and Gareth O Roberts. Kinetic energy choice in Hamiltonian/hybrid Monte Carlo. Biometrika, 106(2):303–319, 2019.
Lu et al. (2016) Xiaoyu Lu, Valerio Perrone, Leonard Hasenclever, Yee Whye Teh, and Sebastian J Vollmer. Relativistic Monte Carlo. arXiv preprint arXiv:1609.04388, 2016.
Maddison et al. (2018) Chris J Maddison, Daniel Paulin, Yee Whye Teh, Brendan O’Donoghue, and Arnaud Doucet. Hamiltonian descent methods. arXiv preprint arXiv:1809.05042, 2018.
Mandelbrot (2013) B. B. Mandelbrot. Fractals and Scaling in Finance:Discontinuity, Concentration, Risk. Springer Science & Business Media, 2013.
Mandt et al. (2016) S. Mandt, M. Hoffman, and D. Blei. A variational analysis of stochastic gradient algorithms. In International Conference on Machine Learning, pages 354–363, 2016.
Martin and Mahoney (2019) Charles H Martin and Michael W Mahoney. Heavy-tailed universality predicts trends in test accuracies for very large pre-trained deep neural networks. arXiv preprint arXiv:1901.08278, 2019.
Montroll and Bendler (1984) Elliott W. Montroll and John T. Bendler. On Lévy (or stable) distributions and the Williams-Watts model of dielectric relaxation. Journal of Statistical Physics, 34(1):129–162, Jan 1984. ISSN 1572-9613.
Montroll and West (1979) Elliott W. Montroll and Bruce J. West. Chapter 2 - On an enriched collection of stochastic processes. In E.W. Montroll and J.L. Lebowitz, editors, Fluctuation Phenomena, pages 61 – 175. Elsevier, 1979. ISBN 978-0-444-85248-9.
Neal (2010) RM Neal. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo (S. Brooks, A. Gelman, G. Jones, and X.-L. Meng, eds.), 2010.
Nguyen et al. (2019) Thanh Huy Nguyen, Umut Şimşekli, and Gaël Richard.
Non-Asymptotic Analysis of Fractional Langevin Monte Carlo for Non-Convex Optimization.
In International Conference on Machine Learning, pages 4810–4819, 2019.
Ortigueira (2006) M. D. Ortigueira. Riesz potential operators and inverses via fractional centred derivatives. International Journal of Mathematics and Mathematical Sciences, 2006, 2006.
Panigrahi et al. (2019) Abhishek Panigrahi, Raghav Somani, Navin Goyal, and Praneeth Netrapalli. Non-Gaussianity of stochastic gradient noise. arXiv preprint arXiv:1910.09626, 2019.
Panloup (2008) Fabien Panloup. Recursive computation of the invariant measure of a stochastic differential equation driven by a Lévy process. Annals of Applied Probability, 18(2):379–426, 2008.
Pascanu et al. (2013) Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio.
On the difficulty of training recurrent neural networks.
In International conference on machine learning, pages 1310–1318, 2013.
Pavliotis (2014) Grigorios A Pavliotis. Stochastic Processes and Applications: Diffusion processes, the Fokker-Planck and Langevin Equations, volume 60. Springer, 2014.
Raginsky et al. (2017) Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky. Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis. In Conference on Learning Theory, pages 1674–1703, 2017.
Riesz (1949) M. Riesz. L’intégrale de Riemann-Liouville et le problème de Cauchy. Acta Mathematica, 81(1):1–222, 1949.
Schertzer et al. (2001) D. Schertzer, M. Larchevêque, J. Duan, V.V. Yanovsky, and S. Lovejoy. Fractional Fokker-Planck equation for nonlinear stochastic differential equations driven by non-Gaussian Lévy stable noises. Journal of Mathematical Physics, 42(1):200–212, 2001.
Şimşekli et al. (2019) Umut Şimşekli, Mert Gürbüzbalaban, Thanh Huy Nguyen, Gaël Richard, and Levent Sagun. On the heavy-tailed theory of stochastic gradient descent for deep neural networks. arXiv preprint arXiv:1912.00018, 2019.
Simsekli et al. (2019) Umut Simsekli, Levent Sagun, and Mert Gurbuzbalaban. A tail-index analysis of stochastic gradient noise in deep neural networks. In International Conference on Machine Learning, pages 5827–5837, 2019.
Sliusarenko et al. (2013) O Yu Sliusarenko, DA Surkov, V Yu Gonchar, and Aleksei V Chechkin. Stationary states in bistable system driven by Lévy noise. The European Physical Journal Special Topics, 216(1):133–138, 2013.
Smith et al. (2018) Samuel L. Smith, Pieter-Jan Kindermans, Chris Ying, and Quoc V. Le. Don’t decay the learning rate, increase the batch size. In International Conference on Learning Representations, 2018.
Sutskever et al. (2013) Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton. On the importance of initialization and momentum in deep learning. In International Conference on Machine Learning, pages 1139–1147, 2013.
Wintner (1941) Aurel Wintner. The singularities of Cauchy’s distributions. Duke Math. J., 8(4):678–681, 12 1941. doi: 10.1215/S0012-7094-41-00857-8.
Ye and Zhu (2018) Nanyang Ye and Zhanxing Zhu. Stochastic fractional Hamiltonian Monte Carlo. In Proceedings of the 27th International Joint Conference on Artificial Intelligence, 2018.
Zhang et al. (2019a) Jingzhao Zhang, Tianxing He, Suvrit Sra, and Ali Jadbabaie. Analysis of gradient clipping and adaptive scaling with a relaxed smoothness condition. arXiv preprint arXiv:1905.11881, 2019a.
Zhang et al. (2019b) Jingzhao Zhang, Sai Praneeth Karimireddy, Andreas Veit, Seungyeon Kim, Sashank J Reddi, Sanjiv Kumar, and Suvrit Sra. Why ADAM beats SGD for attention models. arXiv preprint arXiv:1912.03194, 2019b.
Zhu et al. (2019) Zhanxing Zhu, Jingfeng Wu, Bing Yu, Lei Wu, and Jinwen Ma. The anisotropic noise in stochastic gradient descent: Its behavior of escaping from sharp minima and regularization effects. In Proc. Int. Conf. Mach. Learn., pages 7654–7663, 2019.
Zou et al. (2018) Difan Zou, Pan Xu, and Quanquan Gu. Stochastic variance-reduced Hamilton Monte Carlo methods. In International Conference on Machine Learning, pages 6028–6037, 2018.

Appendix A Proof of Theorem 3

Let $q (x, v, t)$ denote the probability density of $(x_{t}, v_{t})$ . Then it satisfies the fractional Fokker-Planck equation (see Proposition 1 and Section 7 in (Schertzer et al., 2001)):

	$\partial_{t} q (x, v, t)$	$= γ d \sum i = 1 \frac{\partial [(c (v, α))_{i} q (x, v, t)]}{\partial v_{i}} + d \sum i = 1 \frac{\partial [\partial_{x_{i}} f (x) q (x, v, t)]}{\partial v_{i}} - \frac{γ}{β} d \sum i = 1 D_{v_{i}}^{α} q (x, v, t)$
		$- d \sum i = 1 \frac{\partial [(\partial_{v_{i}} g (v) q (x, v, t)]}{\partial x_{i}}$

We can compute that

	$γ d \sum i = 1 \frac{\partial [(c (v, α))_{i} ϕ (x) ψ (v)]}{\partial v_{i}} + d \sum i = 1 \frac{\partial [\partial_{x_{i}} f (x) ϕ (x) ψ (v)]}{\partial v_{i}} - \frac{γ}{β} d \sum i = 1 D_{v_{i}}^{α} ϕ (x) ψ (v) - d \sum i = 1 \frac{\partial [(\partial_{v_{i}} g (v) ϕ (x) ψ (v)]}{\partial x_{i}}$
	$= \frac{γ}{β} ϕ (x) [β d \sum i = 1 \frac{\partial [(c (v, α))_{i} ψ (v)]}{\partial v_{i}} - d \sum i = 1 D_{v_{i}}^{α} ψ (v)] + d \sum i = 1 \frac{\partial [\partial_{x_{i}} f (x) ϕ (x) ψ (v)]}{\partial v_{i}} - d \sum i = 1 \frac{\partial [\partial_{v_{i}} g (v) ϕ (x) ψ (v)]}{\partial x_{i}} .$

We can compute that

	$β d \sum i = 1 \frac{\partial [(c (v, α))_{i} ψ (v)]}{\partial v_{i}} - d \sum i = 1 D_{v_{i}}^{α} ψ (v)$

		$= - d \sum i = 1 \frac{\partial^{2}}{\partial v_{i}^{2}} D_{v_{i}}^{α - 2} ψ (v) - d \sum i = 1 D_{v_{i}}^{α} ψ (v)$
		$= d \sum i = 1 D_{v_{i}}^{2} D_{v_{i}}^{α - 2} ψ (v) - d \sum i = 1 D_{v_{i}}^{α} ψ (v) = 0,$

where we used the property $D^{2} u (x) = - \frac{\partial^{2}}{\partial x^{2}} u (x)$ (Proposition 1 in (Şimşekli, 2017)) and the semi-group property of the Riesz derivative $D^{a} D^{b} u (x) = D^{a + b} u (x)$ .

We can also compute that

	$d \sum i = 1 \frac{\partial [\partial_{x_{i}} f (x) ϕ (x) ψ (v)]}{\partial v_{i}} - d \sum i = 1 \frac{\partial [\partial_{v_{i}} g (v) ϕ (x) ψ (v)]}{\partial x_{i}}$
	$= d \sum i = 1 \partial_{x_{i}} f (x) ϕ (x) \partial_{v_{i}} ψ (v) - d \sum i = 1 \partial_{v_{i}} g (v) ψ (v) \partial_{x_{i}} ϕ (x)$
	$= - β d \sum i = 1 \frac{\partial f (x)}{\partial_{x_{i}}} ϕ (x) \frac{\partial g (v)}{\partial v_{i}} ψ (v) + β d \sum i = 1 \frac{\partial g (v)}{\partial_{v_{i}}} ψ (v) \frac{\partial f (x)}{\partial x_{i}} ϕ (x)$
	$= 0.$

Hence, $π (d x, d v) = \frac{e^{- β (f (x) + g (v))} d x d v}{\int_{R^{d} \times R^{d}} e^{- β (f (x^{'}) + g (v^{'}))} d x^{'} d v^{'}}$ is an invariant probability measure. The proof is complete.

Appendix B Proof of Theorem 3.1

We can compute that

(c (v, α))_{i} = \frac{D_{v_{i}}^{α - 2} (v_{i} e^{- \frac{}{1} 2 ∥ v ∥^{2}})}{e^{- \frac{1}{2} ∥ v ∥^{2}}} = e^{\frac{1}{2} v_{i}^{2}} D_{v_{i}}^{α - 2} (v_{i} e^{- \frac{1}{2} v_{i}^{2}}),

(24)

for every $1 \leq i \leq d$ .

Recall the definition of Fourier transform and its inverse:

F {f (x)} (ω) = \frac{1}{\sqrt{2 π}} \int_{- \infty}^{\infty} e^{- i x ω} f (x) d x, F^{- 1} {f (ω)} (x) = \frac{1}{\sqrt{2 π}} \int_{- \infty}^{\infty} e^{i x ω} f (ω) d ω .

(25)

Notice that the Fourier transform of $e^{- \frac{1}{2} x^{2}}$ is itself, i.e. $F {e^{- \frac{1}{2} x^{2}}} (ω) = e^{- \frac{1}{2} ω^{2}}$ , and moreover, $F {x^{n} f (x)} (ω) = i^{n} \frac{d^{n}}{d ω^{n}} {F {f (x)} (ω)}$ , and therefore,

F {x e^{- \frac{1}{2} x^{2}}} (ω) = - i ω e^{- \frac{1}{} 2 ω^{2}} .

(26)

Hence,

D_{x}^{α - 2} (x e^{- \frac{1}{2} x^{2}}) = F^{- 1} {- i ω | ω |^{α - 2} e^{- \frac{1}{2} ω^{2}}} (x) = \frac{- i}{\sqrt{2 π}} \int_{- \infty}^{\infty} ω | ω |^{α - 2} e^{- \frac{1}{2} ω^{2} + i ω x} d ω .

(27)

Furthermore, we can compute that

	$\frac{- i}{\sqrt{2 π}} \int_{- \infty}^{\infty} ω \| ω \|^{α - 2} e^{- \frac{1}{2} ω^{2} + i ω x} d ω$	$= \frac{- i}{\sqrt{2 π}} \int_{0}^{\infty} ω^{α - 1} e^{- \frac{1}{2} ω^{2} + i ω x} d ω + \frac{- i}{\sqrt{2 π}} \int_{- \infty}^{0} ω (- ω)^{α - 2} e^{- \frac{1}{2} ω^{2} + i ω x} d ω$
		$= \frac{- i}{\sqrt{2 π}} \int_{0}^{\infty} ω^{α - 1} e^{- \frac{1}{2} ω^{2} + i ω x} d ω + \frac{i}{\sqrt{2 π}} \int_{0}^{\infty} ω^{α - 1} e^{- \frac{1}{2} ω^{2} - i ω x} d ω$
		$= \sqrt{\frac{2}{π}} \int_{0}^{\infty} ω^{α - 1} sin (ω x) e^{- \frac{1}{2} ω^{2}} d ω .$

By the Taylor expansion of sine function, we get

	$\sqrt{\frac{2}{π}} \int_{0}^{\infty} ω^{α - 1} sin (ω x) e^{- \frac{1}{2} ω^{2}} d ω$	$= \sqrt{\frac{2}{π}} \int_{0}^{\infty} ω^{α - 1} \infty \sum k = 0 \frac{(- 1)^{k} (ω x)^{2 k + 1}}{(2 k + 1)!} e^{- \frac{1}{2} ω^{2}} d ω$
		$= \sqrt{\frac{2}{π}}$

Fractional Underdamped Langevin Dynamics: Retargeting SGD with Momentum under Heavy-Tailed Gradient Noise

Authors

First Exit Time Analysis of Stochastic Gradient Descent Under Heavy-Tailed Gradient Noise

The Heavy-Tail Phenomenon in SGD

Multiplicative noise and heavy tails in stochastic optimization

Non-Asymptotic Analysis of Fractional Langevin Monte Carlo for Non-Convex Optimization

On the Heavy-Tailed Theory of Stochastic Gradient Descent for Deep Neural Networks

Hausdorff Dimension, Stochastic Differential Equations, and Generalization in Neural Networks

Improved scalability under heavy tails, without strong convexity

1 Introduction

2 Technical Background & Related Work