# Asymmetric Heavy Tails and Implicit Bias in Gaussian Noise Injections

Alexander Camuto  Xiaoyu Wang  Lingjiong Zhu  Mert Gürbüzbalaban  Chris Holmes
Umut Şimşekli

# Asymmetric Heavy Tails and Implicit Bias in Gaussian Noise Injections
SUPPLEMENTARY DOCUMENT

Alexander Camuto  Xiaoyu Wang  Lingjiong Zhu  Mert Gürbüzbalaban  Chris Holmes
Umut Şimşekli

###### Abstract

Gaussian noise injections (GNIs) are a family of simple and widely-used
regularisation methods for training neural networks where one injects additive
or multiplicative Gaussian noise to the network activations at every iteration
of the optimisation algorithm, which is typically chosen as stochastic
gradient descent (SGD). In this paper we focus on the so-called ‘implicit
effect’ of GNIs, which is the effect of the injected noise on the dynamics of
SGD. We show that this effect induces an _asymmetric heavy-tailed noise_ on
SGD gradient updates. In order to model this modified dynamics, we first
develop a Langevin-like stochastic differential equation that is driven by a
general family of _asymmetric_ heavy-tailed noise. Using this model we then
formally prove that GNIs induce an ‘implicit bias’, which varies depending on
the heaviness of the tails and the level of asymmetry. Our empirical results
confirm that different types of neural networks trained with GNIs are well-
modelled by the proposed dynamics and that the implicit effect of these
injections induces a bias that degrades the performance of networks.

[pdftoc]hyperrefToken not allowed in a PDF string [xclr]xcolorIncompatible
color definition [pdftoc] [xclr]

##  1 Introduction

Figure 1:  Illustration of the effect of GNIs added to a network’s
activations. Each colored dot represents a neuron’s activations. We add GNIs,
represented as circles, to each layer’s activations bar the output layer.
Perhaps counter-intuitively, though the forward pass experiences Gaussian
noise, gradient updates in the backward pass experience heavy-tailed
asymmetric noise.

Noise injections are a family of methods that involve adding or multiplying
samples from a noise distribution to the weights and activations of a neural
network during training. The most commonly used distributions are Bernoulli
distributions and Gaussian distributions (Srivastava et al., 2014; Poole et
al., 2014) and the noise is most often inserted at the level of network
activations.

Though the regularisation conferred by Gaussian noise injections (GNIs) can be
observed empirically, and there have been many studies on the benefits of
noising data (Bishop, 1995; Cohen et al., 2019; Webb, 1994), the mechanisms by
which these injections operate are not fully understood. Recently, the
_explicit_ effect of GNIs, which is the added term to the loss function
obtained when marginalising out the injected noise, has been characterised
analytically (Camuto et al., 2020): it corresponds to a penalisation in the
Fourier domain which improves model generalisation.

Here we extend this analysis and focus on the implicit effect of GNIs. This is
the effect of the remaining noise that has been marginalised out when studying
the explicit effect. In particular we focus on the manner in which such noise
alters the dynamics of Stochastic Gradient Descent (SGD) (Wei et al., 2020;
Zhang et al., 2017). We show that the implicit effect is driven by an
_asymmetric heavy-tailed noise_ on the SGD gradient updates, as illustrated in
Figure 1.

To study the effect of this gradient noise, we model the dynamics of SGD for a
network experiencing GNIs by a stochastic differential equation (SDE) driven
by an _asymmetric heavy-tailed_ $\alpha$-stable noise. We demonstrate that
this model captures the dynamics of networks trained with GNIs and we show
that the stationary distribution of this process becomes arbitrarily distant
from the so-called Gibbs measure, whose modes exactly match the local minima
of the loss function, as the gradient becomes increasingly heavy-tailed and
asymmetric. Heavy-tailed and asymmetric gradient noise thus degrades network
performance and this suggests that models trained with the full effect of GNIs
will underperform networks trained solely with the explicit effect. We confirm
this experimentally for a variety of dense and convolutional networks.111See
https://github.com/alexander-camuto/asym-heavy-tails-bias-GNI for all code.

##  2 Background

Stable Distributions. The Generalised Central Limit Theorem (GCLT) (Gnedenko &
Kolmogorov, 1954) states that for a sequence of independent and identically
distributed (i.i.d.) random variables whose distribution has a power-law tail
with index $0<\alpha<2$, the normalised sum converges to a heavy-tailed
distribution called the $\alpha$-stable distribution ($\mathcal{S}_{\alpha}$)
as the number of summands grows. An $\alpha$-stable distributed random
variable $X$ is denoted by $X\sim\mathcal{S}_{\alpha}(\sigma,\theta,\mu)$,
where $\alpha\in(0,2]$ is the _tail-index_ , $\theta\in[-1,1]$ is the
_skewness_ parameter, $\sigma\geq 0$ is the _scale_ parameter, and
$\mu\in\mathbb{R}$ is called the _location_ parameter. The mean of $X$
coincides with $\mu$ if $\alpha>1$, and otherwise the mean of $X$ is
undefined. In this work, we always assume $\mu=0$. The parameter $\theta$ is a
measure of asymmetry. We say that $X$ follows a _symmetric_ $\alpha$-stable
distribution denoted as
$\mathcal{S}\alpha\mathcal{S}(\sigma)=\mathcal{S}_{\alpha}(\sigma,0,0)$ if
$\theta=0$ (and $\mu=0$). The parameter $\alpha\in(0,2]$ determines the tail
thickness of the distribution, and $\sigma>0$ measures the spread of $X$
around its mode. Note that when $\alpha<2$, $\alpha$-stable distributions have
heavy tails such that their moments are finite only up to the order $\alpha$.

The probability density function (p.d.f.) of an $\alpha$-stable random
variable, $\alpha\in(0,2]$, does not have a closed-form expression except for
a few special cases. When $\alpha=1$ and $\alpha=2$, the symmetric
$\alpha$-stable distribution reduces to the Cauchy and the Gaussian
distributions, respectively, (cf. Section 1.1. in (Samorodnitsky & Taqqu,
1994)). By their flexibility, such distributions can model many complex
stochastic phenomena for which exact analytic forms are intractable (Sarafrazi
& Yazdi, 2019; Fiche et al., 2013).

Lévy Processes. A Lévy process (motion) is a stochastic process with
independent, stationary increments. Formally, $\boldsymbol{\mathbf{L}}_{t}$ is
Lévy process if

  1. (i)

$\boldsymbol{\mathbf{L}}_{0}=0$ almost surely;

  2. (ii)

For any $t_{0}<t_{1}<\cdots<t_{N}$, the increments
$\boldsymbol{\mathbf{L}}_{t_{n}}-\boldsymbol{\mathbf{L}}_{t_{n-1}}$ are
independent, $n=1,2,\ldots,N$;

  3. (iii)

The difference $\boldsymbol{\mathbf{L}}_{t}-\boldsymbol{\mathbf{L}}_{s}$ and
$\boldsymbol{\mathbf{L}}_{t-s}$ have the same distribution;

  4. (iv)

$\boldsymbol{\mathbf{L}}_{t}$ is continuous in probability, i.e. for any
$\delta>0$ and $s\geq 0$,
$\mathbb{P}(|\boldsymbol{\mathbf{L}}_{t}-\boldsymbol{\mathbf{L}}_{s}|>\delta)\rightarrow
0$ as $t\rightarrow s$.

The $\alpha$-stable Lévy process is an important class of Lévy processes. In
particular, for $\alpha\in(0,2]$, let
$\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta}$ denote the $d$-dimensional
$\alpha$-stable Lévy process with independent components, i.e. each component
is an independent scalar $\alpha$-stable Levy motion (Duan, 2015) such that
$\boldsymbol{\mathbf{L}}_{t-s}^{\alpha,\theta}$ has the distribution
$\mathcal{S}_{\alpha}((t-s)^{1/\alpha},\theta,0)$ for any $s<t$.

Stochastic Gradient Descent and Differential Equations.  Let $\mathcal{D}$ be
a training dataset composed of data-label pairs of the form
$(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$, and let
$\boldsymbol{\mathbf{w}}\equiv\\{\boldsymbol{\mathbf{W}}_{1},\ldots,\boldsymbol{\mathbf{W}}_{L}\\}\in\mathbb{R}^{d}$
be the $d$ parameters of an $L$ layer neural network in vector form. When a
neural network operates on input data $\boldsymbol{\mathbf{x}}$, we obtain the
activations
${\boldsymbol{\mathbf{h}}}\equiv\\{{\boldsymbol{\mathbf{h}}}_{0},\ldots,{\boldsymbol{\mathbf{h}}}_{L-1}\\}$,
where $\boldsymbol{\mathbf{h}}_{0}=\boldsymbol{\mathbf{x}}$ and
${\boldsymbol{\mathbf{h}}}\in\mathbb{R}^{n_{0}+\dots+n_{L}}$ where $n_{i}$ is
the number of neurons in the $i^{\mathrm{th}}$ layer: we consider a non-
linearity $\kappa:\mathbb{R}\to\mathbb{R}$,
$\boldsymbol{\mathbf{h}}_{i}(\boldsymbol{\mathbf{x}})=\kappa(\boldsymbol{\mathbf{W}}_{i}\boldsymbol{\mathbf{h}}_{i-1}(\boldsymbol{\mathbf{x}}))$,
where $\kappa$ is applied element-wise to each coordinate of its argument. In
supervised settings, our objective is to find the optimal parameters
$\boldsymbol{\mathbf{w}}_{*}$ that minimise the negative log-likelihood $-\log
p_{\boldsymbol{\mathbf{w}}}(\boldsymbol{\mathbf{y}}|\boldsymbol{\mathbf{x}})$
of the labels $\boldsymbol{\mathbf{y}}$, given the parameters
$\boldsymbol{\mathbf{w}}$ and data $\boldsymbol{\mathbf{x}}$:

$\displaystyle\boldsymbol{\mathbf{w}}_{*}=\operatorname*{arg\min}_{\boldsymbol{\mathbf{w}}}\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}}),$
$\displaystyle\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}}):=-\mathbb{E}_{\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}}\sim\mathcal{D}}\left[\log
p_{\boldsymbol{\mathbf{w}}}(\boldsymbol{\mathbf{y}}|\boldsymbol{\mathbf{x}})\right]\,.$
(2.1)

SGD and its variants are the most prevalent optimisation routines that
underpin the training of very large neural networks. Under SGD, we estimate
equation (2.1) by sampling a _random_ mini-batch of data-label pairs
$\mathcal{B}\subset\mathcal{D}$,

$\displaystyle\mathcal{L}(\mathcal{B};\boldsymbol{\mathbf{w}})=-\mathbb{E}_{\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}}\sim\mathcal{B}}\left[\log
p_{\boldsymbol{\mathbf{w}}}(\boldsymbol{\mathbf{y}}|\boldsymbol{\mathbf{x}})\right]\approx\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}}).$
(2.2)

SGD optimises this equation and approximates $\boldsymbol{\mathbf{w}}_{*}$
using iterative parameter updates. At training step $k$

$\displaystyle\boldsymbol{\mathbf{w}}_{k+1}=\boldsymbol{\mathbf{w}}_{k}-\eta\nabla\mathcal{L}(\mathcal{B}_{k+1};\boldsymbol{\mathbf{w}}_{k}),$
(2.3)

where $\eta$ is the step-size for updates (the network’s learning rate)
(Robbins & Monro, 1951; Ruder, 2016).

Studying the dynamics of SGD allows us to understand the subtle effects that
batching may have on neural network training. The similarities between the SGD
algorithm and Langevin diffusions (Roberts & Stramer, 2002) have inspired many
studies modelling the dynamics of SGD using stochastic differential equations
(SDEs) under different noise conditions (Agapiou et al., 2014; Şimşekli et
al., 2019; Raginsky et al., 2017; Gao et al., 2018, 2020; Jastrzȩbski et al.,
2017; Li et al., 2017). In this approach, one models the discrete SGD updates
(2.3), as the discretisation of a continuous-time stochastic process, making
assumptions about the properties of the ‘noise’ that drives this process
(Mandt et al., 2016; Jastrzȩbski et al., 2017). This noise stems from the
stochasticity in approximating the ‘true’ gradient over the dataset
$\nabla\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}}_{k})$ with that of a
mini-batch $\mathcal{B}$,
$\nabla\mathcal{L}(\mathcal{B};\boldsymbol{\mathbf{w}}_{k})$. We denote this
noise as,

$\displaystyle
U_{k+1}(\boldsymbol{\mathbf{w}}):=\nabla\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}}_{k})-\nabla_{\boldsymbol{\mathbf{w}}_{k}}\mathcal{L}(\mathcal{B}_{k+1};\boldsymbol{\mathbf{w}}_{k}).$
(2.4)

The most prevalent assumption is that the gradient noise admits a multi-
variate Gaussian noise:
$U_{k}(\boldsymbol{\mathbf{w}})\sim\mathcal{N}(0,\sigma^{2}\boldsymbol{\mathbf{I}})$.
This is rationalised by the central limit theorem, as the sum of estimation
errors in equation (2.4) is approximately Gaussian for sufficiently large
batches. Under this assumption, we can rewrite the SGD parameter update as:

$\displaystyle\boldsymbol{\mathbf{w}}_{k+1}=\boldsymbol{\mathbf{w}}_{k}-\eta\nabla\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}}_{k})+\eta
U_{k+1}(\boldsymbol{\mathbf{w}}).$ (2.5)

Then, one obtains the following continuous-time SDE to approximate the
gradient updates (Welling & Teh, 2011; Mandt et al., 2016; Jastrzȩbski et al.,
2017):

$\displaystyle
d\boldsymbol{\mathbf{w}}_{t}=-\nabla\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}}_{t})dt+\sqrt{\eta\sigma^{2}}d\boldsymbol{\mathbf{B}}_{t},$
(2.6)

where $\boldsymbol{\mathbf{B}}_{t}$ is the Brownian motion in $\mathbb{R}^{d}$
and $\sigma$ is the assumed noise variance for $U$. However, recent work
suggests that the Gaussian assumption might not be always appropriate
(Gürbüzbalaban et al., 2020; Hodgkinson & Mahoney, 2020), and connectedly, the
gradient noise is observed to be heavy-tailed in different settings (Şimşekli
et al., 2019; Zhou et al., 2020). Under this noise assumption, the
corresponding SDE is such that $\boldsymbol{\mathbf{B}}_{t}$ in (2.6) is
replaced with the _symmetric_ stable process
$\boldsymbol{\mathbf{L}}_{t}^{\alpha,0}$ where $\theta=0$:

$\displaystyle
d\boldsymbol{\mathbf{w}}_{t}=-\nabla\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}}_{t})dt+\eta^{(\alpha-1)/\alpha}\sigma
d\boldsymbol{\mathbf{L}}_{t}^{\alpha,0}.$ (2.7)

Gaussian Noise Injections. GNIs are regularisation methods that consist of
injecting Gaussian noise to the network activations. More precisely, let
$\boldsymbol{\mathbf{\epsilon}}$ be a collection of ‘noise vectors’ injected
to the network activations at each layer except the final layer:
$\boldsymbol{\mathbf{\epsilon}}\equiv\\{{\boldsymbol{\mathbf{\epsilon}}}_{0},\ldots,{\boldsymbol{\mathbf{\epsilon}}}_{L-1}\\}$,
where $\boldsymbol{\mathbf{\epsilon}}_{i}\in\mathbb{R}^{n_{i}}$,
$\boldsymbol{\mathbf{\epsilon}}\in\mathbb{R}^{n_{0}+\dots+n_{L-1}}$ and
$n_{i}$ is the number of neurons in the $i^{\mathrm{th}}$ layer. We have two
values for an activation: the soon-to-be noised value
$\widehat{\boldsymbol{\mathbf{h}}}_{i}$, and the subsequently noised value
$\widetilde{\boldsymbol{\mathbf{h}}}_{i}$. For a multi-layer perceptron (MLP),

$\displaystyle\widehat{\boldsymbol{\mathbf{h}}}_{i}(\boldsymbol{\mathbf{x}})=\kappa\left(\boldsymbol{\mathbf{W}}_{i}\widetilde{\boldsymbol{\mathbf{h}}}_{i-1}(\boldsymbol{\mathbf{x}})\right)\,,\qquad\widetilde{\boldsymbol{\mathbf{h}}}_{i}=\widehat{\boldsymbol{\mathbf{h}}}_{i}\circ\boldsymbol{\mathbf{\epsilon}}_{i}\,,$
(2.8)

where $\circ$ is some element-wise operation (e.g., addition or
multiplication). For additive GNIs typically,
$\boldsymbol{\mathbf{\epsilon}}_{i}\sim\mathcal{N}\left(0,\sigma_{i}^{2}\boldsymbol{\mathbf{I}}\right)$,
and for multiplicative GNIs we often use
$\boldsymbol{\mathbf{\epsilon}}_{i}\sim\mathcal{N}\left(1,\sigma_{i}^{2}\boldsymbol{\mathbf{I}}\right)$,
where $\mathcal{N}$ is the Gaussian distribution. The ‘accumulated’ noise at
each layer, induced by the noise injected to the layer and in previous layers
is

$\displaystyle\mathcal{E}_{i}(\boldsymbol{\mathbf{x}};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})=\widetilde{\boldsymbol{\mathbf{h}}}_{i}(\boldsymbol{\mathbf{x}})-\boldsymbol{\mathbf{h}}_{i}(\boldsymbol{\mathbf{x}}).$
(2.9)

Given that GNIs are commonly used as regularisation methods (Camuto et al.,
2020; Dieng et al., 2018; Srivastava et al., 2014; Poole et al., 2014; Kingma
et al., 2015; Bishop, 1995), our goal is to understand better the mechanisms
by which they affect neural networks.

##  3 The Implicit Effect of GNIs

Recently, Camuto et al. (2020) showed that the effect of GNIs on the cost
function can be expressed as a term
$\Delta\mathcal{L}(\mathcal{B};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})$
that is added to the loss, i.e.,

$\displaystyle\widetilde{\mathcal{L}}(\mathcal{B};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}}):=\mathcal{L}(\mathcal{B};\boldsymbol{\mathbf{w}})+\Delta\mathcal{L}(\mathcal{B};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})\,,$
(3.1)

