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

Introduction

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

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

where $\mathbf{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 $\eta$:

While the Gaussianity assumption can be accurate in certain settings such as small networks (Martin & 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., 2019a; Zhang et al., 2019b)In two recent studies, Gürbüzbalaban et al. (2020) and Hodgkinson & Mahoney (2020) have shown that the stationary distribution of the stochastic gradient descent (SGD) algorithm can be indeed a heavy-tailed distribution depending on the choice of the step-size and the batch-size. On the other hand, in another recent study, Şimşekli et al. (2020) have provided generalization bounds for a general class of SDEs, including heavy- and light-tailed ones.. 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 $\alpha$-stable distribution ($\mathcal{S}\alpha\mathcal{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 ${\cal S}\alpha{\cal S}$, lies the parameter $\alpha\in(0,2]$, which determines the heaviness of the tail of the distribution. The tails get heavier as $\alpha$ gets smaller, the case $\alpha=2$ reduces to the Gaussian random variables. This is illustrated in Figure 1.

Şimşekli et al. (2019b, a) empirically illustrated that, in deep neural networks, the statistical structure of $U_{k}$ can be better captured by using an $\alpha$-stable distribution. With the assumption of $U_{k}$ being ${\cal S}\alpha{\cal S}$ distributed, the choice of Brownian motion will be no longer appropriate and should be replaced with an $\alpha$-stable Lévy motion, which motivates the following Lévy-driven SDE:

A more striking property of (7) was very recently revealed in a statistical physics study (Capała & 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., $\alpha=1$ and $\gamma\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 $\beta$) 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 $\mathbf{v}$, which are then directly transmitted to $\mathbf{x}$ due to the dynamics. Unless ‘tamed’, these updates create an hurling effect on $\mathbf{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.

Technical Background & Related Work

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 often occurs in statistical physics (Eliazar & Klafter, 2003), signal processing (Kuruoglu, 1999), and finance (Mandelbrot, 1997).

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

Here, $\phi(\mathbf{x})=\exp(-f(\mathbf{x}))$ and $\mathcal{D}$ denotes the fractional Riesz derivative (Riesz, 1949):

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

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

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

where $\mathbf{z}_{t}=(\mathbf{x}_{t},\mathbf{v}_{t})$, and $\phi(\mathbf{z})=e^{-f(\mathbf{x})-\frac{1}{2}\|\mathbf{v}\|^{2}}$. They showed that the invariant measure of the process has a density proportional to $\phi(\mathbf{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 & Dybiec, 2019).

Fractional Underdamped Langevin Dynamics

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

Let $c(\mathbf{v},\alpha)=((c(\mathbf{v},\alpha))_{i},1\leq i\leq d)$ has the following form:

One of the main features of FULD is that the fractional Riesz derivatives only appears in the drift $c$, which only depends on $\mathbf{v}$. This is highly in contrast with FHD (13), where the Riesz derivatives are taken over both $\mathbf{x}$ and $\mathbf{v}$, which is the source of intractability. Moreover, FULD enjoys the freedom to choose different kinetic energy functions $g(\mathbf{v})$. In the sequel, we will investigate two options for $g$, such that the drift $c$ can be analytically obtained.

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

In the next result, we will show that in this case, the drift $c$ admits an analytical solution.

Let $g(\mathbf{v})=\frac{1}{2}\|\mathbf{v}\|^{2}$. Then, for any $1\leq i\leq d$,

where $\Gamma$ is the gamma function and ${}_{1}F_{1}$ is the Kummer confluent hypergeometric function. In particular, when $\alpha=2$, we have $(c(\mathbf{v},\alpha))_{i}=v_{i}$.

We observe that the fractional dynamics (16) strictly extends the underdamped Langevin dynamics (6) as $c(\mathbf{v},2)=\mathbf{v}$.

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 & Platen, 1999); however, if we consider the integral form of ${}_{1}F_{1}$ (cf. (Abramowitz & Stegun, 1972)), we observe that the function

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

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 $\mathbf{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 $\mathbf{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 $\mathbf{v}$ to take large values, while still making sure that the drift $c$ in (15) admits an analytical form.

Let $e^{-g_{\alpha}(v)}$ be the probability density function of $\mathcal{S}\alpha\mathcal{S}(\frac{1}{\alpha^{1/\alpha}})$. Choose $\psi(\mathbf{v})=e^{-G_{\alpha}(\mathbf{v})}$ in (15), where $G_{\alpha}(\mathbf{v})=\sum_{i=1}^{d}g_{\alpha}(v_{i})$ for any $\mathbf{v}=(v_{1},\ldots,v_{d})$. Then,

It now follows from Theorem 3 that the FULD with $\alpha$-stable kinetic energy reduces to the following SDE:

It can be easily verified that $\nabla G_{\alpha}(\mathbf{v}_{t})=\mathbf{v}_{t}$ for $\alpha=2$, as $g_{2}(v)=\frac{1}{2}\log 2\pi+\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_{\alpha}$ itself admits a closed-form analytical formula only when $\alpha=1$ or $\alpha=2$, due to the properties of the ${\cal S}\alpha{\cal S}$ densities. Nevertheless, as $\nabla G_{\alpha}$ is based on one-dimensional ${\cal S}\alpha{\cal S}$ densities, it can be very accurately computed by using the recent methods developed in (Ament & O’Neil, 2018). On the other hand, in the next section, we will show that $\nabla G_{\alpha}$ is Lipschitz continuous for all $\alpha\in(0,2]$, which implies that under standard regularity conditions on $f$, the Boltzmann-Gibbs measure is the unique invariant measure of (19).

We visually inspect the behavior of $\nabla G_{\alpha}$ in Figure 2 for dimension one. We observe that, as soon as $\alpha<2$, $\nabla G_{\alpha}$ 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 $\mathbf{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 $\mathbf{v}$ by passing it through the asymptotically vanishing $\nabla G_{\alpha}$.

3 Euler discretization and weak convergence analysis

As visually hinted in Figure 2, the function $\nabla G_{\alpha}$ has strong regularity, which makes (19) to be potentially beneficial for practical implementations. Indeed, it is easy to verify that $\nabla G_{\alpha}$ is Lipschitz continuous for $\alpha=1$ and $2$, and in our next result, we show that this observation is true for any admissible $\alpha$, which is a desired property when discretizing continuous-time dynamics.

For $0<\alpha\leq 2$, the map $v\mapsto g^{\prime}_{\alpha}(v)$ is Lipschitz continuous, hence $\mathbf{v}\mapsto\nabla G_{\alpha}(\mathbf{v})$ is also Lipschitz continuous.

Accordingly we consider the following Euler-Maruyama discretization for (19):

where $S_{K}:=\sum_{k=1}^{K}\eta_{k}$ is the cumulative sum of the step-size sequence.

We note that Langevin-based algorithms have been used in the literature to obtain global convergence guarantees for non-convex optimization, see e.g. (Raginsky et al., 2017; Xu et al., 2018; Gao et al., 2018b; Zou et al., 2019; Nguyen et al., 2019). In particular, Nguyen et al. (2019) used an overdamped fractional Langevin dynamics for non-convex optimizations. The proposed model in our paper can also be used to study the non-convex optimizations and we expect that our underdamped dynamics may have improved theoretical guarantees compared to (Nguyen et al., 2019).

We now present the assumptions that imply our results.

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

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 1, Proposition 1, and (Panloup, 2008) (Theorem 2).

Assume that the gradient $\nabla f$ is Lipschitz continuous and has linear growth i.e., there exists $C>0$ such that $\|\nabla f(\mathbf{x})\|\leq C(1+\|\mathbf{x}\|)$ for all $\mathbf{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

4 Connections to existing approaches

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_{\alpha}$ for $\alpha<2$, as the components of $\mathbf{v}^{k}$ gets larger in magnitude, the update applied on $\mathbf{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 $\alpha$ decreases as opposed to being ‘clipped’ with a threshold.

Numerical Study

In this section, we will illustrate our theory on several experiments which are conducted in both synthetic and real-data settingsWe provide our implementation in https://github.com/umutsimsekli/fuld.. 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 ${\cal S}\alpha{\cal S}$ kinetic energy and from now on we will simply refer to FULD with ${\cal S}\alpha{\cal S}$ kinetic energy as FULD.

We first consider a one-dimensional synthetic setting, similar to the one considered in (Capała & 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 $\alpha$. For $\alpha\notin\{1,2\}$, we used the software given in (Ament & O’Neil, 2018) for computing $\nabla G_{\alpha}$.

Figure 3 illustrates the distribution of the samples generated by simulating the two dynamics. In this setup, we set $\beta=1$, $\eta=0.01$, $\gamma=10$ with number of iterations $K=50000$. We observe that, for $\alpha=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., $\alpha$ 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 $\alpha=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 & Dybiec, 2019)We note that the overdamped dynamics with the uncorrected drift exhibits a similar behavior to the one of the uncorrected underdamped dynamics with sufficiently large $\gamma$. . On the other hand, thanks to the correction brought by $\nabla G_{\alpha}$, 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 $\alpha$, 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 $(\mathbf{x}^{k})_{k=1}^{K}$ which were used for producing Figure 3. We observe that, while the iterates of UD are well-behaved for $\alpha=1.9$, the magnitude range of the iterates gets quite large when $\alpha$ is set to $1$. On the other hand, for both values of $\alpha$, FULD iterates are always kept in a reasonable range, thanks to the clipping-like effect of $\nabla G_{\alpha}$.

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 gradient noise (3) under the assumption that the noise can be well-modeled by using an ${\cal S}\alpha{\cal S}$ vector (see Section 3.4 for the explicit assumption).

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 $\alpha=2$, since $\nabla G_{2}(\mathbf{v})=\mathbf{v}$. Hence in this section, we will refer to SGDm as the special case of (22) with $\alpha=2$. On the other hand, in these experiments, directly computing $\nabla G_{\alpha}$ becomes impractical for $\alpha\notin\{1,2\}$, since the algorithms given in (Ament & O’Neil, 2018) become prohibitively slow with the increased dimension $d$. However, since $\nabla G_{\alpha}$ is based on the derivatives of the one-dimensional ${\cal S}\alpha{\cal S}$ densities $g_{\alpha}(v)$ (see Theorem 3), for $\alpha\in(1,2)$, we first precomputed the values of $g_{\alpha}(v)$ over a fine grid of $v\in$; then, during the SGDm recursion, we approximated $\nabla G_{\alpha}$ by linearly interpolating the values of $g_{\alpha}$ that are precomputed over this grid. We expect that, if the stochastic gradient noise can be well-approximated by using an ${\cal S}\alpha{\cal S}$ 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 $\eta=0.1$, $\gamma=0.1$ for MNIST, and $\gamma=0.9$ for CIFAR10. We run the algorithms for $K=10000$ iterations Since the scale of the gradient noise is proportional to $(\gamma/\beta)^{\frac{1}{\alpha}}$ (see (20)), in this setup, a fixed $\gamma$ implicitly determines $\beta$. . We measure the accuracy and the loss at every 100th iteration and we report the average of the last two measurements. Figures 5 and 6 show the results obtained on the MNIST dataset. We observe that, in most of the cases, setting $\alpha=1.75$ 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 $\approx 2\%$. We obtain a similar result on the CIFAR10 dataset, as illustrated in Figures 7 and 8. In most of the cases $\alpha=1.75$ performs better, with the maximum accuracy difference being $\approx 4.5\%$, implying the gradient noise can be approximated by an ${\cal S}\alpha{\cal S}$ 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. On the other hand, when the width was set to $512$, $\alpha=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 supplementary document.

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 $\alpha$-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 & Zhu (2018), state-dependent noise introduces challenging technical difficulties. Similarly, it has been illustrated that the tail-index $\alpha$ can depend on the state and different components of the noise can have a different $\alpha$ (Şimşekli et al., 2019a). Incorporating such state dependencies would be an important direction of future research. Finally, it has been shown that the heavy-tailed perturbations yield shorter escape times (Nguyen et al., 2019) in the overdamped dynamics, and extending such results to the underdamped case is still an open problem.

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

Proof of Theorem 1

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

where we used the property $\mathcal{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 $\mathcal{D}^{a}\mathcal{D}^{b}u(x)=\mathcal{D}^{a+b}u(x)$.

Therefore, it follows from (24), (25) and (26) that we have

Proof of Theorem 2

Recall the definition of Fourier transform and its inverse:

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

By the Taylor expansion of sine function, we get

where we used the identity $\int_{0}^{\infty}x^{a}e^{-\frac{1}{2}x^{2}}dx=2^{\frac{a-1}{2}}\Gamma(\frac{a+1}{2})$, for any given $a>-1$. Moreover, for any given $x,y>0$, we have the identity:

where ${}_{1}F_{1}$ is the Kummer confluent hypergeometric function. Therefore, we conclude that

By the identity $e^{x}\cdot_{1}F_{1}(a;b;-x)=_{1}F_{1}(b-a;b;x)$, we get

In particular, when $\alpha=2$, by applying the identity

Proof of Theorem 3

Let $\psi_{\alpha}(x)=e^{-g_{\alpha}(x)}$ be the probability density function of the symmetric $\alpha$-stable distribution $\mathcal{S}\alpha\mathcal{S}(\frac{1}{\alpha^{1/\alpha}})$ such that

Proof of Proposition 1

It is straightforward to verify that the result holds for the cases $\alpha=1$ and $\alpha=2$. Assume $\alpha\in(0,1)$ or $\alpha\in(1,2)$. Let $X$ be the unit symmetric $\alpha$-stable random variable defined by its characteristic function

By taking inverse Fourier transformation, its density $\psi_{\alpha}(x)=e^{-g_{\alpha}(x)}$ can be expressed as

Writing $e^{-itx}=\cos(tx)-i\sin(tx)$, we compute

where we used the fact that $\phi_{X}(t)$ and $\cos(tx)$ are even functions of $t$, whereas $\sin(tx)$ is an odd function of $t$. If we define,

In particular, since $|\cos(tx)|\leq 1$ and $|\sin(tx)|\leq 1$ this implies that

It is also well-known that a symmetric $\alpha$ stable random variable has a decay in its density satisfying $\psi_{\alpha}(x)\sim\frac{1}{|x|^{1+\alpha}}$ when $|x|$ is large. In fact, Wintner (1941) derived a large-$x$ expansion for $\psi_{\alpha}(x)$ when $0<\alpha<1$ and $x>0$. This expansion is equivalent to

(see eqn. (11) from (Montroll & Bendler, 1984)) where it can be seen from the Stirling’s approximation of the gamma function and the ratio test that the series converges absolutely. A similar absolutely convergent series sum (with exactly the same leading term) is also available in the literature for $\alpha\in(1,2)$ which says that

(see eqn. (3.58) from (Montroll & West, 1979)). By differentiating the series sum for $\psi_{\alpha}(x)$ with respect to $x$, we can express $\psi^{\prime}_{\alpha}(x)$ and $\psi^{\prime\prime}_{\alpha}(x)$ as a series sum. After a straightforward computation, we obtain

which implies from (40) that $g^{\prime\prime}_{\alpha}(x)\to 0$ as $x\to\infty$. This shows that $g^{\prime\prime}_{\alpha}(x)$ is bounded on the interval $[0,\infty)$. On the other hand, $\psi_{\alpha}(x)$ is an even function and therefore $g^{\prime\prime}_{\alpha}(x)$ is an even function satisfying $g^{\prime\prime}_{\alpha}(x)=g^{\prime\prime}_{\alpha}(-x)$. We conclude that $g^{\prime\prime}_{\alpha}(x)$ is bounded on the real line. This completes the proof. ∎

Proof of Corollary 1

By Proposition 1, we know that $\nabla G_{\alpha}$ is Lipschitz and by our hypthesis $\nabla f$ is also Lipschitz and has linear growth. Then the process (19) admits a unique invariant measure (cf. (Schertzer et al., 2001) Section 9), which is given by Theorem 1. The rest of the proof follows from (Panloup, 2008) (Theorem 2). ∎

Alternative forms of the drift function c𝑐c with the Gaussian kinetic energy

For some special values of $\alpha$, we can get alternative formulas for $(c(\mathbf{v},\alpha))_{i}$, $1\leq i\leq d$.

(1) $\alpha=\frac{3}{2}$. Using the identity ${}_{1}F_{1}(a;2a+1;z)=2^{2a-1}\Gamma(a+\frac{1}{2})e^{\frac{z}{2}}z^{\frac{1}{2}-a}(I_{a-\frac{1}{2}}(\frac{z}{2})-I_{a+\frac{1}{2}}(\frac{z}{2}))$, where $I_{a}(x)$ is the modified Bessel function of the first kind, we get

(2) $\alpha=\frac{1}{2}$. Using the identity ${}_{1}F_{1}(a;2a;z)=2^{2a-1}\Gamma(a+\frac{1}{2})z^{\frac{1}{2}-a}e^{\frac{z}{2}}I_{a-\frac{1}{2}}(\frac{z}{2})$, we get

Visual Illustrations

On the other hand, we visualize the conformal Hamiltonian field generated by this dynamics in Figure 10 for $f(x)=g_{1}(x)=-\log\frac{1}{\pi}\frac{1}{x^{2}+1}$. The figure shows that conformal Hamiltonian generated by the dynamics with $\alpha=2$ has a very slow concentration behavior towards the minimum at the origin, whereas this behavior is alleviated when $\alpha=1.7$ where the field concentrates faster.

Additional Experimental Results

In this section, we provide the additional experimental results that were mentioned in the main document for width $32$, $64$, and $512$.