where $\widetilde{\mathcal{L}}$ is the modified loss that SGD ultimately aims
to minimise. The term $\Delta\mathcal{L}$ can be further broken down into
explicit and implicit effects, as described in Section 1.

We build on the approach of Camuto et al. (2020) and define the explicit
effect as the additional term obtained on the loss when we marginalise out the
noise we have injected. It offers a consistently positive objective for
gradient descent to optimise and we denote it as
$\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}(\Delta\mathcal{L}(\mathcal{B};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}}))$.
The implicit effect is then the remainder of the terms marginalised out in the
explicit effect:

$\displaystyle
E_{\mathcal{L}}(\mathcal{B};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}}):=\Delta\mathcal{L}(\mathcal{B};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})-\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}(\Delta\mathcal{L}(\mathcal{B};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})).$
(3.2)

While the explicit effect focuses on the consistent and non-stochastic
regularisation induced by GNIs, the implicit effect instead studies the effect
of the inherent stochasticity of GNIs. This term does not offer a consistent
objective for SGD to minimise. Rather we show that it affects neural network
training by way of the heavy-tailed and skewed noise it induces on gradient
updates.

Subsequently, we first characterise the tail properties of the noise
accumulated during the forward pass. As this accumulated noise defines the
implicit effect we can then use this result to show that gradients induced by
the implicit effect are heavy-tailed and skewed and apt to be modelled by
heavy-tailed and asymmetric $\alpha$-stable noise.

The Properties of the Accumulated Noise. Before considering the properties of
the implicit effect gradients, we first need to study the noise that is
accumulated during the forward pass of a neural network experiencing GNIs.
Here we use asymptotic analysis to bound the moments of this noise, allowing
us to characterise the tail properties of gradients using the broad class of
sub-Weibull distributions that include a range of heavy-tailed distributions
(Vladimirova et al., 2019, 2020; Kuchibhotla & Chakrabortty, 2018).

######  Definition 3.1.

(Asymptotic order) A positive sequence $a_{m}$ is of the same order as another
positive sequence $b_{m}$ ($a_{m}\lesssim b_{m}$) if $\exists C>0$ such that
$\frac{a_{m}}{b_{m}}\leq C~{}\forall m\in\mathbb{N}$. $a_{m}$ is of the same
order of magnitude as $b_{m}$ ($a_{m}\asymp b_{m}$, i.e. ‘asymptotically
equivalent’) if there exist some $c,C>0$ such that:
$c\leq\frac{a_{m}}{b_{m}}\leq C$ for any $m\in\mathbb{N}$.

######  Definition 3.2.

(Sub-Weibull distributions), : We say that a random variable $X$ is sub-
Weibull (Vladimirova et al., 2020; Kuchibhotla & Chakrabortty, 2018) with tail
parameter $r$ if $\|X\|_{m}\lesssim m^{r},r>0$, where
$\|X\|_{m}:=\operatorname{\mathbb{E}}[|X|^{m}]^{\frac{1}{m}}$. In this case,
we write $X\sim\mathrm{subW}(r)$. Such distributions satisfy the tail bound
$\mathbb{P}(|X|>x)\leq 2e^{-(\frac{x}{C})^{\frac{1}{r}}}$, where $C>0$ is some
constant. Note that for $r=\frac{1}{2}$ and $r=1$ we recover the sub-Gaussian
and sub-exponential families and that if $X\sim\mathrm{subW}(r^{\prime})$,
then $X\sim\mathrm{subW}(r)$ for $r>r^{\prime}$.

As $r$ increases the tail distribution becomes heavier-tailed. To simplify our
analysis we consider activation functions $\phi$ that obey the extended
envelope property.

######  Definition 3.3 (Extended Envelope Property (Vladimirova et al.,
2019):).

A non-linear function $\kappa:\mathbb{R}\to\mathbb{R}$ is said to obey the
extended envelope property if $\exists~{}c_{1},c_{2}\geq 0,~{}d_{1},d_{2}\geq
0$ such that:

  * •

$|\kappa(x)|\geq c_{1}+d_{1}|x|$, for any $x\in\mathbb{R}^{+}$ or
$x\in\mathbb{R}^{-}$;

  * •

$|\kappa(x)|\leq c_{2}+d_{2}|x|$, for any $x\in\mathbb{R}$.

Activation functions that obey this property, such as $\mathrm{ReLU}$, are
broadly moment preserving (Vladimirova et al., 2020). Using this property, we
can characterise the moments for the $l^{\mathrm{th}}$ noised activation in a
layer $i$, $\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})$.

######  Lemma 3.1.

For feed-forward neural networks with an activation function $\phi$ that obeys
the extended envelope property, the noised activations at each layer $i<L-1$,
resulting from additive-GNIs $\boldsymbol{\mathbf{\epsilon}}$ obey

$\displaystyle\left\|\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})\right\|_{m}\lesssim\sqrt{m},\quad\text{for
any $m\geq 1$};\ l=1,\dots,n_{i}\,.$

For multiplicative-GNIs we have,

$\displaystyle\left\|\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})\right\|_{m}\lesssim
m^{\frac{i+1}{2}},\quad\text{for any $m\geq 1$};\ l=1,\dots,n_{i}\,,$

where $n_{i}$ is the dimensionality of the $i^{\mathrm{th}}$ layer.

Lemma 3.1 shows that when the injected noise is additive, the noised
activations at each layer will have (sub)-Gaussian tails. By equation (2.9),
the accumulated noise
$\mathcal{E}_{i}(\boldsymbol{\mathbf{x}};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})$
will also have sub-Gaussian tails as the non-noised activations
$\boldsymbol{\mathbf{h}}(\boldsymbol{\mathbf{x}})$ are deterministic and do
not affect the asymptotic relationships of moments. See Figure A.1 of the
Appendix for a demonstration that the activations experience Gaussian-like
noise for additive-GNIs. For multiplicative-GNIs the noise at each layer,
except the input layer which experiences Gaussian noise, behaves with a sub-
Weibull tail. We demonstrate this behaviour in Figure A.2 of the Appendix.

We can now study the properties of the gradient noise induced by the implicit
effect by taking the gradients of this forward pass noise.

Kurtosis of The Gradient Noise. We characterise the gradient noise
corresponding to $W_{i,l,j}$, the weight that maps from neuron $l$ in layer
$i-1$ to neuron $j$ in layer $i$.

######  Theorem 3.1.

Consider a feed-forward neural network with an activation function $\phi$ that
obeys the extended envelope property and a cross-entropy or mean-squared-error
cost (see Appendix B). The gradient noise from additive-GNIs
$\boldsymbol{\mathbf{\epsilon}}$, has zero mean and has moments that obey for
a pair $(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$:

$\displaystyle\left\|\frac{\partial
E_{\mathcal{L}}((\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}});\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})}{\partial
W_{i,l,j}}\right\|_{m}\lesssim m,\quad\text{for any $m\geq 1$},$

where
$E_{\mathcal{L}}((\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}});\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})$
is defined in (3.2). For multiplicative GNIs, we have

$\displaystyle\left\|\frac{\partial
E_{\mathcal{L}}((\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}});\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})}{\partial
W_{i,l,j}}\right\|_{m}\lesssim m^{\frac{L+i}{2}},\quad\text{for any $m\geq
1$},$ $\displaystyle i=1,\dots,L;\ l=1,\dots,n_{i-1};\ j=1,\dots,n_{i}.$

For the additive noise, these bounds infer that the gradient noise at each
layer will have sub-exponential tails. For the multiplicative case, gradient
noise will be sub-Weibull, with a tail parameter that increases with $i$ the
layer index. Unlike the forward pass, which experienced noise bounded in its
tails by a Gaussian, the backward pass experiences noise that is bounded in
its tails by heavy-tailed Weibull distributions with tail parameter $r\geq 1$.

######  Remark 3.1.

The bounds defined by Lemma 3.1 and Theorem 3.1 become an asymptotic
equivalence ($\asymp$) in the case of 1-D data and 1-neuron-wide neural
networks, i.e. the bounds are maximally tight in this case.

Our result applies to the gradient for a single pair
$(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$. During SGD, we take the
mean gradient across a batch $\mathcal{B}$ of size $B$. To study the tail
distribution of this mean gradient, we restate in a simplified manner
Kuchibhotla & Chakrabortty (2018)’s generalisation of the Bernstein inequality
for zero-mean sub-Weibull random variables in Theorem 3.2.

######  Theorem 3.2 (Theorem 3.1 in Kuchibhotla & Chakrabortty (2018)).

Let $X_{1},\dots,X_{B}$ be independent mean-zero sub-Weibull random variables
with a shared tail parameter $p\geq 1$. Then, for every $x\geq 0$, we have

$\displaystyle\mathbb{P}\left\\{\left|(1/B)\sum\nolimits_{i=1}^{B}X_{i}\right|\geq
x\right\\}$ $\displaystyle\leq
2\exp\left[-\min\left(\frac{Bx^{2}}{C\sum^{B}_{i=1}\|X_{i}\|^{2}_{\psi
2}},\frac{Bx^{\frac{1}{\theta}}}{L\max_{i}\|X_{i}\|_{\psi\frac{1}{p}}}\right)\right],$

where $C,L>0$ are constants that depend on the tail parameter and where

$\|X\|_{\psi\frac{1}{p}}=\inf\left\\{\nu\geq
0:\operatorname{\mathbb{E}}\left[(|X|/\nu)^{\frac{1}{p}}\leq
1\right]\right\\}$ (3.3)

is the ‘sub-Weibull norm’.

This theorem states that the the tails of the mean of zero-mean i.i.d. sub-
Weibull random variables are produced by a single variable, say $X_{i}$, with
the maximal sub-Weibull norm $\|X_{i}\|_{\psi\frac{1}{p}}=\inf\\{\nu\geq
0:\operatorname{\mathbb{E}}[(|X_{i}|/\nu)^{\frac{1}{p}}\leq 1]\\}$, i.e., the
one with the heaviest tails. Assuming gradients are independent across data
points, we can use this inequality to bound the tail probability for the sum
of our zero-mean gradients. If a single gradient is sufficiently heavy-tailed,
then the mean across the batch will also be heavy-tailed.

In Figure 2 we show that the implicit effect gradient, averaged over the
entire dataset $\mathcal{D}$ (i.e., the largest batch-size possible), is
heavy-tailed, unlike the forward pass. To calculate these gradients, we
estimate the explicit regulariser in equation (3.2) using Monte Carlo
sampling,
$\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}(\Delta\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}}))\approx\frac{1}{M}\sum_{m=0}^{M}\Delta\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}}_{m})$,
similarly to (Wei et al., 2020; Camuto et al., 2020). We show these results
for multiplicative-GNIs in Figure A.3 of the Appendix. In this setting as
well, the backward pass experiences heavy-tailed noise from GNIs.

Skewness of The Gradient Noise. In the proof of Theorem 3.1 we decompose
$\partial E_{\mathcal{L}}(\cdot)/\partial W_{i,l,j}$ as $(\partial
E_{\mathcal{L}}(\cdot)/\partial\widetilde{h}_{i,j})\cdot(\partial\widetilde{h}_{i,j}/\partial
W_{i,l,j})$, where $\widetilde{h}_{i,j}$ is the (noised) activation of the
$j^{\mathrm{th}}$ neuron in the $i^{\mathrm{th}}$ layer. Both these
derivatives are likely to be skewed due to the asymmetry of the activation
functions and their gradients. Ignoring potential correlations between
variables, we observe that product of two zero-mean skewed independent
variables $X$ and $Y$ is also skewed
$\left|\mathrm{skew}(XY)\right|=\left|\operatorname{\mathbb{E}}\left[X^{3}\right]\operatorname{\mathbb{E}}\left[Y^{3}\right]\right|>0$.
The skewness of the derivatives in the backward pass will induce skewed
gradient noise, as seen in Figure 2. Though correlations between gradients
could also cause skewness, we show that this is not the case in Appendix C.

##  4 An SDE Model for SGD with GNIs.

In this section, we will analyse the effects of the skewed and heavy-tailed
noise on the SGD dynamics. Recall that the modified loss function by the GNIs
is the sum of the explicit regulariser and the original loss over the dataset
$\mathcal{D}$,

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}\left(\widetilde{\mathcal{L}}(\mathcal{D};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})\right)=\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}})+\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}\left(\Delta\mathcal{L}(\mathcal{D};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})\right),$

and the SGD recursion takes the following form:

$\displaystyle\boldsymbol{\mathbf{w}}_{k+1}=\boldsymbol{\mathbf{w}}_{k}-\eta\nabla\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}\left(\widetilde{\mathcal{L}}(\mathcal{D};\boldsymbol{\mathbf{w}}_{k},\boldsymbol{\mathbf{\epsilon}})\right)+\eta
U_{k}(\boldsymbol{\mathbf{w}}),$ (4.1)

where $U_{k+1}(\boldsymbol{\mathbf{w}})$ is given as follows:

$\displaystyle\nabla\left[\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}\left(\widetilde{\mathcal{L}}(\mathcal{D};\boldsymbol{\mathbf{w}}_{k},\boldsymbol{\mathbf{\epsilon}})\right)\right]$
(4.2)
$\displaystyle\phantom{as}-\nabla\left[\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}\left(\widetilde{\mathcal{L}}(\mathcal{B}_{k+1};\boldsymbol{\mathbf{w}}_{k},\boldsymbol{\mathbf{\epsilon}})\right)+E_{\mathcal{L}}(\mathcal{B}_{k+1};\boldsymbol{\mathbf{w}}_{k},\boldsymbol{\mathbf{\epsilon}})\right].$

To ease the notation, let us denote the modified loss function by
$f(\boldsymbol{\mathbf{w}}):=\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}\left(\widetilde{\mathcal{L}}(\mathcal{D};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})\right)$.

######  Remark 4.1.

Note that when the gradients are computed over the whole dataset
$\mathcal{D}$, the gradient noise solely stems from the implicit effect,
$U_{k}(\boldsymbol{\mathbf{w}})=-\nabla\left[E_{\mathcal{L}}(\mathcal{D};\boldsymbol{\mathbf{w}}_{k},\boldsymbol{\mathbf{\epsilon}})\right]$.

Figure 2: We measure the skewness and kurtosis at initialisation of the
gradients noise accrued on networks weights during the backward pass for
additive-GNIs. The model is a 4-layer-256-unit-wide MLP trained to regress
$\lambda(x)=\sum_{i}\sin(2\pi q_{i}x+\phi(i))$ with
$q_{i}\in(5,10,\dots,45,50),x\in\mathbb{R}$. We plot the probability density
function (p.d.f.) of positive samples, comparing against half-normal and half-
Cauchy distributions. Each blue point represents the mean gradient noise over
the entire dataset $\mathcal{D}$ of an individual weight in a layer $i$. This
gradient noise is skewed and heavy-tailed, with a p.d.f. that is more Cauchy-
like than Gaussian.

Recent studies have proven that heavy-tailed behaviour can already emerge in
stochastic optimisation (without GNIs) (Hodgkinson & Mahoney, 2020;
Gürbüzbalaban et al., 2020; Gürbüzbalaban & Hu, 2020), which can further
result in an overall heavy-tailed behaviour in the gradient noise, as
empirically reported in (Şimşekli et al., 2019; Zhang et al., 2020; Zhou et
al., 2020). Accordingly, (Şimşekli et al., 2019; Zhou et al., 2020) proposed
modelling the gradient noise by using a centred _symmetric_ $\alpha$-stable
noise, which then paved the way for modelling the SGD dynamics by using an SDE
driven by a symmetric $\alpha$-stable process (see (2.7)).

We take a similar route for modelling the trajectories of SGD with GNIs;
however, due to the skewness arising from the GNIs, the symmetric noise
assumption is not appropriate for our purposes. Hence, we propose modelling
the skewed gradient noise (4.2) by using an _asymmetric_ $\alpha$-stable
noise, which aims at modelling both the heavy-tailed behaviour and the
asymmetries at the same time. In particular, we consider an SDE driven by an
asymmetric stable process and its Euler discretisation as follows:

$\displaystyle d\boldsymbol{\mathbf{w}}_{t}=-\nabla
f(\boldsymbol{\mathbf{w}}_{t})dt+\varepsilon
d\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta},$ (4.3)
$\displaystyle\boldsymbol{\mathbf{w}}_{k+1}=\boldsymbol{\mathbf{w}}_{k}-\eta_{k+1}\nabla
f(\boldsymbol{\mathbf{w}}_{k})+\varepsilon\eta_{k+1}^{1/\alpha}\Delta\boldsymbol{\mathbf{L}}^{\alpha,\theta}_{k+1},$
(4.4)

where $\theta=(\theta_{i},1\leq i\leq d)$ is the $d$-dimensional skewness and
each coordinate can have its idiosyncratic skewness $\theta_{i}$.
$\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta}=(L_{t}^{\alpha,\theta_{1}},\ldots,L_{t}^{\alpha,\theta_{d}})$
is a $d$-dimensional asymmetric $\alpha$-stable Lévy process with independent
components, and $\varepsilon$ encapsulates all the scaling parameters.
Furthermore, $(\eta_{k})_{k}$ denotes the sequence of step-sizes, which can be
taken as constant or decreasing, and finally
$(\Delta\boldsymbol{\mathbf{L}}^{\alpha,\theta}_{k})_{k}$ is a sequence of
i.i.d. random vectors where each component of
$\Delta\boldsymbol{\mathbf{L}}^{\alpha,\theta}_{k}$ is i.i.d. with
$\mathcal{S}_{\alpha}(1,\theta_{i},0)$. We then propose the discretised
process (4.4) as a proxy to the original recursion (4.1) and we will directly
analyse the theoretical properties of (4.4). Note that our approach strictly
extends (Şimşekli et al., 2019), which appears as a special case when
$\theta=0$.

Before deriving our theoretical results, we first verify empirically that the
proxy dynamics (4.4) are indeed a good model for representing (4.1). We
ascertain that $\Delta\boldsymbol{\mathbf{L}}^{\alpha,\theta}_{k}$ is
sufficiently general in the sense that it can capture the gradient noise
induced by the implicit effect even when there is no batching noise (when the
batch size approaches the size of the dataset for example, see Remark 4.1). As
a first line of evidence, in Figure A.4 of the Appendix we model $\nabla
E_{\mathcal{L}}(\cdot)$ as being drawn from a univariate
$\mathcal{S}_{\alpha}$. The equivalent $\mathcal{S}_{\alpha}$ distributions
are skewed ($|\theta|>0$) and heavy-tailed ($\alpha<2$), demonstrating that
$\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta}$ captures the core properties of
the implicit effect gradients highlighted in Section 3. To illustrate this
more clearly, in Figure 3 we use symmetric sub-Weibull distributions with
$r=0.8,1,2$, and fit $\mathcal{S}_{\alpha}$ and Gaussians ($\mathcal{N}$)
using maximum likelihood (MLE) from $10^{4}$ samples. We plot the MLE
densities, and clearly MLE $\mathcal{S}_{\alpha}$ ($\alpha<2$) better model
the tails of the sub-Weibulls than a Gaussian distribution, even for $r=0.8$
which is close to Gaussian tails ($r=0.5$). The $\mathcal{S}_{\alpha}$
modelling of the tails improves as $r$ increases. This further illustrates the
appropriateness of our noise model.

Figure 3: [left col.] sub-Weibull pdfs ($r$=$0.8,1.0,2.0$), and pdfs of MLE
fitted $\mathcal{S}_{\alpha}$ and $\mathcal{N}$. [right col.] pdf in the tails
($x>4$ and $x>10$, note log y-axis).

Figure 4:  We train the networks in Figure 2 on $\widetilde{\mathcal{L}}_{M}$
(4.5) for the sinusoidal toy-data with additive GNIs [left] and multiplicative
GNIs [right]. We fit univariate $\mathcal{S}_{\alpha}$ and univariate
$\mathcal{N}$ via maximum likelihood (Nolan, 2001) to $\nabla
E_{\mathcal{L}}(\cdot)$ for $M=1$ models at each training step. We add draws
from these distributions to the gradients of $M=16$ models and plot the
training loss ($\mathcal{L}_{\mathrm{train}}$). Shading is the standard
deviation over 5 random seeds.

In Figure 4, we further show that the gradient noise of the implicit effect
and gradient noise drawn from an equivalent $\mathcal{S}_{\alpha}$
distribution will have similar effects on gradient descent. We sample $M$ GNI
samples and evaluate:

$\displaystyle\widetilde{\mathcal{L}}_{M}(\mathcal{D};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})=(1/M)\sum\nolimits_{m=0}^{M}\widetilde{\mathcal{L}}(\mathcal{D};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}}_{m}).$
(4.5)

The objective is over the entire dataset such that we eliminate noise from the
batching process. $M$ allows us to control the ‘degree’ to which the implicit
effect is marginalised out. $M=1$ corresponds to the usual training with GNI
and larger values of $M$ mimic the effects of marginalising out the implicit
effect. We model $\nabla E_{\mathcal{L}}$ for $M=1$ as being drawn from a
univariate $\mathcal{S}_{\alpha}$ or a univariate normal distribution and
estimate distribution parameters using maximum likelihood estimation, as in
(Nolan, 2001), at each training iteration. We add draws from the estimated
distributions to the gradients of $M=16$ models to mimic the combined implicit
and explicit effects. $M=16$ models with the added $\mathcal{S}_{\alpha}$
noise have the same training path as $M=1$ models, whereas those with Gaussian
gradient noise do not. Thus, $\mathcal{S}_{\alpha}$ distributions are able to
faithfully capture the dynamics induced by the implicit effect on gradient
descent.

In these same experiments $M=16$ models outperform $M=1$ models on training
data. We refine this study for a greater range of $M$ values in Figure 5. As
$M$ increases, performance of models on training data improves gradually,
suggesting that the implicit effect degrades performance. Further, in Figure
4, $M=16$ models trained with Gaussian noise added to gradients outperform
$M=1$ models, suggesting that the heavy-tails and skew of the implicit effect
gradients are responsible for this performance degradation. We now study this
apparent bias.

Figure 5:  We train the networks in Figure 2 on $\widetilde{\mathcal{L}}_{M}$
(4.5) for the sinusoidal toy-data with additive GNIs [left] and MNIST with
multiplicative GNIs [right]. We plot the training loss
($\mathcal{L}_{\mathrm{train}}$). Shading is the standard deviation over 5
random seeds.

Theoretical analysis of implicit bias. Due to their heavy-tailed nature,
stable processes have significantly different statistical properties from
those of their Brownian counterparts. Their trajectories have a countable
number of discontinuities, called _jumps_ , whereas Brownian motion is
continuous almost everywhere. With these jumps the process can escape from
‘narrow’ basins and spend more time in ‘wider’ basins (see Appendix F for the
definition of width). Theoretical results demonstrating this have been
provided for symmetric stable processes ($\theta=0$) (Şimşekli et al., 2019).
By translating the related metastability results from statistical physics
(Imkeller & Pavlyukevich, 2008) to our context, in Appendix F we illustrate
that this also holds for SDEs driven by asymmetric stable processes
($\theta\neq 0$). In this sense, the SDE (4.3) is ‘biased’ towards wider
basins.

While driving SGD iterates towards wider minima could be beneficial, we now
show that heavy tails can also introduce an undesirable bias and that this
bias is magnified by asymmetries. To quantify this bias, we focus on the
invariant measure (i.e., the stationary distribution) of the Markov process
(4.4) and investigate its modes (i.e., its local maxima), around where the
process resides most of the time.

In a statistical physics context, Dybiec et al. (2007) empirically illustrated
that the asymmetric stable noise can cause ‘shifts’, in the sense that the
modes of the stationary distribution of (4.4) can shift away from the true
local minima of $f$, which are of our interest as our aim is to minimise $f$.
They further illustrated that such shifts can be surprisingly large when
$\alpha$ gets smaller and $|\theta|$ gets larger. We illustrate this outcome
by reproducing one of the experiments provided in (Dybiec et al., 2007) in
Figure 6. Here, we consider a one-dimensional problem with the quartic
potential $f(w)=w^{4}/4-w^{2}/2$, and simulate (4.4) for $10$K iterations with
constant step-size $\eta_{k}=0.001$ and $\varepsilon=1$. By using the
generated iterates, we estimate the density of the invariant measure of (4.4)
by using the kernel density estimator provided in
$\mathrm{scikit}$-$\mathrm{learn}$ (Pedregosa et al., 2011), for different
values of $\alpha$ and $\theta$. When $\alpha$ is larger (left), the heavy-
tails cause a shift in the modes of the invariant measure, where these shifts
become slightly larger with increasing asymmetries ($|\theta|>0$). When the
tails are heavier (right), we observe a much stronger interaction between
$\alpha$ and $\theta$, and observe _drastic_ shifts as the asymmetry is
increased.

Figure 6: The stationary distributions of (4.4) with $\alpha=1.9$ (left) and
$\alpha=1.1$ (right). The solid black line represents the density of the Gibbs
measure $\exp(-f(w))$ with $f(w)=w^{4}/4-w^{2}/2$.

From an optimisation perspective, these results are rather unsettling as they
imply that SGD might spend most of its time in regions that are arbitrarily
far from the local minima of the objective function $f$, since the high
probability regions of its stationary distribution might be shifted away from
the local minima of interest. In the symmetric case ($\theta=0$), this
observation has been formally proven in (Sliusarenko et al., 2013) when $f$ is
chosen as the one dimensional quartic potential of Figure 6 and when
$\alpha=1$. A direct quantification of such shifts is non-trivial; and in the
presence of asymmetries ($\theta\neq 0$), even further difficulties emerge.

In a recent study Şimşekli et al. (2020) focused on eliminating the undesired
bias introduced by _symmetric_ stable noise in SGD with momentum (Qian, 1999),
and proposed an indirect way to ensure that modes of the stationary
distribution exactly match the objective function’s local minima. They
developed a ‘modified’ SDE whose invariant distribution can be proven to be
the _Gibbs measure_ , denoted by $\pi(d\boldsymbol{\mathbf{w}})$, which is a
probability measure that has a density proportional to
$\exp(-Cf(\boldsymbol{\mathbf{w}}))$ for some $C>0$. Clearly, all the local
maxima of this density coincide with the local minima of the function $f$;
hence, their approach eliminates the possibility of a shift in the modes by
imposing a stronger condition which controls the entire invariant
distribution.

By following a similar approach, we will bound the gap between the invariant
measure of (4.4) and the Gibbs measure in terms of the tail index $\alpha$ and
the skewness $\theta$, using it as a quantification of the bias induced by the
asymmetric heavy-tailed noise. In particular, for any sufficiently regular
test function $g$, we consider its expectation under the Gibbs measure
$\nu(g):=\int g(\boldsymbol{\mathbf{w}})\pi(d\boldsymbol{\mathbf{w}})$, and
its sample average computed over (4.4), i.e.,
${\nu}_{N}(g):=\frac{1}{H_{N}}\sum_{k=1}^{N}\eta_{k}g({\boldsymbol{\mathbf{w}}}_{k})$,
where $H_{N}=\sum_{k=1}^{N}\eta_{k}$. We then bound the _weak error_ :
$|\nu(g)-\lim_{N\to\infty}{\nu}_{N}(g)|$, whose convergence to zero is
sufficient for ensuring the modes do not shift.

We derive this bound for our case in three steps: (i) We first link the
discrete-time process (4.4) to its continuous-time limit (4.3) by directly
using the results of (Panloup, 2008). (ii) We then design a modified SDE that
has the unique invariant measure as the Gibbs measure, with all the modes
matching those of the loss function (Theorem 4.1). (iii) Finally, we show that
the SDE (4.3) is a poor numerical approximation to the modified SDE, and we
develop an upper-bound for the approximation error (Theorem 4.2).

To address (ii), we introduce a modification to (4.3), coined asymmetric
fractional Langevin dynamics:

$\displaystyle
d\boldsymbol{\mathbf{w}}_{t}=b(\boldsymbol{\mathbf{w}}_{t-},\alpha,\theta)dt+\varepsilon
d\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta}\,.$ (4.6)

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

$\displaystyle(b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}=\frac{\varepsilon^{\alpha}}{\varphi(\boldsymbol{\mathbf{w}})}\mathcal{D}^{\alpha-2,-\theta_{i}}_{w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})),$
(4.7)

where $\theta_{i}\in(-1,1)$, $1\leq i\leq d$, $1<\alpha<2$, and
$\varphi(\boldsymbol{\mathbf{w}}):=e^{-\varepsilon^{-\alpha}f(\boldsymbol{\mathbf{w}})}$.
Here, the operator $\mathcal{D}^{\alpha-2,-\theta_{i}}$ denotes a _Riesz-
Feller type fractional derivative_ (Gorenflo & Mainardi, 1998; Mainardi et
al., 2001) whose exact (and rather complicated) definition is not essential in
our problematic, and is given in Appendix E in order to avoid obscuring the
main results. The next theorem states that the SDE (4.6) targets the Gibbs
measure.

######  Theorem 4.1.

The Gibbs measure
$\pi(d\boldsymbol{\mathbf{w}})\propto\exp(-\varepsilon^{-\alpha}f(\boldsymbol{\mathbf{w}}))d\boldsymbol{\mathbf{w}}$
is an invariant distribution of (4.6). If
$b(\boldsymbol{\mathbf{w}},\alpha,\theta)$ is Lipschitz continuous in
$\boldsymbol{\mathbf{w}}$, then $\pi(d\boldsymbol{\mathbf{w}})$ is the unique
invariant distribution of (4.6).

This theorem states that the use of the modified drift $b$ in place of
$-\nabla f$, prevents any potential shifts in the modes of the invariant
measure.

The fractional derivative in (4.7) is a non-local operator that requires the
knowledge of the full function, and does not admit a closed-form expression.
In the next step, we develop an approximation scheme for the drift $b$ in
(4.6), and show that the gradient $-\nabla f$ in (4.4) appears as a special
case of this scheme. To simplify notation, we consider the one-dimensional
case ($d=1$); however, our results can be easily extended to multivariate
settings by applying the same approach to each coordinate. Hence, in the
general case the bounds will scale linearly with $d$. We define the following
approximation for $b$:

$\displaystyle
b_{h,K}(w,\alpha,\theta):=\frac{\varepsilon^{\alpha}}{\varphi(w)}\Delta_{h,K}^{\alpha-2,-\theta}(\partial_{w}\varphi(w)),$
(4.8)

where, for an arbitrary function $\psi$, we have

$\displaystyle\Delta_{h,K}^{\gamma,-\theta}\psi(w):=\frac{c_{\gamma}}{h^{\gamma}}\sum_{k=-K}^{K}(1+\theta\sgn(k))\tilde{g}_{\gamma,k}\psi(w-kh).$

Here, $\sgn$ denotes the sign function, $h>0$, $K\in\mathbb{N}\cup\\{0\\}$,
$c_{\gamma}:=1/(2\cos(\gamma\pi/2))$, and
$\tilde{g}_{\gamma,k}:=(-1)^{k}\Gamma(-\gamma+k)/\Gamma(k+1)\Gamma(-\gamma)$.
This approximation is designed in the way that we recover the original drift
$b$ as $h\to 0$ and $K\to\infty$ for sufficiently regular $\varphi$. It is
clear that when we set $K=0$ and
$h=h_{0}:=[2\varepsilon^{-\alpha}\cos((\alpha-2)\pi/2)]^{1/(2-\alpha)}$, we
have $b_{h,K}(w,\alpha,\theta)=-\partial_{w}f(w)$. In the multidimensional
case, where we apply this approximation to each coordinate, the same choice of
$K$ and $h$ gives us the original gradient $-\nabla f$; hence, we fall back to
the original recursion (4.4). By considering the recursion with this
approximate drift

$\tilde{\boldsymbol{\mathbf{w}}}_{n+1}=\tilde{\boldsymbol{\mathbf{w}}}_{n}+\eta_{n+1}b_{h,K}(\tilde{\boldsymbol{\mathbf{w}}}_{n},\alpha,\theta)+\varepsilon\eta_{n+1}^{1/\alpha}\Delta\boldsymbol{\mathbf{L}}^{\alpha,\theta}_{n+1},$

and the corresponding sample averages
$\tilde{\nu}_{N}(g):=\frac{1}{H_{N}}\sum_{k=1}^{N}\eta_{k}g(\tilde{\boldsymbol{\mathbf{w}}}_{k})$222Note
that, with the choice of $K=0$ and $h=h_{0}$, $\tilde{\nu}_{N}(g)$ reduces to
the original sample average $\nu_{N}(g)$., we are ready to state our error
bound. We believe this result is interesting on its own, and would be of
further interest in statistical physics and applied probability. To avoid
obscuring the result, we state the required assumptions in the Appendix, which
mainly require decreasing step-size and ergodicity.

######  Theorem 4.2.

Let $\gamma:=\alpha-2\in(-1,0)$. Suppose that the assumptions stated in the
Appendix hold. Then, the following bound holds almost surely:

$\displaystyle\left\lvert\nu(g)-\lim\nolimits_{N\rightarrow\infty}\tilde{\nu}_{N}(g)\right\rvert$
(4.9)
$\displaystyle\leq\frac{\tilde{C}}{4\pi(|\gamma|+2)}\left[|\theta||\gamma|+\left|\tan\left(\gamma\pi/2\right)\right||\gamma|\right]h$
$\displaystyle\qquad+\left((1+\theta)C^{\prime}_{0}+(1-\theta)C^{\prime\prime}_{0}\right)\frac{1}{hK}+\mathcal{O}(h^{2})\,,$

where $\tilde{C},\,C^{\prime}_{0},\,C^{\prime\prime}_{0}>0$ are constants.

We note our result extends the case $\alpha=2$, $\theta=0$ in Durmus &
Moulines (2017) and the case $\alpha\neq 2$, $\theta=0$ in Şimşekli (2017);
whereas we cover the case $\alpha\neq 2$, $\theta\in(-1,1)$. The right-hand-
side of (4.9) contains two main terms. The second term shows that the error
increases linearly with decreasing $K$, indicating that the error can be
_arbitrarily large_ when $K=0$, and the gap cannot be controlled without
imposing further assumptions on $f$.

More interestingly, even when $K$ goes to infinity (i.e., the second term
vanishes), the first term stays unaffected. Note that for large enough
$\varepsilon$, $h_{0}$ increases as $\alpha\in(1,2)$ decreases. In this
regime, the first term indicates that the error increases with decreasing
$\alpha$, and an additional error term appears whenever $\theta\neq 0$, which
is further amplified with the heaviness of the tails (measured by $|\gamma|$).
This outcome provides a theoretical justification to the empirical
observations stated in Figures 6 and 8.

Figure 7:  We show the test-set loss for SVHN [bottom] and CIFAR10 [top], for
2-layer convolutional (CONV) and 4-layer MLPs with 512 units per layer trained
with the explicit regulariser approximation $R(\cdot)$ of Camuto et al. (2020)
(Exp Reg), with additive-GNIs ($\sigma^{2}=0.1$) (Noise), and no
regularisation (Baseline).

###  4.1 Further Experiments

We have already ascertained that the bias implied by Theorem 4.2 has a visible
impact on training performance in Figures 4 and 5. We have also already shown
that the asymmetry and heavy-tails of the implicit-effect gradient noise are
responsible for this performance degradation in Figure 4: models trained with
Gaussian noise on gradients outperform models trained with
$\mathcal{S}_{\alpha}$ gradient noise, and those trained with the implicit
effect, _on training data without batching_. We corroborate these findings
with experiments with mini-batching and results on test data. In Figure 7 we
use the approximation of the explicit regulariser
$R(\mathcal{B};\boldsymbol{\mathbf{w}})$ derived by Camuto et al. (2020) for
computational efficiency. Convolutional networks trained with $R$ consistently
outperform those trained with GNIs and mini-batching on held-out data,
supporting that the implicit effect degrades performance. In Figure 8, we
sample $M$ multiplicative-GNI samples and marginalise out the implicit effect
as before. We model the gradients of the implicit effect, $\nabla
E_{\mathcal{L}}(\cdot)$, as an $\mathcal{S}_{\alpha}$ distribution.
Empirically, we found that when increasing the variance ($\sigma^{2}$) of the
injected noise, the gradient noise $\nabla E_{\mathcal{L}}(\cdot)$ becomes
increasingly heavy-tailed and skewed, i.e. $\alpha$ decreases and $|\theta|$
increases, and in tandem larger $M$ models begin to outperform smaller $M$
models on _held-out data_ , when trained with mini-batches. These results
support that GNIs induce bias in SGD because of the asymmetric heavy-tailed
noise they induce on gradient updates.

Figure 8: [first row] We train 2-dense-layer-256-unit-per-layer ELU networks
on the objective
$\frac{1}{M}\sum_{m=0}^{M}\widetilde{\mathcal{L}}(\mathcal{B};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}}_{m})$
with a cross-entropy loss (see Appendix B) for SVHN. We use multiplicative
noise of variance $\sigma^{2}$ and batch size of 512. We plot the test-set
loss ($\mathcal{L}_{\mathrm{test}}$). [second row] We fit univariate
$\mathcal{S}_{\alpha}$ via maximum likelihood (Nolan, 2001) to $\nabla
E_{\mathcal{L}}(\cdot)$ and show KDE plots of parameters’ estimates.

##  5 Conclusion

Our work lays the foundations for the study of regularisation methods from the
perspective of SDEs. We have shown that Gaussian Noise Injections (GNIs),
though they inject Gaussian noise in the forward pass, induce asymmetric
heavy-tailed noise on gradient updates by way of the implicit effect. By
modelling the overall induced noise using an asymmetric $\alpha$-stable noise,
we demonstrate that the stationary distribution of this process gets
arbitrarily distant from the so-called Gibbs measure, whose modes exactly
match the local minima of the loss function, shedding light on why neural
networks trained with GNIs underperform networks trained solely with the
explicit effect. Given the deleterious effects of asymmetric gradient noise on
gradient descent, extensions of this work could focus on methods that
symmetrise gradient noise, stemming from batching or noise injections, so as
to limit these negative effects.

## Acknowledgements

This research was directly funded by the Alan Turing Institute under
Engineering and Physical Sciences Research Council (EPSRC) grant EP/N510129/1.
Alexander Camuto was supported by an EPSRC Studentship. Xiaoyu Wang and
Lingjiong Zhu are partially supported by the grant NSF DMS-2053454 from the
National Science Foundation. Lingjiong Zhu is also grateful to the partial
support from a Simons Foundation Collaboration Grant. Mert Gürbüzbalaban’s
research is supported in part by the grants Office of Naval Research Award
Number N00014-21-1-2244, National Science Foundation (NSF) CCF-1814888, NSF
DMS-1723085, NSF DMS-2053485. Umut Şimşekli’s research is partly supported by
the French government under management of Agence Nationale de la Recherche as
part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001
(PRAIRIE 3IA Institute).

## References

  * Agapiou et al. (2014) Agapiou, S., Stuart, A. M., and Zhang, Y. X.  Bayesian posterior contraction rates for linear severely ill-posed inverse problems.  _Journal of Inverse and Ill-Posed Problems_ , 22(3):297–321, 2014.
  * Bishop (1995) Bishop, C. M.  Training with Noise is Equivalent to Tikhonov Regularization.  _Neural Computation_ , 7(1):108–116, 1995.
  * Camuto et al. (2020) Camuto, A., Willetts, M., Şimşekli, U., Roberts, S., and Holmes, C.  Explicit Regularisation in Gaussian Noise Injections.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2020.
  * Cohen et al. (2019) Cohen, J., Rosenfeld, E., and Kolter, J. Z.  Certified adversarial robustness via randomized smoothing.  _36th International Conference on Machine Learning, ICML 2019_ , 2019-June:2323–2356, 2019.
  * Dieng et al. (2018) Dieng, A. B., Ranganath, R., Altosaar, J., and Blei, D. M.  Noisin: Unbiased regularization for recurrent neural networks.  In _Proceedings of the 35th International Conference on Machine Learning_ , pp. 1252–1261, 2018.
  * Duan (2015) Duan, J.  _An Introduction to Stochastic Dynamics_.  Cambridge University Press, New York, 2015.
  * Durmus & Moulines (2017) Durmus, A. and Moulines, E.  Nonasymptotic convergence analysis for the unadjusted Langevin algorithm.  _Annals of Applied Probability_ , 27(3):1551–1587, 2017.
  * Dybiec et al. (2007) Dybiec, B., Gudowska-Nowak, E., and Sokolov, I.  Stationary states in Langevin dynamics under asymmetric Lévy noises.  _Physical Review E_ , 76(4):041122, 2007.
  * Fiche et al. (2013) Fiche, A., Cexus, J. C., Martin, A., and Khenchaf, A.  Features modeling with an $\alpha$-stable distribution: Application to pattern recognition based on continuous belief functions.  _Information Fusion_ , 14(4):504–520, 2013.
  * Gao et al. (2018) Gao, X., Gürbüzbalaban, M., and Zhu, L.  Global Convergence of Stochastic Gradient Hamiltonian Monte Carlo for Non-Convex Stochastic Optimization: Non-Asymptotic Performance Bounds and Momentum-Based Acceleration.  _arXiv:1809.04618_ , 2018.
  * Gao et al. (2020) Gao, X., Gürbüzbalaban, M., and Zhu, L.  Breaking reversibility accelerates Langevin dynamics for global non-convex optimization.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2020.
  * Gnedenko & Kolmogorov (1954) Gnedenko, B. V. and Kolmogorov, A.  _Limit Distributions for Sums of Independent Random Variables_.  Addison-Wesley, Cambridge, MA, 1954.  Translated by Kai Lai Chung.
  * Gorenflo & Mainardi (1998) Gorenflo, R. and Mainardi, F.  Random walk models for space-fractional diffusion processes.  _Fractional Calculus & Applied Analysis_, 1:167–191, 1998\.
  * Gürbüzbalaban & Hu (2020) Gürbüzbalaban, M. and Hu, Y.  Fractional moment-preserving initialization schemes for training fully-connected neural networks.  _arXiv preprint arXiv:2005.11878_ , 2020.
  * Gürbüzbalaban et al. (2020) Gürbüzbalaban, M., Şimşekli, U., and Zhu, L.  The heavy-tail phenomenon in SGD.  _arXiv preprint arXiv:2006.04740_ , 2020.
  * Hodgkinson & Mahoney (2020) Hodgkinson, L. and Mahoney, M. W.  Multiplicative noise and heavy tails in stochastic optimization.  _arXiv preprint arXiv:2006.06293_ , 2020.
  * Imkeller & Pavlyukevich (2008) Imkeller, P. and Pavlyukevich, I.  Metastable behavior of small noise Lévy-driven diffusions.  _ESAIM: Probability and Statistics_ , 12:412–437, 2008\.
  * Jastrzȩbski et al. (2017) Jastrzȩbski, S., Kenton, Z., Arpit, D., Ballas, N., Fischer, A., Bengio, Y., and Storkey, A.  Three Factors Influencing Minima in SGD.  _arXiv:1711.04623_ , 2017.
  * Kingma et al. (2015) Kingma, D. P., Salimans, T., and Welling, M.  Variational dropout and the local reparameterization trick.  In _Advances in Neural Information Processing Systems_ , volume 2015-Janua, pp. 2575–2583, 2015.
  * Kroneburg (2011) Kroneburg, M.  The binomial coefficient for negative arguments.  _arXiv preprint arXiv:1105.3689_ , 2011.
  * Kuchibhotla & Chakrabortty (2018) Kuchibhotla, A. K. and Chakrabortty, A.  Moving beyond sub-Gaussianity in high dimensional statistics: Applications in covariance estimation and linear regression.  _arXiv:1804.02605_ , 2018.
  * Li et al. (2017) Li, Q., Tai, C., and E, W.  Stochastic modified equations and adaptive stochastic gradient algorithms.  In _Proceedings of the 34th International Conference on Machine Learning_ , pp. 2101–2110, 06–11 Aug 2017.
  * Mainardi et al. (2001) Mainardi, F., Luchko, Y., and Pagnini, G.  The fundamental solution of the space-time fractional diffusion equation.  _Fractional Calculus & Applied Analysis_, 4:153–192, 2001\.
  * Mandt et al. (2016) Mandt, S., Hoffman, M. D., and Blei, D. M.  A variational analysis of stochastic gradient algorithms.  _33rd International Conference on Machine Learning, ICML 2016_ , 1:555–566, 2016.
  * Meerschaert & Tadjeran (2004) Meerschaert, M. and Tadjeran, C.  Finite difference approximation for fractional advection-dispersion flow equations.  _Journal of Computational and Applied Mathematics_ , 172:65–77, 2004.
  * Nadarajah & Pogány (2016) Nadarajah, S. and Pogány, T. K.  On the distribution of the product of correlated normal random variables.  _Comptes Rendus Mathematique_ , 354(2):201–204, 2016.
  * Nolan (2001) Nolan, J. P.  Maximum likelihood estimation and diagnostics for stable distributions.  In Barndorff-Nielsen, O. E., Resnick, S. I., and Mikosch, T. (eds.), _Lévy Processes: Theory and Applications_ , pp. 379–400. Birkhäuser Boston, Boston, MA, 2001.
  * Oliveira et al. (2016) Oliveira, A., Oliveira, T. A., and Seijas-Macias, A.  Skewness into the product of two normally distributed variables and the risk consequences.  _Revstat Statistical Journal_ , 14(2):119–138, 2016.
  * Ortigueira (2006) Ortigueira, M. D.  Riesz potential operators and inverses via fractional centred derivatives.  _International Journal of Mathematics and Mathematical Sciences_ , 2006(Article ID 48391):1–12, 2006.
  * Ortiguera (2006b) Ortiguera, M. D.  Fractional central differences and derivatives.  _IFAC Proceedings Volumes_ , 39(11):58–63, 2006b.
  * Panloup (2008) Panloup, F.  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.
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Poole et al. (2014) Poole, B., Sohl-Dickstein, J., and Ganguli, S.  Analyzing noise in autoencoders and deep networks.  _arXiv:1406.1831_ , 2014.
  * Qian (1999) Qian, N.  On the momentum term in gradient descent learning algorithms.  _Neural networks_ , 12(1):145–151, 1999.
  * Raginsky et al. (2017) Raginsky, M., Rakhlin, A., and Telgarsky, M.  Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis.  In _Conference on Learning Theory_ , pp. 1674–1703, 2017.
  * Robbins & Monro (1951) Robbins, H. and Monro, S.  A Stochastic Approximation Method.  _The Annals of Mathematical Statistics_ , 22(3):400–407, 1951.
  * Roberts & Stramer (2002) Roberts, G. O. and Stramer, O.  Langevin diffusions and Metropolis-Hastings algorithms.  _Methodology and Computing in Applied Probability_ , 4(4):337–357, 2002.
  * Ruder (2016) Ruder, S.  An overview of gradient descent optimization algorithms.  _arXiv:1609.04747_ , 2016.
  * Samorodnitsky & Taqqu (1994) Samorodnitsky, G. and Taqqu, M. S.  _Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance_.  Chapman & Hall, New York, 1994.
  * Sarafrazi & Yazdi (2019) Sarafrazi, K. and Yazdi, M.  Skewed alpha-stable distribution for natural texture modeling and segmentation in contourlet domain.  _Eurasip Journal on Image and Video Processing_ , 2019(1):1–12, 2019.
  * Schertzer et al. (2001) Schertzer, D., Larchevêque, M., Duan, J., Yanovsky, V., and Lovejoy, S.  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 (2017) Şimşekli, U.  Fractional Langevin Monte Carlo: Exploring Lévy driven stochastic differential equations for Markov Chain Monte Carlo.  In _International Conference on Machine Learning_ , pp. 3200–3209, 2017.
  * Şimşekli et al. (2019) Şimşekli, U., Sagun, L., and Gürbüzbalaban, M.  A tail-index analysis of stochastic gradient noise in deep neural networks.  In _Proceedings of the 36th International Conference on Machine Learning_ , pp. 5827–5837, 2019.
  * Şimşekli et al. (2020) Şimşekli, U., Zhu, L., Teh, Y. W., and Gürbüzbalaban, M.  Fractional Underdamped Langevin Dynamics: Retargeting SGD with Momentum under Heavy-Tailed Gradient Noise.  In _Proceedings of the 37th International Conference on Machine Learning_ , pp. 8970–8980, 2020.
  * Şimşekli et al. (2019) Şimşekli, U., Gürbüzbalaban, M., Nguyen, T. H., Richard, G., and Sagun, L.  On the Heavy-Tailed Theory of Stochastic Gradient Descent for Deep Neural Networks.  _arXiv:1912.00018_ , 2019.
  * Sliusarenko et al. (2013) Sliusarenko, O. Y., Surkov, D., Gonchar, V. Y., and Chechkin, A. V.  Stationary states in bistable system driven by Lévy noise.  _The European Physical Journal Special Topics_ , 216(1):133–138, 2013.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15:1929–1958, 2014\.
  * Tian et al. (2015) Tian, W., Zhou, H., and Deng, W.  A class of second order difference approximations for solving space fractional diffusion equations.  _Mathematics of Computation_ , 84(294):1703–1727, 2015.
  * Vershynin (2018) Vershynin, R.  _High-Dimensional Probability: An Introduction with Applications in Data Science_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2018.
  * Vladimirova et al. (2019) Vladimirova, M., Verbeek, J., Mesejo, P., and Arbel, J.  Understanding priors in Bayesian neural networks at the unit level.  _36th International Conference on Machine Learning, ICML 2019_ , 2019-June:11248–11257, 2019.
  * Vladimirova et al. (2020) Vladimirova, M., Girard, S., Nguyen, H., and Arbel, J.  Sub‐Weibull distributions: generalizing sub‐Gaussian and sub‐Exponential properties to heavier‐tailed distributions.  _Stat_ , 9(1):1–10, 2020.
  * Webb (1994) Webb, A. R.  Functional Approximation by Feed-Forward Networks: A Least-Squares Approach to Generalization.  _IEEE Transactions on Neural Networks_ , 5(3):363–371, 1994.
  * Wei et al. (2020) Wei, C., Kakade, S., and Ma, T.  The Implicit and Explicit Regularization Effects of Dropout.  In _Proceedings of the 37th International Conference on Machine Learning_ , pp. 10181–10192, 2020.
  * Welling & Teh (2011) Welling, M. and Teh, Y. W.  Bayesian Learning via Stochastic Gradient Langevin Dynamics.  In _Proceedings of the 28th International Conference on Machine Learning_ , ICML’11, pp. 681–688, Madison, WI, USA, 2011. Omnipress.
  * Zhang et al. (2017) Zhang, C., Recht, B., Bengio, S., Hardt, M., and Vinyals, O.  Understanding deep learning requires rethinking generalization.  In _5th International Conference on Learning Representations, ICLR 2017_ , 2017.
  * Zhang et al. (2020) Zhang, J., Karimireddy, S. P., Veit, A., Kim, S., Reddi, S., Kumar, S., and Sra, S.  Why are adaptive methods good for attention models?  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 33, 2020.
  * Zhou et al. (2020) Zhou, P., Feng, J., Ma, C., Xiong, C., Hoi, S., and E, W.  Towards theoretically understanding why SGD generalizes better than ADAM in deep learning.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , volume 33, 2020.

The supplementary document is organised as follows.

  1. 1.

The supplementary document begins first with a presentation of additional
experiments that are referenced directly in the main text (Section A).

  2. 2.

We then cover the cost-functions used to train neural networks in Section B;
and give an overview in Section C of the other potential sources of the
implicit effect gradient noise skew which we explored.

  3. 3.

In Section D, we provide an overview of the assumptions we will be making in
our analysis. We then describe in Section E the numerical method we use to
approximate the drift term $b(\boldsymbol{\mathbf{w}},\alpha,\theta)$ defined
in (4.7).

  4. 4.

We end with metastability analysis of asymmetric stable processes (Section F);
followed by the technical proofs of the lemmas, theorems, and corollaries that
we present in the main body and the supplementary document of the paper
(Section G).

Before beginning the supplementary document we make a quick note of network
architectures and training hyper-parameters.

#### Network Architectures

Networks were trained using stochastic gradient descent with a learning rate
of 0.0003 and batch sizes specified in text. MLP network architectures are
specified in text. Convolutional (CONV) networks are 2 hidden layer networks.
The first layer has 32 filters, a kernel size of 4, and a stride length of 2.
The second layer has 128 filters, a kernel size of 4, and a stride length of
2. The final output layer is a dense layer.

##  Appendix A Additional Experimental Results

Figure A.1: We measure the skewness and kurtosis at initialisation of the
noise accumulated on network activations at each layer $i$ for a 4-layer
256-unit wide MLP trained to regress the function
$\lambda(\boldsymbol{\mathbf{x}})=\sum_{i}\sin(2\pi
q_{i}\boldsymbol{\mathbf{x}}+\phi(i))$ with
$q_{i}\in(5,10,\dots,45,50),~{}\boldsymbol{\mathbf{x}}\in\mathbb{R}$ and
experiencing additive-GNIs. We plot the probability density function of
positive samples, comparing against half-normal (non-heavy-tailed) and half-
Cauchy (heavy-tailed) distributions, where
$\mathcal{E}_{i}(\boldsymbol{\mathbf{x}};\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})$
is defined in (2.9). Each blue point represents the noise on an individual
activation in a layer $i$ for a point $\boldsymbol{\mathbf{x}}$. This noise is
Gaussian (low skew and kurtosis) with a p.d.f. that tracks that of a half-
normal.  Figure A.2:  Here we show the same plots as in Figure A.1 but for
multiplicative-GNIs. The forward pass here experiences symmetric heavy-tailed
noise for all layers past the data layer.  Figure A.3:  Here we show the same
plots as in Figure 2 but for multiplicative-GNIs. The gradient noise is skewed
and heavy-tailed, with a p.d.f. that is more Cauchy-like than Gaussian. The
kurtosis decays as the gradients approach the input layer, as predicted by
Theorem 3.1.

(a) CIFAR10 ($\times$) (b) SVHN ($\times$) (c) BHP ($+$)

Figure A.4: We model $\nabla E_{\mathcal{L}}(\cdot)$ as being drawn from some
$\alpha$-stable distribution $\mathcal{S}_{\alpha}$ and estimate the tail-
index $\alpha$ and skewness $\theta$ using maximum likelihood estimation as in
Nolan (2001). We plot the results as a scatter for a batch of size $B=512$ for
CIFAR10 and SVHN and $B=32$ for Boston House Prices. Additive ($+$) and
multiplicative ($\times$) GNIs have $\sigma^{2}=0.1$.

##  Appendix B Cost Functions

###  B.1 Mean Square Error

In the case of regression the most commonly used loss is the mean-square
error.

$\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})=\left(\boldsymbol{\mathbf{y}}-\boldsymbol{\mathbf{h}}_{L}(\boldsymbol{\mathbf{x}})\right)^{2}\,.$

###  B.2 Cross Entropy Loss

In the case of classification, we use the cross-entropy loss. If we consider
our network outputs $\boldsymbol{\mathbf{h}}_{L}$ to be the
pre-$\mathrm{softmax}$ of logits of the final layer then the loss is for a
data-label pair $(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$

$\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})=-\sum_{c=0}^{C}\boldsymbol{\mathbf{y}}_{c}\log\left(\mathrm{softmax}(\boldsymbol{\mathbf{h}}_{L}(\boldsymbol{\mathbf{x}}))_{c}\right)\,,$
(B.1)

where $c$ indexes over the $C$ possible classes of the classification problem.

##  Appendix C Other Potential Sources of the Skewness in the Gradient Noise

The product of correlated random variables can be skewed (Oliveira et al.,
2016; Nadarajah & Pogány, 2016). Our first hypothesis was that the skew came
from the correlation of $(\partial E_{\mathcal{L}}(\cdot)/\partial h_{i}^{m})$
and $(\partial h_{i}^{m}/\partial W_{i,l,j})$. As a test, Figure C.5 of the
Appendix reproduces Figure 2 with linear $\kappa$, isolating gradient
correlation as a potential source of skew. Gradients are not skewed,
demonstrating that the asymmetry stems from non-linear $\kappa$.

Figure C.5: Here we reproduce Figure 2 but with no non-linearities. This
gradient noise is clearly heavy-tailed but not skewed.

##  Appendix D Overview of the Assumptions

Due to space limitations and to avoid obscuring the main take home messages of
our theoretical results, we did not present the two assumption required for
Theorem 4.2. These two assumptions are properly presented in their respective
sections (Sections E and G.6), where we first provide the required technical
context for defining them in each section. In this section, we will shortly
discuss the semantics of these assumptions from a higher-level perspective for
the convenience of the reader.

  * •

Assumption E.1. This assumption is essentially an assumption of the tails of
the function $\partial\varphi$, with
$\varphi({w})=e^{-\varepsilon^{-\alpha}f({w})}$. In particular, in order to
make our approximation scheme (to the fractional derivatives) convergent, this
assumption makes sure that outside of a compact region, the function
$\partial\varphi$ exponentially decays.

  * •

Assumption G.1. This assumption enforces a certain structure on the Euler-
Maruyama discretisation given in Section 4:

$\tilde{\boldsymbol{\mathbf{w}}}_{n+1}=\tilde{\boldsymbol{\mathbf{w}}}_{n}+\eta_{n+1}b_{h,K}(\tilde{\boldsymbol{\mathbf{w}}}_{n},\alpha,\theta)+\varepsilon\eta_{n+1}^{1/\alpha}\Delta\boldsymbol{\mathbf{L}}^{\alpha,\theta}_{n+1}.$

As a first condition, we make sure that the step-sizes are decreasing while
their sum is diverging, which is a standard assumption. The second condition
is essentially a Lyapunov condition that requires the modified drift $b$
behaves well, so that we can control the weak error of the sample averages by
using (Panloup, 2008). The final condition is similar to the second condition
in nature, and requires ergodicity of an SDE defined through the approximate
drift $b_{h,K}$, in order to enable us link the weak error to the error
induced by the approximation scheme used for the fractional derivatives.

##  Appendix E Fractional Differentiation and the Approximation Scheme

In this section, we provide the details of the Riesz-Feller type fractional
derivative $\mathcal{D}^{\alpha-2,-\theta}$, whose definition was omitted in
the main document for clarity. We then present the details of the
approximation method for the drift term
$b(\boldsymbol{\mathbf{w}},\alpha,\theta)$ defined in (4.7).

The building block of our analysis is a first-order approximation of
$\mathcal{D}^{\alpha-2,-\theta}\partial_{w}\varphi$ for any
$\partial_{w}\varphi\in L^{1}(\mathbb{R})\cap\mathcal{C}^{4}(\mathbb{R})$. We
consider the one-dimensional case for simplicity since the Lévy motion we
consider has independent coordinates, and the multi-dimensional numerical
approximation can be reduced to the one-dimensional case. Assume the tail
index $1<\alpha<2$ and the skewness parameter satisfies $-1<\theta<1$.

When $\theta=0$, Şimşekli (2017) developed the numerical approximation method
for the drift term $b(\boldsymbol{\mathbf{w}},\alpha,0)$ by approximately
computing the Riesz potential333Note that when $\theta=0$,
$\mathcal{D}^{\alpha-2}=\mathcal{D}^{\alpha-2,0}$ recovers the Riesz
potential. $\mathcal{D}^{\gamma}$ via the fractional centred difference method
provided by Ortigueira (2006); Ortiguera (2006b). It is shown that for any
$-1<\gamma<0$, we have the following numerical error,

$\displaystyle\left|\mathcal{D}^{\gamma}\partial_{w}\varphi(w)-\Delta_{h,K}^{\gamma}\partial_{w}\varphi(w)\right|=\mathcal{O}\left(h^{2}+1/(hK)\right)\,,$

as $h\rightarrow 0$, where $K\in\mathbb{N}\cup\\{0\\}$ is the truncation
parameter and $\partial_{w}\varphi(w)$ satisfies some regularity conditions,
and the operator $\Delta_{h,K}^{\gamma}$ is given by

$\displaystyle\Delta_{h,K}^{\gamma}f(w)=\frac{1}{h^{\gamma}}\sum_{k=-K}^{K}g_{\gamma,k}f(w-kh),$

for any test function $f$ satisfying some regularity conditions, where

$g_{\gamma,k}:=\frac{(-1)^{k}\Gamma(\gamma+1)}{\Gamma\left(\frac{\gamma}{2}-k+1\right)\Gamma\left(\frac{\gamma}{2}+k+1\right)}.$

We study the numerical error when approximating the drift term
$b(\boldsymbol{\mathbf{w}},\alpha,\theta)$ with the skewness parameter
$-1<\theta<1$ and provide the truncation error with a truncation parameter $K$
in Corollary E.1. Based on this result, Theorem 4.2 follows, which quantifies
the bias induced by $\alpha$-stable noise on gradient updates using the Euler-
Maruyama scheme.

Instead of using the centred difference method to implement the approximation
for the $\theta=0$ case, we tackle the more general $\theta\neq 0$ case by
using shifted Grünwald-Letnikov difference operators to approach the left and
right fractional derivative respectively. Let us define the parameter
$-1<\gamma:=\alpha-2<0$. Then, we can now formally define the Riesz-Feller
type fractional derivative operator as follows:

$\displaystyle\mathcal{D}^{\gamma,-\theta}f(w):=\frac{1}{2\cos(\gamma\pi/2)}\left[(1-\theta)\mathcal{I}^{-\gamma}_{+}f(w)+(1+\theta)\mathcal{I}^{-\gamma}_{-}f(w)\right]\,,$
(E.1)

with

$\mathcal{I}^{-\gamma}_{\pm}f(w):=\frac{1}{\Gamma(-\gamma)}\int_{0}^{\infty}\frac{f(w\pm\xi)}{\xi^{\gamma+1}}d\xi\,.$
(E.2)

Before we proceed, we first introduce difference operators
$\mathcal{A}_{h,p}^{\gamma}$ and $\mathcal{B}_{h,q}^{\gamma}$, where $p$ and
$q$ are two non-negative integers chosen to be the shifted parameters,

$\displaystyle\mathcal{A}^{\gamma}_{h,p}f(w)=\frac{1}{h^{\gamma}}\sum_{k=0}^{\infty}\tilde{g}_{\gamma,k}f(w-(k-p)h),$
(E.3)
$\displaystyle\mathcal{B}_{h,q}^{\gamma}f(w)=\frac{1}{h^{\gamma}}\sum_{k=0}^{\infty}\tilde{g}_{\gamma,k}f\left(w+(k-q)h\right).$
(E.4)

Essentially, we defined a forward shifted difference operator
$\mathcal{A}_{h,p}^{\gamma}$ to approximate the left fractional derivatives,
and a backward shifted difference operator $\mathcal{B}_{h,p}^{\gamma}$ to
approximate the right one. The coefficients
$\tilde{g}_{\gamma,k}:=\frac{(-1)^{k}\Gamma(-\gamma+k)}{\Gamma(k+1)\Gamma(-\gamma)}$
are from the coefficients of of the power series $(1-z)^{\gamma}$ with
$-1<\gamma<0$ and $|z|\leq 1$. For any negative fractional number
$-1<\gamma<0$ and $|z|\leq 1$, we have

$\displaystyle(1-z)^{\gamma}=\sum_{k=0}^{\infty}(-1)^{k}\binom{-\gamma+k-1}{k}z^{k},$
(E.5)

where the binomial coefficient $\binom{-\gamma+k-1}{k}$ is well-defined and
the binomial series converges for any complex number $\lvert z\rvert\leq 1$;
see e.g. Kroneburg (2011). Indeed, when $-1<\gamma<0$, we get

$\displaystyle\binom{-\gamma+k-1}{k}=\frac{\Gamma(-\gamma+k)}{\Gamma(k+1)\Gamma(-\gamma)}.$
(E.6)

We first present the following first-order approximation result of the
fractional derivative $\mathcal{D}^{\gamma,-\theta}$.

######  Theorem E.1.

Let $\mathcal{D}^{\gamma,-\theta}$ denote the fractional derivative for
$-1<\gamma<0$ and $-1<\theta<1$ as in (E.1). Suppose the function $f\in
L^{1}(\mathbb{R})\cap\mathcal{C}^{4}(\mathbb{R})$. Define

$\displaystyle\Delta_{h,p,q}^{\gamma,-\theta}f(w)=\frac{1}{2\cos(\gamma\pi/2)}\left[(1+\theta)\mathcal{A}_{h,p}^{\gamma}f(w)+(1-\theta)\mathcal{B}_{h,q}^{\gamma}f(w)\right].$
(E.7)

Then $\Delta_{h,p,q}^{\gamma,-\theta}f(w)$ is an approximation of
$\mathcal{D}^{\gamma,-\theta}f(w)$ with the first-order accuracy:

$\displaystyle\left\lvert\mathcal{D}^{\gamma,-\theta}f(w)-\Delta_{h,p,q}^{\gamma,-\theta}f(w)\right\rvert$
$\displaystyle\qquad\leq\left[|p-q|+|\theta|(p+q-\gamma)+\left|\tan\left(\frac{\gamma\pi}{2}\right)\right|\left(p+q-\gamma+|\theta||p-q|\right)\right]\,\frac{C}{4\pi(|\gamma|+2)}h+\mathcal{O}\left(h^{2}\right),$
(E.8)

as $h\rightarrow 0$, uniformly for all $w\in\mathbb{R}$, where $C>0$ is a
constant that may depend on $f$ and $\mathcal{O}(\cdot)$ hides the dependence
on $p$, $q$ and $\gamma$.

Next, we provide an error bound for numerically computing the drift term
$b(\boldsymbol{\mathbf{w}},\alpha,\theta)$ by truncating the approximation
series in Theorem E.1 as follows. Let us first define the operators
$\mathcal{A}^{\gamma}_{h,p,K}$ and $\mathcal{B}^{\gamma}_{h,q,K}$:

$\displaystyle\mathcal{A}^{\gamma}_{h,p,K}f(w):=\frac{1}{h^{\gamma}}\sum_{k=0}^{K}\tilde{g}_{\gamma,k}f\left(w-(k-p)h\right),$
(E.9)
$\displaystyle\mathcal{B}_{h,q,K}^{\gamma}f(w):=\frac{1}{h^{\gamma}}\sum_{k=0}^{K}\tilde{g}_{\gamma,k}f\left(w+(k-q)h\right),$
(E.10)

with
$\tilde{g}_{\gamma,k}:=\frac{(-1)^{k}\Gamma(-\gamma+k)}{\Gamma(k+1)\Gamma(-\gamma)}$,
and $K\in\mathbb{N}\cup\\{0\\}$.

Before we state the next result, let us first introduce the following
assumption.

######  Assumption E.1.

Suppose the function $\partial_{w}\varphi\in
L^{1}(\mathbb{R})\cap\mathcal{C}^{4}(\mathbb{R})$. In addition, there exist
constants $C_{p},C_{q}>0$ satisfying

$\lvert\partial_{w}\varphi(w-|k-p|h)\rvert\leq
C_{p}e^{-|k-p|h},\quad\lvert\partial_{w}\varphi(w+|k-q|h)\rvert\leq
C_{q}e^{-|k-q|h},$ (E.11)

and $\min\left\\{|k-p|,\,|k-q|\right\\}>K$ for the constant
$K\in\mathbb{N}\cup\\{0\\}$.

We have the following result.

######  Corollary E.1.

Suppose Assumption E.1 holds for $\partial_{w}\varphi$, and recall the
truncated series $\mathcal{A}_{h,p,K}^{\gamma}$ and
$\mathcal{B}_{h,q,K}^{\gamma}$ with $K\in\mathbb{N}\cup\\{0\\}$ defined in
(E.9) and (E.10). Let us also define the operator:

$\Delta_{h,p,q,K}^{\gamma,-\theta}=\frac{1}{2\cos(\gamma\pi/2)}\left[(1-\theta)\mathcal{B}_{h,q,K}^{\gamma}+(1+\theta)\mathcal{A}_{h,p,K}^{\gamma}\right].$
(E.12)

Then the truncation error is bounded in first-order accuracy as follows,

$\displaystyle\left|\mathcal{D}^{\gamma,-\theta}\partial_{w}\varphi(w)-\Delta_{h,p,q,K}^{\gamma,-\theta}\partial_{w}\varphi(w)\right|$
$\displaystyle\qquad\leq\frac{C}{4\pi(|\gamma|+2)}\left[|p-q|+|\theta|(p+q-\gamma)+\left|\tan\left(\frac{\gamma\pi}{2}\right)\right|\left(p+q-\gamma+|\theta||p-q|\right)\right]h$
$\displaystyle\qquad\qquad\qquad+\left((1+\theta)C_{p}+(1-\theta)C_{q}\right)\frac{1}{hK}+\mathcal{O}(h^{2}),$
(E.13)

where $C,C_{p},C_{q}>0$ are constants that may depend on $\partial_{w}\varphi$
and $\mathcal{O}(\cdot)$ hides the dependence on $p$ and $q$.

In particular, by taking $p=q=0$, Corollary E.1 implies that

$\displaystyle\left|\mathcal{D}^{\gamma,-\theta}\partial_{w}\varphi(w)-\Delta_{h,p=0,q=0,K}^{\gamma,-\theta}\partial_{w}\varphi(w)\right|$
$\displaystyle\qquad\leq\frac{C}{4\pi}\left(|\theta|+\left|\tan\left(\frac{\gamma\pi}{2}\right)\right|\right)h+\left((1+\theta)C_{p=0}+(1-\theta)C_{q=0}\right)\frac{1}{hK}+\mathcal{O}(h^{2}),$
(E.14)

where $\Delta_{h,K}^{\gamma,-\theta}=\Delta_{h,p=0,q=0,K}^{\gamma,-\theta}$.

Corollary E.1 implies that one can approximate $\mathcal{D}^{\gamma,-\theta}$
by the truncated $\Delta_{h,K}^{\gamma,-\theta}$ instead of
$\Delta_{h}^{\gamma,-\theta}$. Based on this result, we are able to quantify
in Theorem 4.2 the bias induced when implementing Euler-Maruyama scheme to
approximate the expectation of a test function $g$ with respect to the target
distribution $\pi$, where $\nu(g)=\int
g(\boldsymbol{\mathbf{w}})\pi(d\boldsymbol{\mathbf{w}})$.

##  Appendix F Metastability Analysis

In this section, we will focus on the metastability properties of the process

$\displaystyle dw_{t}=-\nabla_{w}f(w_{t})dt+\varepsilon
dL_{t}^{\alpha,\theta}.$ (F.1)

We will be interested in the first exit time, which is, roughly speaking, the
expected time required for the process to exit a neighborhood of a local
minimum. We will summarise the related theoretical results, which show that
the first exit time behaviour of systems driven by asymmetric stable processes
are similar to the ones of driven by symmetric stable processes. This
informally implies that the process will quickly escape from narrow minima
regions and will spend more time (in fact will get stuck) in wide minima
regions. In this section, we make this argument rigorous.

For simplicity of the presentation, we consider the one-dimensional case where
$L_{t}^{\alpha,\theta}$ is an asymmetric $\alpha$-stable Lévy process with
Lévy measure

$\displaystyle\nu(dy)=\left(\frac{1-\theta}{2}c_{\alpha}1_{y<0}+\frac{1+\theta}{2}c_{\alpha}1_{y>0}\right)\frac{dy}{|y|^{1+\alpha}},$
(F.2)

where $\theta\in(-1,1)$ and $\alpha\in(0,2)$ and
$c_{\alpha}:=\frac{\alpha}{\Gamma(1-\alpha)\cos(\pi\alpha/2)}$. Then, the left
and right tails of the Lévy measure are given by

$\displaystyle
H_{-}(-u):=\int_{(-\infty,-u)}\nu(dy)=\frac{1-\theta}{2}C_{\alpha}u^{-\alpha},$
$\displaystyle
H_{+}(u):=\int_{(u,+\infty)}\nu(dy)=\frac{1+\theta}{2}C_{\alpha}u^{-\alpha},$

where $C_{\alpha}:=\frac{1-\alpha}{\Gamma(2-\alpha)\cos(\pi\alpha/2)}$, and

$H(u):=H_{-}(-u)+H_{+}(u)=C_{\alpha}u^{-\alpha},\qquad\text{for any $u>0$}.$

Let us assume that the function $w\mapsto f(w)$ satisfies the following
conditions:

######  Assumption F.1.

(i) $f\in\mathcal{C}^{1}(\mathbb{R})\cap\mathcal{C}^{3}([-K,K])$ for some
$K>0$;

(ii) $f$ has exactly $n$ local minima $m_{i}$, $1\leq i\leq n$ and $n-1$ local
maxima $s_{i}$, $1\leq i\leq n-1$, enumerated in increasing order with
$s_{0}=-\infty$ and $s_{n}=+\infty$:

$\displaystyle-\infty<m_{1}<s_{1}<m_{2}<\cdots<s_{n-1}<m_{n}<+\infty.$ (F.3)

All extrema of $f$ are non-degenerate, i.e. $\partial_{w}^{2}f(m_{i})>0$,
$1\leq i\leq n$, and $\partial_{w}^{2}f(s_{i})<0$, $1\leq i\leq n-1$.

(iii) $|\partial_{w}f(w)|>c_{1}|w|^{1+c_{2}}$ as $w\rightarrow\pm\infty$ for
some $c_{1},c_{2}>0$.

First, we consider the first exit time from a single well. For $\varepsilon>0$
and $\gamma>0$, define

$\Omega_{\varepsilon}^{i}:=\left[s_{i-1}+2\varepsilon^{\gamma},s_{i}-2\varepsilon^{\gamma}\right],$
(F.4)

with the convention that
$\Omega_{\varepsilon}^{1}:=(-\infty,s_{1}-2\varepsilon^{\gamma}]$ and
$\Omega_{\varepsilon}^{n}:=[s_{n-1}+2\varepsilon^{\gamma},+\infty)$. The first
exit time from the $i$-th well is defined as

$\displaystyle\sigma^{i}(\varepsilon;\theta):=\inf\\{t\geq
0:w_{t}\notin[s_{i-1}+\varepsilon^{\gamma},s_{i}-\varepsilon^{\gamma}]\\},$
(F.5)

for $i=1,\ldots,n$. Let us also define

$\displaystyle\lambda^{i}(\varepsilon;\theta):=\frac{1-\theta}{2}C_{\alpha}\left|\frac{s_{i-1}-m_{i}}{\varepsilon}\right|^{-\alpha}+\frac{1+\theta}{2}C_{\alpha}\left|\frac{s_{i}-m_{i}}{\varepsilon}\right|^{-\alpha},$
(F.6)

for $i=1,\ldots,n$. We have the following first exit time result from Imkeller
& Pavlyukevich (2008).

######  Proposition F.1 (Proposition 3.1. in (Imkeller & Pavlyukevich, 2008)).

There exists $\gamma_{0}>0$ such that for any $0<\gamma\leq\gamma_{0}$,
$i=1,2,\ldots,n$,

$\displaystyle\lambda^{i}(\varepsilon;\theta)\sigma^{i}(\varepsilon;\theta)\rightarrow\exp(1),\qquad\text{in
distribution as $\varepsilon\rightarrow 0$},$ (F.7)

where $\exp(1)$ denotes the exponential distribution with mean $1$, and

$\displaystyle\lim_{\varepsilon\rightarrow
0}\mathbb{E}_{w}\left[\lambda^{i}(\varepsilon;\theta)\sigma^{i}(\varepsilon;\theta)\right]=1,$
(F.8)

where the limit holds uniformly over $w\in\Omega_{\varepsilon}^{i}$.

The above result implies that as $\varepsilon\rightarrow 0$,

$\displaystyle\mathbb{E}_{w}\left[\sigma^{i}(\varepsilon;\theta)\right]\sim\Bigg{(}\frac{1-\theta}{2}C_{\alpha}\left|s_{i-1}-m_{i}\right|^{-\alpha}+\frac{1+\theta}{2}C_{\alpha}\left|s_{i}-m_{i}\right|^{-\alpha}\Bigg{)}^{-1}\varepsilon^{-\alpha}.$
(F.9)

If $|s_{i}-m_{i}|>|s_{i-1}-m_{i}|$, i.e. the $i$-th well is asymmetric and the
local minimum $m_{i}$ is closer to the saddle point on the left $s_{i-1}$ than
the saddle point on the right $s_{i}$, then
$\lambda^{i}(\varepsilon;\theta)<\lambda^{i}(\varepsilon;0)$ and
$\mathbb{E}_{w}[\sigma^{i}(\varepsilon;\theta)]>\mathbb{E}_{w}[\sigma^{i}(\varepsilon;0)]$
for positive $\theta$ and
$\lambda^{i}(\varepsilon;\theta)>\lambda^{i}(\varepsilon;0)$ and
$\mathbb{E}_{w}[\sigma^{i}(\varepsilon;\theta)]<\mathbb{E}_{w}[\sigma^{i}(\varepsilon;0)]$
for negative $\theta$. Similarly, if $|s_{i}-m_{i}|<|s_{i-1}-m_{i}|$, i.e. the
$i$-th well is asymmetric and the local minimum $m_{i}$ is closer to the
saddle point on the right $s_{i}$ than the saddle point on the left $s_{i-1}$,
then $\lambda^{i}(\varepsilon;\theta)>\lambda^{i}(\varepsilon;0)$ and
$\mathbb{E}_{w}[\sigma^{i}(\varepsilon;\theta)]<\mathbb{E}_{w}[\sigma^{i}(\varepsilon;0)]$
for positive $\theta$ and
$\lambda^{i}(\varepsilon;\theta)<\lambda^{i}(\varepsilon;0)$ and
$\mathbb{E}_{w}[\sigma^{i}(\varepsilon;\theta)]>\mathbb{E}_{w}[\sigma^{i}(\varepsilon;0)]$
for negative $\theta$. The intuition is that when the well is asymmetric, the
dynamics can exit the well faster when there is a skewness $\theta$ towards
the the saddle point closer to the minimum of the well.

Next, we consider transitions between the wells. For any
$0<\Delta<\Delta_{0}:=\min_{1\leq i\leq n}\\{|m_{i}-s_{i-1}|,|m_{i}-s_{i}|\\}$
and $w\in\mathbb{R}$ denote $B_{\Delta}(w):=\\{v:|w-v|\leq\Delta\\}$. Define

$\tau^{i}(\varepsilon;\theta):=\inf\left\\{t\geq 0:w_{t}\in\cup_{k\neq
i}B_{\Delta}(m_{k})\right\\}.$ (F.10)

Then, we have the following result about transitions between the wells from
Imkeller & Pavlyukevich (2008).

######  Proposition F.2 (Proposition 4.3. in (Imkeller & Pavlyukevich, 2008)).

For any $0<\Delta<\Delta_{0}$ and $j\neq i$

$\displaystyle\lim_{\varepsilon\rightarrow
0}\mathbb{P}_{w}\left(w_{\tau^{i}(\varepsilon;\theta)}\in
B_{\Delta}(m_{j})\right)=\frac{q_{ij}}{q_{i}},$ (F.11)

uniformly for $w\in B_{\Delta}(m_{i})$, $i=1,\ldots,n$, where

$\displaystyle
q_{ij}=\left(\frac{1-\theta}{2}1_{j<i}+\frac{1+\theta}{2}1_{j>i}\right)\cdot\left\|s_{j-1}-m_{i}|^{-\alpha}-|s_{j}-m_{i}|^{-\alpha}\right|,\quad
i\neq j,$ (F.12) $\displaystyle-q_{ii}=q_{i}=\sum_{j\neq
i}q_{ij}=\frac{1-\theta}{2}|s_{i-1}-m_{i}|^{-\alpha}+\frac{1+\theta}{2}|s_{i}-m_{i}|^{-\alpha}.$
(F.13)

From (F.11)-(F.13), we can compute that

$\displaystyle\frac{q_{ij}}{q_{i}}=\begin{cases}\frac{\left\|s_{j-1}-m_{i}|^{-\alpha}-|s_{j}-m_{i}|^{-\alpha}\right|}{\frac{1-\theta}{1+\theta}|s_{i-1}-m_{i}|^{-\alpha}+|s_{i}-m_{i}|^{-\alpha}}&\text{if
$j>i$},\\\
\frac{\left\|s_{j-1}-m_{i}|^{-\alpha}-|s_{j}-m_{i}|^{-\alpha}\right|}{|s_{i-1}-m_{i}|^{-\alpha}+\frac{1+\theta}{1-\theta}|s_{i}-m_{i}|^{-\alpha}}&\text{if
$j<i$}.\end{cases}$ (F.14)

Therefore, $q_{ij}/q_{i}$ is increasing in $\theta$ for $j>i$ and decreasing
in $\theta$ for $j<i$. This is consistent with the intuition that when
$\theta>0$, it is more likely for the dynamics to transit to a well on the
right side, and when $\theta<0$, it is more likely for the dynamics to transit
to a well on the left side.

Next, we consider the following metastability result due to Theorem 1.1. in
Imkeller & Pavlyukevich (2008). It describes the metastability phenomenon,
which basically says that there exists a time scale under which the system
behaves like a continuous time Markov process with a finite state space
consisting of values in the set of stable attractors.

######  Theorem F.1 (Theorem 1.1. in (Imkeller & Pavlyukevich, 2008)).

If $w_{0}=w\in(s_{i-1},s_{i})$ for some $i=1,2,\ldots,n$, then for any $t>0$,
in the sense of finite-dimensional distributions,

$\displaystyle w_{t/H(1/\varepsilon)}\rightarrow Y_{t}(m_{i}),\qquad\text{as
$\varepsilon\rightarrow 0$},$ (F.15)

where $w_{t}$ is defined in (F.1) and
$H(1/\varepsilon)=C_{\alpha}\varepsilon^{\alpha}$, where $Y_{t}(m_{i})$ that
starts at $m_{i}$ is a continuous-time Markov process on a finite states space
$\\{m_{1},\ldots,m_{n}\\}$ with the infinitesimal generator
$Q=(q_{ij})_{i,j=1}^{n}$, where $q_{ij}$ is defined in (F.12).

The Markov process $Y_{t}(m_{i})$ admits a unique invariant distribution $\pi$
satisfying $Q^{T}\pi=0$. In the case of double well, i.e. $n=2$ and
$m_{1}<s_{1}=0<m_{2}$ separated by a local maximum at $s_{1}=0$, where without
loss of generality we assume that $m_{2}>|m_{1}|$, i.e. the second local
minimum lies in a wider valley. A simple calculation yields that

$q_{12}=\frac{1+\theta}{2}\frac{1}{|m_{1}|^{\alpha}}=-q_{11},\quad\text{and}\quad
q_{21}=\frac{1-\theta}{2}\frac{1}{m_{2}^{\alpha}}=-q_{22},$ (F.16)

so that it follows from $Q^{T}\pi=0$ and $\pi_{1}+\pi_{2}=1$ that

$\displaystyle\pi_{1}=\frac{(1+\theta)^{-1}|m_{1}|^{\alpha}}{(1+\theta)^{-1}|m_{1}|^{\alpha}+(1-\theta)^{-1}m_{2}^{\alpha}},$
(F.17)
$\displaystyle\pi_{2}=\frac{(1-\theta)^{-1}m_{2}^{\alpha}}{(1+\theta)^{-1}|m_{1}|^{\alpha}+(1-\theta)^{-1}m_{2}^{\alpha}}.$
(F.18)

In particular, the ratio
$\frac{\pi_{2}}{\pi_{1}}=\frac{1+\theta}{1-\theta}\left(\frac{m_{2}}{|m_{1}|}\right)^{\alpha}$
is increasing in $\frac{m_{2}}{|m_{1}|}$ and $\theta$. That reveals that in
the equilibrium the process will spend more time in the second valley if the
second valley is wide and there is a drift towards to the right. In the
symmetric case, i.e. $\theta=0$, $\pi_{2}>\pi_{1}$ since $m_{2}>|m_{1}|$ so
that in the equilibrium the process spends more time in the wider valley. In
the asymmetric case, i.e. $\theta\neq 0$, if there is a strong skewness
towards the left, i.e. $\theta<0$ and $|\theta|$ is large, then in the
equilibrium the process may spend more time in the narrower valley. Indeed
$\pi_{2}>\pi_{1}$ if and only if
$\theta>\frac{|m_{1}|^{\alpha}-m_{2}^{\alpha}}{|m_{1}|^{\alpha}+m_{2}^{\alpha}}$.

##  Appendix G Postponed Proofs

###  G.1 Proof of Lemma 3.1

Before we proceed to the proof of Lemma 3.1 we present some intermediary
results that are required for the proof.

######  Definition G.1.

(Asymptotic order of magnitude) A positive sequence $a_{m}$ is of the same
order of magnitude as another positive sequence $b_{m}$ ($a_{m}\asymp b_{m}$,
i.e. ‘asymptotically equivalent’) if there exist some $c,C>0$ such that:
$c\leq\frac{a_{m}}{b_{m}}\leq C$ for any $m\in\mathbb{N}$.

######  Lemma G.1 (Lemma A.1 in (Vladimirova et al., 2019)).

Let $X$ be a normal random variable such that
$X\sim\mathcal{N}(0,\sigma^{2})$. Then the following asymptotic equivalence
holds

$\|X\|_{m}\asymp\sqrt{m}.$

We know that the centering of variables does not change their tail properties
(Vershynin, 2018; Kuchibhotla & Chakrabortty, 2018). As such Lemma G.1 also
applies to $X\sim\mathcal{N}(\mu,\sigma^{2})$, as
$\|X\|\asymp\|X-\mu\|\asymp\sqrt{m}$.

######  Lemma G.2 (Lemma 3.1 of Vladimirova et al. (2019)).

Let $\kappa:\mathbb{R}\to\mathbb{R}$ be a non-linearity that obeys the
extended envelope property. And let $X$ be a variable for which
$\|X_{+}\|_{m}\asymp\|X_{-}\|_{m}$ where $X_{-}$ and $X_{+}$ denote the left
and right tail of the variable respectively 444We weaken Vladimirova et al.
(2019)’s requirement for $X$ to be symmetric as the proof they give still
holds here.. Then we have:

$\displaystyle\|\kappa(X)\|_{m}\asymp\|X\|_{m}\,,\qquad\text{for any $m\geq
1$}\,.$ (G.1)

######  Lemma G.3.

Let $X_{1},\dots,X_{N}$ be variables that each obeys $\|X_{i}\|_{m}\lesssim
m^{r},p\in R,i=1,\dots,N$ and $(W_{i},\dots,W_{N})\in\mathbb{R}^{N}$.

$\left\|\sum_{i=1}^{N}W_{i}X_{i}\right\|_{m}\lesssim m^{r}\,.$

###### Proof of Lemma G.3.

By Minkowski’s inequality we have that

$\displaystyle\left\|\sum_{i=1}^{N}W_{i}X_{i}\right\|_{m}$
$\displaystyle\leq\sum_{i=1}^{N}\|W_{i}X_{i}\|_{m}\leq\sum_{l=1}^{N}|W_{i}A_{i}|m^{r},\
(A_{1},\dots,A_{N})\in\mathbb{R}^{N}$
$\displaystyle\Leftrightarrow\left\|\sum_{i=1}^{N}W_{i}X_{i}\right\|_{m}$
$\displaystyle\lesssim m^{r}\,.$

The $A_{i}$ here are constants that upper bound the asymptotics of each norm
$\|X_{i}\|_{m}$ in the sum. ∎

###### Proof of Lemma 3.1.

Additive Noise. Consider first the noised data,
$\widetilde{\boldsymbol{\mathbf{h}}}_{0}(\boldsymbol{\mathbf{x}})=\boldsymbol{\mathbf{x}}+\boldsymbol{\mathbf{\epsilon}}_{0},\boldsymbol{\mathbf{\epsilon}}_{0}\sim\mathcal{N}(0,\sigma_{0}^{2})$.
As Lemma G.1 shows, Gaussian random variables have an $m^{\mathrm{th}}$ norm
that is asymptotically equivalent to $\sqrt{m}$,

$\left\|\widetilde{h}_{0,l}(\boldsymbol{\mathbf{x}})\right\|_{m}\asymp\sqrt{m},\qquad\text{for
any $l=1,\dots,n_{0}$}\,,$

where $n_{0}$ is the dimensionality of data.

Let us now assume that
$\left\|\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})\right\|_{m}\lesssim\sqrt{m}$,
for any $l=1,\dots,n_{i}$, for some layer $i$. The pre-non-linearity at this
layer is given by
$\tilde{\boldsymbol{\mathbf{g}}}_{i}=\boldsymbol{\mathbf{W}}_{i+1}\widetilde{\boldsymbol{\mathbf{h}}}_{i}$.
The $j^{\mathrm{th}}$ element of $\tilde{\boldsymbol{\mathbf{g}}}_{i}$ is
defined as a sum,

$g_{i,j}(\boldsymbol{\mathbf{x}})=\sum_{l=1}^{n_{i}}W_{i+1,l,j}\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})\,,$

where $W_{i+1,l,j}$ is the weight that maps from the $l^{\mathrm{th}}$ neuron
in layer $i$ to the $j^{\mathrm{th}}$ in layer $i+1$. By Lemma G.3,

$\displaystyle\left\|g_{i,j}(\boldsymbol{\mathbf{x}})\right\|_{m}$
$\displaystyle\lesssim\sqrt{m},\qquad m=1,\dots,n_{i}\,.$

As such if we assume the non-linearities $\phi$ at each layer obey the
extended envelope property, then we have by Lemma G.2:

$\displaystyle\left\|\phi(g_{i,j}(\boldsymbol{\mathbf{x}}))\right\|_{m}=\left\|\widehat{h}_{i+1,j}(\boldsymbol{\mathbf{x}})\right\|_{m}$
$\displaystyle\asymp\left\|\kappa(g_{i,j}(\boldsymbol{\mathbf{x}}))\right\|_{m}$
$\displaystyle\Leftrightarrow\left\|\widehat{h}_{i+1,j}(\boldsymbol{\mathbf{x}})\right\|_{m}$
$\displaystyle\lesssim\sqrt{m},\qquad j=1,\dots,n_{i+1}\,.$

Note that
$\widetilde{\boldsymbol{\mathbf{h}}}_{i+1}=\widehat{\boldsymbol{\mathbf{h}}}_{i+1}+\boldsymbol{\mathbf{\epsilon}}_{i+1},\boldsymbol{\mathbf{\epsilon}}_{i+1}\sim\mathcal{N}(0,\sigma_{i+1}^{2})$.
By Lemma G.3, once again
$\|\widetilde{h}_{i+1,j}(\boldsymbol{\mathbf{x}})\|_{m}\lesssim\sqrt{m}$, is
Gaussian in its tails. By recursion, with
$\widetilde{\boldsymbol{\mathbf{h}}}_{0}$ as the base case, we have that

$\displaystyle\left\|\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})\right\|_{m}\lesssim\sqrt{m},\qquad\text{for
any $m\geq 1$};\ i=1,\dots,L-1;\ l=1,\dots,n_{i}\,.$

Multiplicative Noise. Consider first the noised data,
$\widetilde{\boldsymbol{\mathbf{h}}}_{0}(\boldsymbol{\mathbf{x}})=\boldsymbol{\mathbf{x}}\circ\boldsymbol{\mathbf{\epsilon}}_{0},\boldsymbol{\mathbf{\epsilon}}_{0}\sim\mathcal{N}(1,\sigma_{0}^{2})$.
As Lemma G.1 shows, Gaussian random variables have an $m^{\mathrm{th}}$ norm
that is asymptotically equivalent to $\sqrt{m}$,

$\left\|\widetilde{h}_{0,l}(\boldsymbol{\mathbf{x}})\right\|_{m}\asymp\sqrt{m},\qquad\text{for
any $l=1,\dots,n_{0}$}\,,$

where $n_{0}$ is the dimensionality of data.

Let us now assume that
$\left\|\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})\right\|_{m}\lesssim
m^{r},\ \forall l=1,\dots,n_{i}$, for some layer $i$. The pre-non-linearity at
this layer is given by
$\tilde{\boldsymbol{\mathbf{g}}}_{i}=\boldsymbol{\mathbf{W}}_{i+1}\widetilde{\boldsymbol{\mathbf{h}}}_{i}$.
The $j^{\mathrm{th}}$ element of $\tilde{\boldsymbol{\mathbf{g}}}_{i}$ is
defined as a sum,

$g_{i,j}(\boldsymbol{\mathbf{x}})=\sum_{l=1}^{n_{i}}W_{i+1,l,j}\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})\,,$

where $W_{i+1,l,j}$ is the weight that maps from the $l^{\mathrm{th}}$ neuron
in layer $i$ to the $j^{\mathrm{th}}$ in layer $i+1$. By Lemma G.3,

$\displaystyle\left\|g_{i,j}(\boldsymbol{\mathbf{x}})\right\|_{m}$
$\displaystyle\lesssim m^{r},\qquad m=1,\dots,n_{i}\,.$

As such if we assume the non-linearities $\phi$ at each layer obey the
extended envelope property, then we have by Lemma G.2:

$\displaystyle\left\|\phi(g_{i,j}(\boldsymbol{\mathbf{x}}))\right\|_{m}=\left\|\widehat{h}_{i+1,j}(\boldsymbol{\mathbf{x}})\right\|_{m}$
$\displaystyle\asymp\left\|\kappa(g_{i,j}(\boldsymbol{\mathbf{x}}))\right\|_{m}$
$\displaystyle\Leftrightarrow\left\|\widehat{h}_{i+1,j}(\boldsymbol{\mathbf{x}})\right\|_{m}$
$\displaystyle\lesssim m^{r},\qquad j=1,\dots,n_{i+1}\,.$

Note that
$\widetilde{\boldsymbol{\mathbf{h}}}_{i+1}=\widehat{\boldsymbol{\mathbf{h}}}_{i+1}\circ\boldsymbol{\mathbf{\epsilon}}_{i+1},\boldsymbol{\mathbf{\epsilon}}_{i+1}\sim\mathcal{N}(1,\sigma_{i+1}^{2})$.
By Hölder’s inequality we have that,

$\left\|\widetilde{\boldsymbol{\mathbf{h}}}_{i+1,j}(\boldsymbol{\mathbf{x}})\right\|_{m}\lesssim
m^{r+\frac{1}{2}},\qquad j=1,\dots,n_{i+1}\,.$

By recursion, with $\widetilde{\boldsymbol{\mathbf{h}}}_{0}$ as the base case,
we have that

$\displaystyle\left\|\widetilde{h}_{i,l}(\boldsymbol{\mathbf{x}})\right\|_{m}\lesssim
m^{\frac{i+1}{2}},\qquad\text{for any $m\geq 1$};\ i=1,\dots,L-1;\
l=1,\dots,n_{i}\,.$

∎

###  G.2 Proof of Theorem 3.1

Before we proceed to the proof of Theorem 3.1 we present some intermediary
results that are required for the proof.

######  Lemma G.4.

Let $X$ be a bounded random variable such that $|X|\leq C$, then $X$ is sub-
Weibull with parameter $\theta=0$,

$\displaystyle\left\|X\right\|_{m}\lesssim m^{0},\qquad\text{for every $m\geq
1$}.$ (G.2)

###### Proof of Lemma G.4.

The moments of $X$ obey

$\displaystyle\operatorname{\mathbb{E}}\left[|X|^{m}\right]\leq C^{k}\,.$

Taking the root of this we find that
$\|X\|_{m}:=\operatorname{\mathbb{E}}[|X|^{m}]^{\frac{1}{m}}$ is not dependent
on $m$ and scales as a constant, $\left\|X\right\|_{m}\asymp m^{0}$. ∎

###### Proof of Theorem 3.1.

We first consider the gradient for a single data-label pair of $W_{i,l,j}$ the
weight that maps from neuron $l$ in layer $i-1$ to neuron $j$ in layer $i$. We
study this gradient using the chain rule, where we decompose $\partial
E_{\mathcal{L}}(\cdot)/\partial W_{i,l,j}$ as

$\frac{\partial
E_{\mathcal{L}}(\cdot)}{\partial\widetilde{h}_{i,j}}\cdot\frac{\partial\widetilde{h}_{i,j}}{\partial
W_{i,l,j}},$

where $\widetilde{h}_{i,j}$ is the (noised) activation of the
$j^{\mathrm{th}}$ neuron in the $i^{\mathrm{th}}$ layer. Thus, the gradient
noise on the weights can be described as the product of two random variables.

Additive Noise. Let us first consider the properties of $\partial
E_{\mathcal{L}}(\cdot)/\partial\widetilde{h}_{i,j}$ for the additive case.

Regression. In the case of regression we use a mean-square-error (MSE) we have
that:

$\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})=2(\boldsymbol{\mathbf{y}}-\boldsymbol{\mathbf{h}}_{L}(\boldsymbol{\mathbf{x}}))\mathcal{E}_{L}(\boldsymbol{\mathbf{x}})+\left(\mathcal{E}_{L}(\boldsymbol{\mathbf{x}})\right)^{2}\,,$

where we imply all terms’ dependence on
$\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}}$ for brevity of
notation. One can already see that the derivative of this object with respect
to each element of $\widetilde{\boldsymbol{\mathbf{h}}}_{L}$ will have tail
properties that are asymptotically equivalent to those of $\mathcal{E}_{L}$,
which we know by Lemma 3.1.

$\left\|\frac{\partial\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial\widetilde{h}_{L,j}}\right\|_{m}\lesssim\sqrt{m},\qquad
m=1,\dots,n_{L}\,.$

If we center this distribution the tail properties of this variable are
unchanged (Vershynin, 2018; Kuchibhotla & Chakrabortty, 2018). In particular
the asymptotic behaviour of $\|\cdot\|_{m}$ is unchanged and we have that:

$\displaystyle\left\|\frac{\partial\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial\widetilde{h}_{L,j}}-\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}\left[\frac{\partial\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial\widetilde{h}_{L,j}}\right]\right\|_{m}$
$\displaystyle\lesssim\sqrt{m},$
$\displaystyle\Leftrightarrow\left\|\frac{\partial
E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial\widetilde{h}_{L,j}}\right\|_{m}$
$\displaystyle\lesssim\sqrt{m},\qquad m=1,\dots,n_{L}\,.$

Classification. In the case of classification we use a cross-entropy (CE)
error. There is no easy closed-form for
$\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$ here, but
we can infer the properties of
$\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{L}}\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$
by studying the properties of the gradient
$\nabla_{\boldsymbol{\mathbf{h}}_{L}}\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$.
For CE we know that:

$\nabla_{\boldsymbol{\mathbf{h}}_{L}}\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})=\mathrm{sigmoid}(\boldsymbol{\mathbf{h}}_{L}(\boldsymbol{\mathbf{x}}))-\boldsymbol{\mathbf{y}}$

in the binary label case. In the multi-label classification case we typically
use a $\mathrm{softmax}$, which is also a bounded function. We can already see
that any noise $\mathcal{E}_{L}$ added to
$\boldsymbol{\mathbf{h}}_{L}(\boldsymbol{\mathbf{x}})$ will induce a change in
the gradient that is inherently bounded by the $\mathrm{sigmoid}$ non-
linearity, meaning that
$\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$ will be
bounded. As such the centered variable
$E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$ will also
be bounded and zero mean,. By Lemma G.4 any bounded and zero mean distribution
will be sub-Weibull with parameter $\theta=0$, and will thus also be sub-
Gaussian

$\left\|\frac{\partial
E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial\widetilde{h}_{L,j}}\right\|_{m}\lesssim\sqrt{m},\qquad
m=1,\dots,n_{L}\,.$

Synthesizing the regression and classification settings we can conclude that
each constitutive element of
$\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{L}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$
will be sub-Gaussian and will have zero mean.

We can now turn to the partial derivatives of the form $\partial
E_{\mathcal{L}}(\cdot)/\partial\widetilde{h}_{i}^{m}$. Assume
$\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})$
is of order $m^{r}$, with,

$\left\|\frac{\partial
E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial\widetilde{h}_{i,j}}\right\|_{m}\lesssim
m^{r},\qquad m=1,\dots,n_{i}\,,$

which entails for gradients at the previous $(i-1)^{\mathrm{th}}$ layer we
have

$\frac{\partial
E_{\mathcal{L}}(\cdot)}{\partial\widetilde{h}_{i-1,l}}=\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})\nabla_{\widetilde{h}_{i-1,l}}\widetilde{\boldsymbol{\mathbf{h}}}_{i}(\boldsymbol{\mathbf{x}})\,,$

where

$\nabla_{\widetilde{h}_{i-1,l}}\widetilde{\boldsymbol{\mathbf{h}}}_{i}(\boldsymbol{\mathbf{x}})=\kappa^{\prime}\left(\boldsymbol{\mathbf{W}}_{i,l}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\circ(\boldsymbol{\mathbf{W}}_{i,l})\,,$

where $\circ$ denotes the element wise product and
$\boldsymbol{\mathbf{W}}_{i,l}$ is the $l^{\mathrm{th}}$ column of the weight
matrix $\boldsymbol{\mathbf{W}}_{i}$. By definition, activation functions that
obey the extended envelope property will have gradients that are bounded in
norm, by some constant $d_{2}$. As such, by Lemma G.4 $\kappa^{\prime}$ will
be sub-Weibull with $r=0$. By Hölder’s inequality we have that

$\left\|\left(\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})\right)_{z}\left(\kappa^{\prime}\left(\boldsymbol{\mathbf{W}}_{i,l}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\right)_{j}\right\|_{m}\leq\left\|\left(\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})\right)_{j}\right\|_{2m}\left\|\left(\kappa^{\prime}\left(\boldsymbol{\mathbf{W}}_{i,l}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\right)_{j}\right\|_{2m}\,,$

where $j$ indexes over the elements of both Jacobians. By definition, we know
that there exists $A>0$ and $B>0$ such that
$\|(\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}}))_{z}\|_{2m}\leq
A(2m)^{p}$ and
$\|(\kappa^{\prime}(\boldsymbol{\mathbf{W}}_{i,l}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})))_{j}\|_{2m}\leq
B$. As such

$\left\|\left(\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})\right)_{z}\left(\kappa^{\prime}\left(\boldsymbol{\mathbf{W}}_{i,l}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\right)_{j}\right\|_{m}\leq\left\|\left(\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})\right)_{j}\right\|_{2m}\left\|\left(\kappa^{\prime}\left(\boldsymbol{\mathbf{W}}_{i,l}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\right)_{j}\right\|_{2m}\leq
AB2^{p}m^{r}\,.$

Thus we know that the product of these two variables will be asymptotically
upper-bounded by

$\left\|\left(\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})\right)_{j}\left(\kappa^{\prime}\left(\boldsymbol{\mathbf{W}}_{i,l}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\right)_{j}\right\|_{m}\lesssim
m^{r}\,.$

We now need to take into account the weighted sum across rows and columns
(i.e. over indices $j$) that occurs. By Lemma G.3 we know that

$\displaystyle\left\|\sum_{j=1}^{n_{i}}\left(\nabla_{\widetilde{\boldsymbol{\mathbf{h}}}_{i}}E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})\right)_{j}\left(\kappa^{\prime}\left(\boldsymbol{\mathbf{W}}_{i,l}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\circ(\boldsymbol{\mathbf{W}}_{i,l})\right)_{j}\right\|_{m}\lesssim
m^{r},$ $\displaystyle\Leftrightarrow\left\|\frac{\partial
E_{\mathcal{L}}(\cdot)}{\partial\widetilde{h}_{i-1,l}}\right\|_{m}\lesssim
m^{r},\qquad m=1,\dots,n_{i}\,.$

By recursion, with $\widetilde{\boldsymbol{\mathbf{h}}}_{L}$ as the base case,
gradients at layer $i$ bounded in norm by $m^{r}$ induce gradients at layer
$i-1$, also bounded in norm by $m^{r}$. By recursion with the
$L^{\mathrm{th}}$ layer as the base case we can claim that,

$\displaystyle\left\|\frac{\partial
E_{\mathcal{L}}(\cdot)}{\partial\widetilde{h}_{i,j}}\right\|_{m}\lesssim\sqrt{m}\,.$

We have now defined the first constitutive term of $\partial
E_{\mathcal{L}}(\cdot)/\partial W_{i,l,j}$. Defining
$\partial\widetilde{h}_{i,j}/\partial W_{i,l,j}$ is much simpler:

$\frac{\partial\widetilde{h}_{i,j}}{\partial
W_{i,l,j}}=\kappa^{\prime}\left(W_{i,l,j}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\left(\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\,.$

Here $\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})$, which we know is sub-
Gaussian by Lemma 3.1, is once again multiplied to a bounded variable,
$\kappa^{\prime}$. Thus reapplying Hölder’s inequality we obtain that

$\left\|\frac{\partial\widetilde{h}_{i,j}}{\partial
W_{i,l,j}}\right\|_{m}\lesssim\sqrt{m}\,.$

We can now bring together the characterisations of the gradients that
constitute $\partial E_{\mathcal{L}}(\cdot)/\partial W_{i,l,j}$. We can re-use
Hölder’s inequality to show that the product of these variables will be sub-
exponential

$\left\|\frac{\partial
E_{\mathcal{L}}((\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}});\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})}{\partial
W_{i,l,j}}\right\|_{m}\lesssim m,\qquad\text{for every $m\geq 1$}\,.$

Multiplicative Noise. In the case of multiplicative noise we know that by
Lemma 3.1, the accumulated noise at layer $L$ will be of the same order as
that at layer $L-1$, because we are not multiplying noise to the final layer,
thus

$\left\|\frac{\partial
E_{\mathcal{L}}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial\widetilde{h}_{L,j}}\right\|_{m}\lesssim
m^{\frac{{L}}{2}},\qquad m=1,\dots,n_{L}\,.$

Repeating the analysis done for the additive case we can claim that,

$\displaystyle\left\|\frac{\partial
E_{\mathcal{L}}(\cdot)}{\partial\widetilde{h}_{i}^{m}}\right\|_{m}\lesssim
m^{\frac{{L}}{2}}.$

We have now defined the first constitutive term of $\partial
E_{\mathcal{L}}(\cdot)/\partial W_{i,l,j}$. We now define
$\partial\widetilde{h}_{i,j}/\partial W_{i,l,j}$.

$\frac{\partial\widetilde{h}_{i}^{m}}{\partial
W_{i,l,j}}=\kappa^{\prime}\left(W_{i,l,j}\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\left(\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})\right)\,.$

Here $\widetilde{h}_{i-1,l}(\boldsymbol{\mathbf{x}})$, which we know is sub-
Weibull with a parameter $p=\frac{i}{2}$ by Lemma 3.1, is once again
multiplied to a bounded variable, $\kappa^{\prime}$. Thus reapplying Hölder’s
inequality we obtain that

$\left\|\frac{\partial\widetilde{h}_{i}^{m}}{\partial
W_{i,l,j}}\right\|_{m}\lesssim m^{\frac{i}{2}}.$

We can re-use Hölder’s inequality to show that the product of these variables
will be sub-Weibull with $p=\frac{L+i}{2}$

$\left\|\frac{\partial
E_{\mathcal{L}}((\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}});\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})}{\partial
W_{i,l,j}}\right\|_{m}\lesssim m^{\frac{L+i}{2}},\qquad\text{for every $m\geq
1$}\,.$

Mean of Noise. Finally, by definition these gradients are zero mean,

$\frac{\partial
E_{\mathcal{L}}((\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}});\boldsymbol{\mathbf{w}},\boldsymbol{\mathbf{\epsilon}})}{\partial
W_{i,l,j}}=\frac{\partial\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial
W_{i,l,j}}-\operatorname{\mathbb{E}}_{\boldsymbol{\mathbf{\epsilon}}}\left[\frac{\partial\Delta\mathcal{L}(\boldsymbol{\mathbf{x}},\boldsymbol{\mathbf{y}})}{\partial
W_{i,l,j}}\right].$

∎

###  G.3 Proof of Theorem 4.1

Before we proceed to the proof of Theorem 4.1, we present some technical
results that will be used in the proof of Theorem 4.1 later.

For the $d$-dimensional asymmetric fractional Langevin dynamics
$\boldsymbol{\mathbf{w}}_{t}$, its infinitesimal generator is given in the
following proposition.

######  Proposition G.1.

The asymmetric fractional Langevin dynamics $\boldsymbol{\mathbf{w}}_{t}$ has
the infinitesimal generator:

$\displaystyle\mathcal{L}f(\boldsymbol{\mathbf{w}})=\sum_{i=1}^{d}\left((b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}-\varepsilon^{\alpha}\frac{\theta_{i}}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}\right)\frac{\partial
f(\boldsymbol{\mathbf{w}})}{\partial
w_{i}}+\varepsilon^{\alpha}\sum_{i=1}^{d}\mathcal{H}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}}),$
(G.3)

where

$\displaystyle\mathcal{H}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}})$
$\displaystyle:=\left(\frac{1}{2}+\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{f(\boldsymbol{\mathbf{w}}+\xi\mathbf{e}_{i})-f(\boldsymbol{\mathbf{w}})-\partial_{w_{i}}f(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi$
$\displaystyle\qquad+\left(\frac{1}{2}-\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{f(\boldsymbol{\mathbf{w}}-\xi\mathbf{e}_{i})-f(\boldsymbol{\mathbf{w}})+\partial_{w_{i}}f(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi,$
(G.4)

where $\mathbf{e}_{i}$ is the $i$-th basis vector in $\mathbb{R}^{d}$, i.e. a
$d$-dimensional unit vector with $i$-th coordinate being $1$ and all the other
coordinates being $0$.

###### Proof of Proposition G.1.

Since the asymmetric fractional Langevin dynamics is driven by the
$d$-dimensional $\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta}$, it suffices to
show that the infinitesimal generator of $d$-dimensional
$\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta}$ is given by

$\displaystyle\sum_{i=1}^{d}\mathcal{G}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}})=\sum_{i=1}^{d}\mathcal{H}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}})-\sum_{i=1}^{d}\frac{\theta_{i}}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}\frac{\partial}{\partial
w_{i}}f(\boldsymbol{\mathbf{w}}).$ (G.5)

We start the proof by considering the dimension $d=1$ first. The one-
dimensional $\alpha$-stable Lévy motion with tail-index $1<\alpha<2$ and
skewness $\theta\in(-1,1)$ has the infinitesimal generator given by:

$\displaystyle\mathcal{G}^{\alpha,\theta}f(w)$
$\displaystyle:=\frac{1+\theta}{2}\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{z>0}[f(w+z)-f(w)-1_{|z|\leq
1}f^{\prime}(w)z]\frac{dz}{|z|^{1+\alpha}}$
$\displaystyle\quad+\frac{1-\theta}{2}\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{z<0}[f(w+z)-f(w)-1_{|z|\leq
1}f^{\prime}(w)z]\frac{dz}{|z|^{1+\alpha}}+af^{\prime}(w),$ (G.6)

where $a\in\mathbb{R}$ is chosen so that $\mathcal{G}^{\alpha,\theta}w=0$ to
be consistent with $\mu=0$ in
$\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta}-\boldsymbol{\mathbf{L}}_{s}^{\alpha,\theta}\sim\mathcal{S}_{\alpha}((t-s)^{1/\alpha},\theta,\mu)$
for any $t>s$. Thus, we can compute that

$\displaystyle\mathcal{G}^{\alpha,\theta}w=\frac{1+\theta}{2}\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{z>1}\frac{dz}{z^{\alpha}}-\frac{1-\theta}{2}\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{z>1}\frac{dz}{z^{\alpha}}+a=0,$
(G.7)

which yields that

$\displaystyle
a=-\theta\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}.$
(G.8)

Therefore, with $1<\alpha<2$,

$\displaystyle\mathcal{G}^{\alpha,\theta}f(w)=\mathcal{H}^{\alpha,\theta}f(w)-\theta\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}f^{\prime}(w),$
(G.9)

where

$\displaystyle\mathcal{H}^{\alpha,\theta}f(w)$
$\displaystyle:=\left(\frac{1}{2}+\frac{\theta}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{f(w+\xi)-f(w)-f^{\prime}(w)\xi}{\xi^{\alpha+1}}d\xi$
$\displaystyle\qquad\qquad+\left(\frac{1}{2}-\frac{\theta}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{f(w-\xi)-f(w)+f^{\prime}(w)\xi}{\xi^{\alpha+1}}d\xi.$
(G.10)

Similarly, for the multi-dimensional case, we can show that the infinitesimal
generator for $\boldsymbol{\mathbf{L}}_{t}^{\alpha,\theta}$ is given by:

$\displaystyle\sum_{i=1}^{d}\mathcal{G}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}})=\sum_{i=1}^{d}\mathcal{H}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}})-\sum_{i=1}^{d}\frac{\theta_{i}}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}\frac{\partial}{\partial
w_{i}}f(\boldsymbol{\mathbf{w}}),$ (G.11)

where

$\displaystyle\mathcal{H}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}})$
$\displaystyle:=\left(\frac{1}{2}+\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{f(\boldsymbol{\mathbf{w}}+\xi\mathbf{e}_{i})-f(\boldsymbol{\mathbf{w}})-\partial_{w_{i}}f(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi$
$\displaystyle\qquad\qquad+\left(\frac{1}{2}-\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{f(\boldsymbol{\mathbf{w}}-\xi\mathbf{e}_{i})-f(\boldsymbol{\mathbf{w}})+\partial_{w_{i}}f(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi,$
(G.12)

where $\mathbf{e}_{i}$ is the $i$-th basis vector in $\mathbb{R}^{d}$, i.e. a
$d$-dimensional unit vector with $i$-th coordinate being $1$ and all the other
coordinates being $0$. ∎

We recall that for any $\theta_{i}\in(-1,1)$, $1\leq i\leq d$, and
$1<\alpha<2$, we have

$\displaystyle(b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}=\frac{\varepsilon^{\alpha}}{\varphi(\boldsymbol{\mathbf{w}})}\mathcal{D}^{\alpha-2,-\theta_{i}}_{w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})),\qquad\varphi(\boldsymbol{\mathbf{w}})=e^{-\varepsilon^{-\alpha}f(\boldsymbol{\mathbf{w}})},$
(G.13)

where

$\displaystyle\mathcal{D}^{\alpha-2,-\theta_{i}}_{w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})):=\frac{-1}{2\cos(\alpha\pi/2)}\left[(1-\theta_{i})\mathcal{I}^{2-\alpha}_{+,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}))+(1+\theta_{i})\mathcal{I}^{2-\alpha}_{-,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}))\right],$
(G.14)

and

$\displaystyle\mathcal{I}^{2-\alpha}_{\pm,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})):=\frac{1}{\Gamma(2-\alpha)}\int_{0}^{\infty}\frac{\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}\pm\xi\mathbf{e}_{i})}{\xi^{\alpha-1}}d\xi.$
(G.15)

In the next result, we provide an alternative formula for
$b(\boldsymbol{\mathbf{w}},\alpha,\theta)=((b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i},1\leq
i\leq d)$ that is defined in (4.7).

######  Proposition G.2.

For any $\theta_{i}\in(-1,1)$, $1\leq i\leq d$, and $1<\alpha<2$, we have

$\displaystyle(b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}$
$\displaystyle:=\frac{\varepsilon^{\alpha}}{\varphi(\boldsymbol{\mathbf{w}})}\left(\frac{1}{2}-\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{\int_{w_{i}}^{w_{i}+\xi}\varphi(\boldsymbol{\mathbf{w}}+(y-w_{i})\mathbf{e}_{i})dy-\varphi(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi$
$\displaystyle\qquad+\frac{\varepsilon^{\alpha}}{\varphi(\boldsymbol{\mathbf{w}})}\left(\frac{1}{2}+\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{\int_{w_{i}}^{w_{i}-\xi}\varphi(\boldsymbol{\mathbf{w}}+(y-w_{i})\mathbf{e}_{i})dy+\varphi(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi$
$\displaystyle\qquad\qquad+\varepsilon^{\alpha}\theta_{i}\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}\,,$
(G.16)

with
$\varphi(\boldsymbol{\mathbf{w}}):=\exp(-\varepsilon^{-\alpha}f(\boldsymbol{\mathbf{w}}))$.

###### Proof of Proposition G.2.

Let us first consider the case when $\theta_{i}=0$, $1\leq i\leq d$. We have

$\displaystyle\mathcal{D}^{\alpha-2}_{w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}))$
$\displaystyle:=-\mathcal{I}^{2-\alpha}_{w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}))$
$\displaystyle=\frac{-1}{2\cos(\alpha\pi/2)}\left[\mathcal{I}^{2-\alpha}_{+,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}))+\mathcal{I}^{2-\alpha}_{-,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}))\right],$
(G.17)

where

$\displaystyle\mathcal{I}^{2-\alpha}_{+,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})):=\frac{1}{\Gamma(2-\alpha)}\int_{0}^{\infty}\frac{\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}+\xi\mathbf{e}_{i})}{\xi^{\alpha-1}}d\xi,$
(G.18)
$\displaystyle\mathcal{I}^{2-\alpha}_{-,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})):=\frac{1}{\Gamma(2-\alpha)}\int_{0}^{\infty}\frac{\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}-\xi\mathbf{e}_{i})}{\xi^{\alpha-1}}d\xi.$
(G.19)

Similarly, when $\theta_{i}\in(-1,1)$, $1\leq i\leq d$, and $1<\alpha<2$, we
have

$\displaystyle(b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}=\frac{\varepsilon^{\alpha}}{\varphi(\boldsymbol{\mathbf{w}})}\mathcal{D}^{\alpha-2,-\theta_{i}}_{w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})),\qquad\varphi(\boldsymbol{\mathbf{w}})=e^{-\varepsilon^{-\alpha}f(\boldsymbol{\mathbf{w}})},$
(G.20)

where

$\displaystyle\mathcal{D}^{\alpha-2,-\theta_{i}}_{w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})):=\frac{-1}{2\cos(\alpha\pi/2)}\left[(1-\theta_{i})\mathcal{I}^{2-\alpha}_{+,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}))+(1+\theta_{i})\mathcal{I}^{2-\alpha}_{-,w_{i}}(\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}}))\right].$
(G.21)

∎

Now, we are ready to prove Theorem 4.1.

###### Proof of Theorem 4.1.

We recall from Proposition G.1 that the asymmetric fractional Langevin
dynamics $\boldsymbol{\mathbf{w}}_{t}$ has the infinitesimal generator:

$\displaystyle\mathcal{L}f(\boldsymbol{\mathbf{w}})=\sum_{i=1}^{d}\left((b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}-\varepsilon^{\alpha}\frac{\theta_{i}}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}\right)\frac{\partial
f(\boldsymbol{\mathbf{w}})}{\partial
w_{i}}+\varepsilon^{\alpha}\sum_{i=1}^{d}\mathcal{H}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}}),$
(G.22)

where $\mathcal{H}^{\alpha,\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}})$ is
given in (G.4).

It follows that the adjoint operator $\mathcal{L}^{\ast}$ of $\mathcal{L}$ is
given by:

$\displaystyle\mathcal{L}^{\ast}f(\boldsymbol{\mathbf{w}})=-\sum_{i=1}^{d}\frac{\partial}{\partial
w_{i}}\left(\left((b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}-\varepsilon^{\alpha}\frac{\theta_{i}}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}\right)f(\boldsymbol{\mathbf{w}})\right)+\varepsilon^{\alpha}\sum_{i=1}^{d}\mathcal{H}^{\alpha,-\theta_{i}}_{w_{i}}f(\boldsymbol{\mathbf{w}})\,.$
(G.23)

The probability density function $p(\boldsymbol{\mathbf{w}},t)$ of the Lévy-
driven SDE satisfies the Fokker-Planck equation (Schertzer et al., 2001):

$\displaystyle\partial_{t}p(\boldsymbol{\mathbf{w}},t)$
$\displaystyle=\mathcal{L}^{\ast}p(\boldsymbol{\mathbf{w}},t)$
$\displaystyle=-\sum_{i=1}^{d}\frac{\partial}{\partial
w_{i}}\left[\left((b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}-\varepsilon^{\alpha}\theta_{i}\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}\right)p(\boldsymbol{\mathbf{w}},t)\right]+\varepsilon^{\alpha}\sum_{i=1}^{d}\mathcal{H}_{w_{i}}^{\alpha,-\theta_{i}}p(\boldsymbol{\mathbf{w}},t).$
(G.24)

We can compute that

$\displaystyle\sum_{i=1}^{d}\frac{\partial}{\partial
w_{i}}\left[\left((b(\boldsymbol{\mathbf{w}},\alpha,\theta))_{i}-\varepsilon^{\alpha}\theta_{i}\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{(\alpha-1)\Gamma(1-\alpha)}\right)\varphi(\boldsymbol{\mathbf{w}})\right]+\varepsilon^{\alpha}\sum_{i=1}^{d}\mathcal{H}_{w_{i}}^{\alpha,-\theta_{i}}\varphi(\boldsymbol{\mathbf{w}})$
$\displaystyle=-\varepsilon^{\alpha}\sum_{i=1}^{d}\frac{\partial}{\partial
w_{i}}\left[\left(\frac{1}{2}-\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{\int_{w_{i}}^{w_{i}+\xi}\varphi(\boldsymbol{\mathbf{w}}+(y-w_{i})\mathbf{e}_{i})dy-\varphi(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi\right]$
$\displaystyle\quad-\varepsilon^{\alpha}\sum_{i=1}^{d}\frac{\partial}{\partial
w_{i}}\left[\left(\frac{1}{2}+\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{\int_{w_{i}}^{w_{i}-\xi}\varphi(\boldsymbol{\mathbf{w}}+(y-w_{i})\mathbf{e}_{i})dy+\varphi(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi\right]$
$\displaystyle\qquad+\varepsilon^{\alpha}\sum_{i=1}^{d}\mathcal{H}_{w_{i}}^{\alpha,-\theta_{i}}\varphi(\boldsymbol{\mathbf{w}})$
$\displaystyle=-\varepsilon^{\alpha}\sum_{i=1}^{d}\left(\frac{1}{2}-\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{\varphi(\boldsymbol{\mathbf{w}}+\xi\mathbf{e}_{i})-\varphi(\boldsymbol{\mathbf{w}})-\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi$
$\displaystyle\qquad-\varepsilon^{\alpha}\sum_{i=1}^{d}\left(\frac{1}{2}+\frac{\theta_{i}}{2}\right)\frac{1}{\cos(\alpha\pi/2)}\frac{\alpha}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{\varphi(\boldsymbol{\mathbf{w}}-\xi\mathbf{e}_{i})-\varphi(\boldsymbol{\mathbf{w}})+\partial_{w_{i}}\varphi(\boldsymbol{\mathbf{w}})\xi}{\xi^{\alpha+1}}d\xi$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\varepsilon^{\alpha}\sum_{i=1}^{d}\mathcal{H}_{w_{i}}^{\alpha,-\theta_{i}}\varphi(\boldsymbol{\mathbf{w}})=0.$

Hence, we conclude that
$\pi(d\boldsymbol{\mathbf{w}})=\exp(-\varepsilon^{-\alpha}f(\boldsymbol{\mathbf{w}}))d\boldsymbol{\mathbf{w}}/\int_{\mathbb{R}^{d}}\exp(-\varepsilon^{-\alpha}(\boldsymbol{\mathbf{w}}))d\boldsymbol{\mathbf{w}}$
is an invariant distribution of the asymmetric fractional Langevin dynamics
(4.6). Finally, if $b(\boldsymbol{\mathbf{w}},\alpha,\theta)$ is Lipschitz
continuous in $\boldsymbol{\mathbf{w}}$, then $\pi(d\boldsymbol{\mathbf{w}})$
is the unique invariant distribution of (4.6), see e.g. Schertzer et al.
(2001). ∎

###  G.4 Proof of Theorem E.1

Theorem E.1 provides a first-order approximation of the fractional derivative
$\mathcal{D}^{\gamma,-\theta}$ when $d=1$.

Based on the work of Meerschaert & Tadjeran (2004), we will show a first-order
approximation for the asymmetric fractional derivative
$\mathcal{D}^{-\gamma,-\theta}$ when $-1<\gamma<0$ by using the shifted
Grünwald-Letnikov difference operators defined in (E.3) and (E.4). Before we
proceed to the proof of Theorem E.1, we will first present the Fourier
transform property from equations (1) and (12) in Tian et al. (2015).

######  Property G.1 ((Tian et al., 2015)).

Let $-1<\gamma<0$ and $f\in L^{1}(\mathbb{R})$. The Fourier transform of
$\mathcal{I}_{-}^{-\gamma}f$ and $\mathcal{I}_{+}^{-\gamma}f$ satisfy the
following identities:

$\mathcal{F}\left[\mathcal{I}^{-\gamma}_{-}f(w)\right](\zeta)=(i\zeta)^{\gamma}\hat{f}(\zeta),\quad\mathcal{F}\left[\mathcal{I}^{-\gamma}_{+}f(w)\right](\zeta)=(-i\zeta)^{\gamma}\hat{f}(\zeta),$
(G.25)

where $\hat{f}(\zeta)$ denotes the Fourier transform of $f$, such that
$\hat{f}(\zeta)=\int_{-\infty}^{\infty}e^{-i\zeta w}f(w)dw.$

Now, we are ready to prove Theorem E.1.

###### Proof of Theorem E.1.

The main idea for the proof of Theorem E.1 is to use the Fourier transform to
estimate the difference between
$\mathcal{F}\left[\left((1+\theta)\mathcal{A}^{\gamma}_{h,p}+(1-\theta)\mathcal{B}^{\gamma}_{h,q}\right)f(w)\right](\zeta)$
and
$\mathcal{F}\left[\left((1+\theta)\mathcal{I}^{-\gamma}_{-}+(1-\theta)\mathcal{I}^{-\gamma}_{+}\right)f(w)\right](\zeta)$,
and then apply the inverse Fourier transform to complete the proof. By the
linearity of Fourier transforms, we can apply Fourier transform to (E.3) to
obtain

$\displaystyle\mathcal{F}\left[\mathcal{A}^{\gamma}_{h,p}f(w)\right](\zeta)$
$\displaystyle=\frac{1}{h^{\gamma}}\sum_{k=0}^{\infty}(-1)^{k}\binom{-\gamma+k-1}{k}e^{-i\zeta(k-p)h}\hat{f}(\zeta)$
$\displaystyle=\frac{1}{h^{\gamma}}e^{i\zeta ph}\left(1-e^{-i\zeta
h}\right)^{\gamma}\hat{f}(\zeta)$
$\displaystyle=(i\zeta)^{\gamma}\left(\frac{1-e^{-i\zeta h}}{i\zeta
h}\right)^{\gamma}e^{i\zeta ph}\hat{f}(\zeta)$ $\displaystyle=W_{p}(i\zeta
h)(i\zeta)^{\gamma}\hat{f}(\zeta).$ (G.26)

Similarly, we can compute that

$\displaystyle\mathcal{F}\left[\mathcal{B}^{\gamma}_{h,q}f(w)\right](\zeta)$