My sincere thanks also go to all the members of the LTCI lab at Télécom Paris, who provided me an opportunity to join a big team with great researchers, and who gave access to a welcoming and professional research environment.

I would like to offer my special thanks to my family for their unwavering support and belief in me. Without their precious support, I would not have made it through my PhD journey.

Abstract

In this thesis, we are concerned with the stochastic gradient descent (SGD) algorithm, which has been widely used in machine learning due to its computational efficiency and favorable generalization properties. Specifically, we perform theoretical and empirical analysis of the behavior of the stochastic gradient noise (GN) in deep neural networks, which is defined as the difference between the true gradient and the stochastic gradient. Based on these results, we bring an alternative perspective to the existing approaches for investigating SGD. The GN in SGD is often considered to be Gaussian for mathematical convenience. This assumption enables SGD to be studied as a stochastic differential equation (SDE) driven by a Brownian motion. We argue that the Gaussianity assumption might fail to hold in deep learning settings and hence render the Brownian motion-based analyses inappropriate. Inspired by non-Gaussian natural phenomena, we consider the GN in a more general context that suggests that the GN is better approximated by a heavy-tailed α-stable random vector, where tail-index α decides the heavy-tailedness of the distribution. Accordingly, we propose to analyze SGD as a discretization of an SDE driven by a Lévy motion. Firstly, to justify the α-stable assumption, we conduct experiments on common deep learning scenarios and show that in all settings, the GN is highly non-Gaussian and exhibits heavy-tails. We investigate the tail behavior in various network architectures and sizes, loss functions, and datasets.

Secondly, under the heavy-tailed GN assumption, we provide a non-asymptotic analysis for the discrete-time dynamics to converge to the global minimum in terms of suboptimality. This finite-time guarantee is then extended to application in posterior sampling.

Finally, we investigate the metastability nature of the SDE driven by Lévy motion that can then be exploited for clarifying the behavior of SGD, especially in terms of 'preferring wide minima'. More precisely, we provide formal theoretical analysis where we derive explicit conditions for the step-size such that the metastability behavior of SGD, viewed as a discrete-time SDE, is similar to its continuous-time limit. We show that the behaviors of the two systems are indeed similar for small step-sizes and we describe how the error depends on the algorithm and problem parameters. We illustrate our metastability results with simulations on a synthetic model and neural networks. Our results open up a different perspective and shed more light on the view that SGD prefers wide minima.

5

Résumé

Dans cette thèse, nous nous intéressons à l'algorithme du gradient stochastique (SGD), qui est largement utilisé en apprentissage automatique en raison de son efficacité de calcul et de ses propriétés de généralisation. SGD est applicable à un large ensemble de problèmes d'optimisation convexe et non convexe survenant dans l'apprentissage automatique, y compris l'apprentissage profond où ils sont particulièrement réussis. Dans cette étude, nous effectuons une analyse théorique et empirique du comportement du bruit de gradient stochastique (GN), qui est défini comme la différence entre le gradient réel et le gradient stochastique, dans les réseaux de neurones profonds. Sur la base de ces résultats, nous apportons une perspective alternative aux approches existantes pour étudier SGD. Le bruit de gradient stochastique dans SGD est souvent considéré comme gaussien dans le régime des données volumineuses en supposant que le théorème limite central classique (CLT) entre en jeu. Cette hypothèse est souvent faite pour des raisons de commodité mathématique, car elle permet d'étudier SGD comme une équation différentielle stochastique (SDE) pilotée par un mouvement brownien. Nous soutenons que l'hypothèse de la gaussianité pourrait ne pas tenir dans les contextes d'apprentissage profond et donc rendre inappropriées les analyses basées sur le mouvement brownien. Inspiré de phénomènes naturels non gaussiens, nous considérons le bruit de gradient stochastique dans un contexte plus général et invoquons le théorème limite central généralisé, qui suggère que le bruit de gradient stochastique est mieux approché par un vecteur aléatoire à queue lourde α-stable, où l'indice de queue α décide de la lourdeur de la distribution. En conséquence, nous proposons d'analyser SGD comme une discrétisation d'une équation différentielle stochastique pilotée par un mouvement Lévy.

Cette équation différentielle stochastique peut subir des "sauts", qui obligent l'équation différentielle stochastique et son transition de discrétisation de minima étroits vers des minima plus larges, comme le prouve la théorie existante de la métastabilité.

Premièrement, pour justifier l'hypothèse α-stable, nous menons des expériences sur des scénarios communs d'apprentissage en profondeur et montrons que dans tous les contextes, le bruit de gradient stochastique est hautement non gaussien et présente des queues lourdes. Nous étudions le comportement de la queue dans diverses architectures et tailles de réseau, fonctions de perte et ensembles de données.

Deuxièmement, sous l'hypothèse du bruit de gradient stochastique à queue lourde, nous fournissons une analyse non asymptotique pour que la dynamique en temps discret SGD converge vers le minimum global en termes de sous-optimalité. Nos résultats montrent que l'erreur faible sous notre hypothèse augmente plus rapidement que sous l'hypothèse du bruit gaussien, ce qui suggère d'utiliser des tailles de pas plus petits dans notre cas. Cette garantie à temps fini est ensuite étendue à l'application en postérieur échantillonnage.

Enfin, nous étudions la nature de métastabilité de l'équation différentielle stochastique pilotée par le mouvement de Lévy qui peut ensuite être exploitée pour clarifier le comportement de SGD, notamment en termes de "préférence de larges minima". Bien que notre approche apporte une nouvelle perspective pour l'analyse de SGD, elle est limitée en ce sens qu'en raison de la discrétisation temporelle, SGD peut admettre un comportement sensiblement différent de sa limite de temps continu. Dans ce sujet, nous fournissons une analyse théorique formelle où nous dérivons des conditions explicites pour la taille de pas de sorte que le comportement de métastabilité de SGD, considéré comme une équation différentielle stochastique en temps discret, est similaire à sa limite de temps continu. Nous montrons que les comportements des deux systèmes sont en effet similaires pour les petits tailles de pas et nous décrivons comment l'erreur dépend de l'algorithme et des paramètres du problème. Nous illustrons nos résultats de métastabilité par des simulations sur un modèle synthétique et des réseaux de neurones. Nos résultats ouvrent une perspective différente et éclairent davantage l'idée selon laquelle SGD préfère les minima larges.

Chapter 1

Introduction

Machine learning is the study of computer algorithms that instruct a system how to improve from experience. The primary goal is to allow the system to learn automatically to perform predictions, decisions or classifications without being explicitly programmed.

In order to carry out these assigned tasks, machine learning algorithms construct mathematical models using sample data. The more complex the tasks are, the more data is needed for a good performance of the learning process. With the accessibility of large amounts of data (images, speech, video) in recent years, the interest in machine learning has been increasing rapidly. Another even more important reason that makes machine learning successful nowadays is the availability of modern technologies such as parallelprocessing power, high-performance graphics processing units, which can be combined in clusters to significantly reduce the training time for a learning model.

Machine learning systems are present all around us. They are exploited in a broad variety of applications, for example: computer vision [START_REF] Khan | Machine learning in computer vision[END_REF], audio recognition [START_REF] Purwins | Deep learning for audio signal processing[END_REF], [START_REF] Cunningham | Supervised machine learning for audio emotion recognition: Enhancing film sound design using audio features, regression models and artificial neural networks[END_REF], speech recognition Padmanabhan and Johnson Premkumar [2015], natural language processing [START_REF] Shetty | Natural language processing (NLP) for machine learning[END_REF], machine translation [START_REF] Popel | Transforming machine translation: a deep learning system reaches news translation quality comparable to human professionals[END_REF], bioinformatics [START_REF] Larranaga | Machine learning in bioinformatics[END_REF], material inspection [START_REF] Sacco | Machine learning in composites manufacturing: A case study of automated fiber placement inspection[END_REF], medical image analysis [START_REF] Lundervold | An overview of deep learning in medical imaging focusing on MRI[END_REF], drug design [START_REF] Vamathevan | Applications of machine learning in drug discovery and development[END_REF] or board game programs [START_REF] Xenou | Deep reinforcement learning in strategic board game environments[END_REF]. In some cases, they have achieved outcomes as good as and in some circumstances surpassing human expert performance.

Deep learning is a machine learning technique that employs neural networks with a large number of layers for training, using huge amounts of data. Deep neural networks have revolutionized machine learning and have ubiquitous use in many application domains [START_REF] Lecun | Deep learning[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], Hinton et al., 2012a]. There are different types of neural networks that are suitable for different tasks: convolutional neural networks (CNNs) are appropriate for image recognition, while recurrent neural networks (RNNs) are suited for language processing and speech recognition due to their better capabilities to model time series. The architecture of neural networks is also developing through time: researchers design a more efficient type of RNN model called long short-term memory (LSTM), making it run fast enough to be applied in on-demand applications like Google Translate.

In full generality, many key tasks in deep learning reduce to solving the following optimization problem:

w = arg min w∈R d f (w) 1 n n i=1 f (i) (w) (1.1)
where w ∈ R d denotes the weights of the neural network, f : R d → R denotes the loss function that is typically non-convex in w, each f (i) denotes the (instantaneous) loss function that is contributed by the data point i ∈ {1, . . . , n}, and n denotes the total number of data points. Stochastic gradient descent (SGD) is one the most popular approaches for attacking this problem in practice and is based on the following iterative updates:

w k+1 = w k -η∇ fk (w k ) (1.2)
where k ∈ {1, . . . , K} denotes the iteration number, η is the step-size (or the learning rate), and ∇ fk denotes the stochastic gradient at iteration k, that is defined as follows:

∇ fk (w) ∇ fΩ k (w) 1 b i∈Ω k ∇f (i) (w).

(1.3)

Here, Ω k ⊂ {1, . . . , n} is a random subset that is drawn with or without replacement at iteration k, and b = |Ω k | denotes the number of elements in Ω k .

SGD is widely used in deep learning with a great success in its computational efficiency [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF][START_REF] Bottou | The tradeoffs of large scale learning[END_REF][START_REF] Daneshmand | Escaping saddles with stochastic gradients[END_REF]. Beyond efficiency, understanding how SGD performs better than its full batch counterpart in terms of test accuracy remains a major challenge. Even though SGD seems to find perfect training performance at (near-) zero loss solutions on the training landscape (at least in certain regimes [Zhang et al., 2017a[START_REF] Sagun | Explorations on high dimensional landscapes[END_REF][START_REF] Keskar | On largebatch training for deep learning: Generalization gap and sharp minima[END_REF], Geiger et al., 2018]), it appears that the algorithm finds solutions with different properties depending on how it is tuned [START_REF] Sutskever | On the importance of initialization and momentum in deep learning[END_REF][START_REF] Keskar | On largebatch training for deep learning: Generalization gap and sharp minima[END_REF][START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF][START_REF] Hoffer | Train longer, generalize better: closing the generalization gap in large batch training of neural networks[END_REF][START_REF] Masters | Revisiting small batch training for deep neural networks[END_REF][START_REF] Smith | Don't decay the learning rate, increase the batch size[END_REF]. Despite the fact that the impact of SGD on generalization has been studied [START_REF] Advani | High-dimensional dynamics of generalization error in neural networks[END_REF][START_REF] Wu | How SGD selects the global minima in over-parameterized learning: A dynamical stability perspective[END_REF][START_REF] Neyshabur | Exploring generalization in deep learning[END_REF], a satisfactory theory that can explain its success in a way that encompasses such peculiar empirical properties is still lacking.

A popular approach for investigating the behavior of SGD is based on considering SGD as a discretization of a continuous-time process [START_REF] Mandt | A variational analysis of stochastic gradient algorithms[END_REF][START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF], Li et al., 2017, Hu et al., 2017[START_REF] Ye | Stochastic fractional Hamiltonian Monte Carlo[END_REF][START_REF] Chaudhari | Stochastic gradient descent performs variational inference, converges to limit cycles for deep networks[END_REF].

This approach models the stochastic gradient noise, i.e., ∇ fk (w)-∇f (w) as a Gaussian distribution. From this perspective, the SGD recursion can be seen as a first-order Euler-Maruyama discretization of the Langevin dynamics [START_REF] Li | Stochastic modified equations and adaptive stochastic gradient algorithms[END_REF][START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF], Hu et al., 2017], which is often referred to as the Unadjusted Langevin Algorithm (ULA) [START_REF] Roberts | Langevin diffusions and Metropolis-Hastings algorithms[END_REF][START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion: the case of a weakly mean reverting drift[END_REF][START_REF] Durmus | Non-asymptotic convergence analysis for the unadjusted Langevin algorithm[END_REF], Durmus et al., 2016]. Based on this observation, [START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF] focus on the relation between this invariant measure and the algorithm parameters, namely the step-size η and minibatch size. They conclude that the ratio of step-size divided by the batch size is the control parameter that determines the width of the minima found by SGD. Furthermore, they revisit the famous wide minima folklore [START_REF] Hochreiter | Flat minima[END_REF]]: Among the minima found by SGD, the wider it is, the better it performs on the test set.

However, we will identify several fundamental issues with this approach S ¸imşekli et al.

[2019], which show that the gradient noise in SGD can be highly non-Gaussian. Inspired by non-Gaussian natural phenomena, our main approach will be to model the gradient noise by using a more general family of heavy-tailed distributions, which contains the Gaussian distribution as a special case.

This thesis aims to investigate the behavior of SGD by inheriting the idea of stochastic differential equations (SDE). Our general approach can be summarized by the following points:

• By assuming the gradient noise to be Gaussian, SGD is often studied as an approximation of an SDE driven by a Brownian motion. However, by providing empirical evidence on deep learning settings, we show that the gradient noise is highly non-Gaussian. Also according to these experiments, we come up with a better-suited hypothesis for the gradient noise -the α-stable noise. This noise assumption will enable us to use a Lévy-driven SDE as a proxy to SGD. • Next, we investigate the global convergence of SGD by using our new hypothesis. In particular, by using heavy-tailed assumption, we establish a non-asymptotic guarantee for SGD to converge to the global optimum.

• We conclude the thesis by studying the metastable behavior of SGD, which is often considered as an aspect for explaining its generalization property in deep learning settings. Especially, we derive explicit conditions for the step-size such that the metastability of discrete-time SGD is similar to its continuous-time limit.

Summary of the thesis

The main contributions of this thesis are twofold: (i) we perform theoretical and empirical analysis of the tail-index of the stochastic gradient noise in deep neural networks and (ii) based on these results, we bring an alternative perspective to the existing approaches for analyzing SGD and shed more light on the folklore that SGD prefers wide minima by establishing a bridge between SGD and the related theoretical results from statistical physics and stochastic analysis.

The thesis is organized as follows. In Chapter 2 we provide the technical background for Brownian motions, the α-stable distributions as well as the Lévy motions. We briefly introduce the notion of stochastic differential equations (SDEs) and the basic tools such as Wasserstein distance and total variation, which are needed for the theoretical analysis of Fractional Langevin Monte Carlo methods and the metastability analysis.

In Chapter 3, we present the related work that concerns some important aspects in deep learning and diffusion-based Markov chain Monte Carlo.

In Chapter 4 we first formalize the framework in which we analyze SGD by using such SDEs as a proxy. We then describe in Section 4.4 the metastability and first exit time properties of such SDEs and their discretization, and discuss their connection with the wide minima phenomenon. In Section 4.5 we describe our experimental methodology. In Section 4.6 we provide our empirical results which validate our theory. We conduct experiments on the most common deep learning architectures. In particular, we investigate the tail behavior under fully-connected and convolutional models using negative log likelihood (NLL) and linear hinge loss functions on MNIST, CIFAR10, and CIFAR100 datasets. For each configuration, we scale the size of the network and batch size used in SGD and monitor the effect of each of these settings on the tail index α.

In particular, we present results on stability tests (Section 4.6.1), finer-grained layerwise tail-index estimation (sections 4.6.2 and 4.6.3), and an investigation of the relation between the tail-index and the generalization properties of the network (4.6.4).

In Chapter 5, we investigate the global convergence property of SGD for non-convex optimization via a stochastic process, which can be seen as a perturbed version of the gradient descent algorithm with heavy-tailed α-stable noise. In Section 5.2, we state the assumptions that imply the main result presented in Theorem 7, which provides a finitetime guarantee for the discrete-time dynamics of SGD, in terms of suboptimality with respect to the global minimum, as a function of the step-size and the scale parameter. We then describe in Section 5.3 the proof strategy of Theorem 7. Finally in Section 5.5, we extend our results to the case where exact gradients are replaced by stochastic gradients and show that similar results hold in this setting as well. All the proofs of the results in the chapter are given in appendix.

In Chapter 6, we derive explicit conditions for the step-size such that the discretetime SGD (6.3) can inherit the metastability behavior of its continuous-time limit (6.2).

In sections 6.2 and 6.3, we in turn present the main theorem of the chapter and describe the proof overview of the theorem. Our theoretical result is illustrated by numerical experiments on a synthetic model and neural networks in Section 6.4.

Our experiments in Chapter 4 reveal several remarkable results:

• In all our configurations, the stochastic gradient noise turns out to be highly non-Gaussian and possesses a heavy-tailed behavior.

• There is a strong interaction between the network architecture, network size, dataset, and the tail-index, which ultimately determine the dynamics of SGD on the training surface. This observation supports the view that, the geometry of the problem and the dynamics induced by the algorithm cannot be separated from each other.

• In almost all configurations, we observe two distinct phases of SGD throughout iterations. During the first phase, the tail-index rapidly decreases and SGD possesses a clear jump when the tail-index is at its lowest value and causes a sudden jump in the accuracy. This behavior strengthens the view that SGD crosses barriers at the very initial phase.

Our approach also opens up several interesting future directions and open questions, as we discuss in Chapter 7.

List of publications

In this section, we present the list of works that have been realized during this PhD thesis. First, we specify the publications associated directly to the central topic of this document: This work is mainly presented in Chapter 6 of this document. Some other problems have also been dealt with within the time frame of the PhD however they are not included in this document in order to maintain a coherent manuscript.

• S ¸imşekli, U.,
These works are listed below. In this work, we develop an asynchronous-parallel stochastic L-BFGS algorithm for non-convex optimization. The proposed algorithm is suitable for both distributed and shared-memory settings. We provide formal theoretical analysis and show that the proposed algorithm provides a significant speedup over the recently proposed synchronous distributed L-BFGS algorithm.

Chapter 2

Elements of notations and definitions

In this chapter, we introduce basic notations and the definitions of the Wasserstein distance and the total variation. Then, we provide some technical backgrounds for stochastic differential equations driven by Brownian motions as well as by the Lévy motions.

Basic notations

For z > 0, we denote Γ(z) as the gamma function at z:

Γ(z) ∞ 0 x z-1 e -x dx.
For any Borel probability measures µ and ν with domain Ω, their total variation (TV) distance is defined as follows:

µ -ν T V 2 sup A∈B(Ω) |µ(A) -ν(A)|,
where B(Ω) denotes the Borel subsets of Ω.

We use < •, • > to denote the inner product between two vectors, • denotes the Euclidean norm, E ω [•] denotes the expectation with respect to the random variable ω,

and E[•] denotes the expectation with respect to all the random sources. We will use the Wasserstein metric to quantify the distance between two probability measures.

Definition 1 (Wasserstein distance). Let µ and ν be two probability measures. For λ ≥ 1, we define the λ-Wasserstein distance between µ and ν as follows:

W λ (µ, ν) (inf{E V -W λ : V ∼ µ, W ∼ ν}) 1/λ ,
where the infimum is taken over all the couplings of µ and ν (i.e. the joint probability distributions whose marginal distributions are µ and ν).

Notational convenience: In this thesis, depending on certain circumstances, we use f (t) (variable t is in parentheses) as well as f t (variable t is a subscript) to denote a function of t. Accordingly, we use Z(t) as well as Z t to denote a stochastic process.

Stochastic differential equations

Stochastic differential equations (SDEs) play an important role in modeling various phenomena in physics, biology and finance, such as thermal fluctuations, fluid flow or stock prices. An SDE is a differential equation in which one or more coefficients are stochastic processes or stochastic functions of stochastic processes.

SDEs driven by Brownian motions

Consider a (real) differential equation of the form:

dx(t) dt = a(t)x(t), x(0) = x 0 . (2.1)
Suppose that a(t) is a stochastic function and is given as:

a(t) = f (t) + h(t)Z(t),
where f (t) and h(t) are some deterministic functions and Z(t) denotes some stochastic process. Then, equation (2.1) can be written as follows:

dx(t) = f (t)x(t)dt + h(t)x(t)Z(t)dt.
If we further assume that Z(t)dt is in fact the differential form of a Brownian motion, then we obtain the following stochastic differential equation:

dx(t) = f (t)x(t)dt + h(t)x(t)dB(t),
where B(t) denotes the Brownian motion (or Wiener process), which is defined as:

Definition 2 (Brownian motion). A Brownian motion B(t) is a stochastic process satisfying the following properties:

(i) B(0) = 0.
(ii) Independent increments: for every t > 0, s ≥ 0, u ≥ 0 with s ≤ t, the random variables B(t + u) -B(t), B(s) are independent.

(iii) Gaussian increments: for all 0 ≤ u ≤ t, the random variable B(t + u) -B(t) has the same distribution as Gaussian random variable N (0, u).

(iv) Continuous paths: B(t) is continuous in t.

Notational convenience:

In this thesis, we also use B t to denote the Brownian motion.

More generally, an SDE is given as:

dX(t) = f (t, X(t))dt + h(t, X(t))dB(t).
(2.2) In the following, we consider examples of SDEs driven by Brownian motions.

Examples:

1. Ornstein-Uhlenbeck process. It is defined by the following SDE:

dX(t) = -θX(t)dt + σdB(t),
where θ, σ > 0. The solution of this SDE can be written as follows.

X(t) = X(0) exp(-θt) + t 0 σ exp(-θ(t -s))dB(s).
2. Langevin equation. The Langevin diffusion is described by the following SDE:

dX(t) = -∇f (X(t))dt + √ 2dB(t),
where f is a smooth function. The above Langevin diffusion possesses an invariant measure π, whose density is proportional to exp(-f (x)).

In order to illustrate the random behaviors of different phenomena, it is possible to use other types of stochastic processes, such as jump processes like Lévy process, instead of Brownian motion. In the next section, we introduce the notion of stable distributions and the SDEs driven by Lévy motions.

Stable distributions and Lévy motions

The central limit theorem (CLT) states that the sum of i.i.d. random variables with a finite second moment converges to a normal distribution if the number of summands grows. However, if the variables have heavy-tails, the second moment may not exist. For instance, if their density p(x) has a power-law tail decreasing as 1/|x| α+1 (Figure 2.1) where 0 < α < 2; only r-th moment exists with r < α. In this case, the generalized central limit theorem (GCLT) says that the sum of such variables will converge to a distribution called the α-stable distribution instead as the number of summands grows (see e.g. [START_REF] Fischer | A history of the central limit theorem: From classical to modern probability theory[END_REF]). In this work, we focus on the centered symmetric α-stable (SαS) distribution, which is a special case of α-stable distributions that are symmetric around the origin.

We can view the SαS distribution as a heavy-tailed generalization of a centered Gaussian distribution. The SαS distributions are defined through their characteristic function:

Definition 3 (Symmetric α-stable random variables). The α-stable distribution appears as the limiting distribution in the generalized CLT Samorodnitsky and Taqqu [1994]. A

scalar random variable X ∈ R is called symmetric α-stable if its characteristic function has the following form:

E[e iωX ] = exp(-σ|ω| α )
where α ∈ (0, 2] and σ > 0. We denote X ∼ SαS(σ).

Even though their probability density function does not admit a closed-form formula in general except in special cases, their density decays with a power law tail like 1/|x| α+1

where α ∈ (0, 2] is called the tail-index which determines the behavior of the distribution:

as α gets smaller, the distribution has a heavier tail. In fact, the parameter α also determines the moments: when α < 2, E[|X| r ] < ∞ if and only if r < α; implying X has infinite variance when α = 2. The parameter σ ∈ R + is the scale parameter and controls the spread of X around 0. We recover the Gaussian distribution N (0, 2σ 2 ) as a special case when α = 2 and the Cauchy distribution when α = 1.

For the scalar case, α-stable Lévy motion is defined as follows for α ∈ (0, 2] ( [2015]):

Definition 4 (Symmetric α-stable Lévy motion). A scalar symmetric α-stable Lévy motion L α (t), with 0 < α ≤ 2, is a stochastic process satisfying the following properties:

(i) L α (0) = 0, almost surely.

(ii) Independent increments: for 0 ≤ t 1 < . . . < t n , the random variables L α (t 2 ) -L α (t 1 ),..., L α (t n ) -L α (t n-1 ) are independent.

(iii) Stationary increments: for all 0 ≤ s < t, the random variables L α (t) -L α (s) and L α (t -s) have the same distribution as SαS((t -s) 1/α ).

(iv) Continuity in probability: for any δ > 0 and s ≥ 0,

P(|L α (s) -L α (t)| > δ) → 0, as t → s.
Notational convenience: In this thesis, we also use L t to denote the Lévy motion.

By using this definition, one can define the Lévy-driven SDEs by replacing the Brownian motion in SDE (2.2) by a Lévy motion:

dX(t) = f (t, X(t))dt + h(t, X(t))dL(t). (2.3)
To conclude this chapter, let us define the notion of invariant measure for equation (2.3).

Definition 5. Let q(X, t) be the probability density function of X(t). Then q(X, t) is called invariant measure for (2.3) if ∂ t q(X, t) = 0.

In the next chapter, we will give an example of Lévy-driven SDE, which is called the Fractional Langevin Monte Carlo, along with its invariant measure.

Chapter 3

Related work

In this chapter, we present some lines of works and contributions that have connections with our work: Important aspects of deep learning theory and diffusion-based Markov

Chain Monte Carlo algorithms.

Deep learning theory

Despite the success of deep neural networks in many important machine learning applications [START_REF] Carleo | Machine learning and the physical sciences[END_REF], [START_REF] Barbier | Optimal errors and phase transitions in high-dimensional generalized linear models[END_REF], [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF], a satisfactory theory of deep learning is still lacking. The theoretical understanding of deep learning includes the following areas [START_REF] Poggio | Theoretical issues in deep networks: Approximation, optimization and generalization[END_REF]: 1) representation power of deep networks, 2) optimization of the empirical risk, 3) generalization properties of gradient descent techniques.

Representation power of deep networks

We start with the representation power of deep networks. Even though both deep and shallow networks have a universal property [START_REF] Poggio | Theoretical issues in deep networks: Approximation, optimization and generalization[END_REF], that is, they can approximate any continuous function of finite variables on a compact domain, using deep networks for approximation can achieve much better performance than using shallow networks. In contrast to shallow networks, deep networks can avoid the curse of dimensionality on certain classes of problems. For instance, with the same degree of approximation, one can represent compositional functions with much smaller number of parameters for the deep networks than for the shallow networks [START_REF] Poggio | Theoretical issues in deep networks: Approximation, optimization and generalization[END_REF].

That is to say, the expressive power of deep neural networks can be considered from the view of the depth of a network. Recent works [START_REF] Cohen | On the expressive power of deep learning: A tensor analysis[END_REF], [START_REF] Telgarsky | Benefits of depth in neural networks[END_REF] show the depth-efficiency of deep neural networks by proving that there exist classes of deep networks that cannot be approximated by any shallow network whose width is smaller than or equal to an exponential bound. In [START_REF] Lu | The expressive power of neural networks: A view from the width[END_REF], the authors deal with the dual problem on the width-efficiency of neural networks. By using width-bounded ReLU (rectified linear unit) networks, they show that there exist classes of wide networks that cannot be represented by any narrow network whose depth is smaller than or equal to an polynomial bound. These results suggest that increasing the depth may be more effective than the width for increasing expressive power of neural networks.

Recently, several new approaches are proposed to understand the representation power of deep neural networks. In [START_REF] Raghu | On the expressive power of deep neural networks[END_REF], the authors introduce a set of measures of expressivity, which determines how the output of a network changes as the input moves along an one-dimensional trajectory. By using these measures, they prove that the complexity of the computed function grows exponentially with depth.

Furthermore, they find that the weights of trained networks are not equal in the sense that the networks tend to be more sensitive to the weights in initial layers. This result suggests that optimizing these weights is particularly important. Another approach to investigate the representation problem is based on the mathematical theory of quiver representation [START_REF] Armenta | The representation theory of neural networks[END_REF], which is used for exploring the combinatorial and algebraic nature of neural networks. In [START_REF] Armenta | The representation theory of neural networks[END_REF], the authors establish an explicit connection between neural networks and quiver representations and show that quiver representations are able to adapt common concepts of neural networks such as fully-connected layers, convolution operations, residual connections, batch normalization and pooling operations. The quiver representations also help understand how neural networks create representations from the data. The representation power alters when dealing with different types of neural networks. In [START_REF] Dehmamy | Understanding the representation power of graph neural networks in learning graph topology[END_REF], the authors study the expressiveness of graph convolutional networks (GCNs), which are capable of distinguishing graphs from different graph generation models. They conclude that GCNs with different propagation rules could improve the representation power significantly.

Optimization techniques

Deep neural networks have proven successful in a wide variety of applications. However, as neural networks grow deeper and datasets become bigger, training these networks becomes more difficult. Therefore, a large number of optimization techniques have been developed over the past years to improve the training performance of these networks.

Although traditional gradient-based methods using a full-gradient approach may be effective for small-scale learning problems, stochastic-gradient-based methods (SG), first proposed by [START_REF] Robbins | A stochastic approximation method[END_REF], are the dominant techniques for large-scale learning problems. Stochastic gradient methods gain a great success in the study of perceptual tasks such as speech or image recognition [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF], in which, due to the use of deep neural networks, highly nonlinear and non-convex problems are involved. Over the years, a large variety of stochastic-gradient-based algorithms have been proposed. These algorithms can be classified into two categories: noise reduction methods and second-order methods.

In optimization problems, methods with noise reduction capabilities have been developed to improve the accuracy of the outcomes as well as the convergence performance.

The first type of noise reduction is dynamic sampling [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF], in which one attempts to reduce the stochastic gradient noise by gradually increasing the mini-batch size, thus improving the accuracy of the gradient estimates as the algorithm proceeds.

On the other hand, gradient aggregation methods, such as stochastic variance reduced gradient (SVRG), stochastic average gradient (SAG) and SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF]), achieve the noise reduction by improving the quality of the search directions using the information from computed gradient estimates in previous iterations [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF].

These noise reduction techniques allow the optimization process to attain a linear convergence rate to the optimal value using constant step-size [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF].

Another important family of stochastic-gradient-based methods, known as secondorder methods, consists of algorithms that aim at addressing the negative effects of high non-linearity and ill-conditioning of the objective function with the help of second-order information [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF] such as second derivative, Hessian matrix, second-order Taylor series, Fisher information matrix, etc. Popular representatives of these methods are inexact Newton and quasi-Newton methods, (generalized) Gauss-Newton methods [START_REF] Bertsekas | Incremental least squares methods and the extended Kalman filter[END_REF], [START_REF] Schraudolph | Fast curvature matrix-vector products for second-order gradient descent[END_REF], the natural gradient method [START_REF] Amari | Natural gradient works efficiently in learning[END_REF],

and scaled gradient iterations [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF], [START_REF] Tieleman | Divide the gradient by a running average of its recent magnitude. coursera: Neural networks for machine learning[END_REF].

Some other well-known adaptive SG methods are not well classified within the two categories: noise reduction and second-order methods, yet show great potential in theoretical and/or practical problems. We start with SG methods with momentum, which aim at improving learning performance by solving the problem of poor conditioning on the Hessian matrix and the variance in stochastic gradient [START_REF] Soydaner | A comparison of optimization algorithms for deep learning[END_REF]. The idea is to take a moving average between the stochastic gradient of the current and previous iterations according to a momentum parameter [START_REF] Alpaydin | Introduction to machine learning[END_REF]. A variant of the standard momentum is Nesterov momentum [START_REF] Sutskever | On the importance of initialization and momentum in deep learning[END_REF], in which the gradient of the objective function is measured slightly ahead of the current update, in the direction of the momentum. Next, we move to AdaGrad [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF], which is an optimization algorithm that adapts the learning rates of the model parameters, namely, at each iteration the parameters with larger (smaller) partial derivative are assigned with larger (smaller) learning rates. The learning rates are computed using all the squared values of the gradient from the previous iteration. One drawback of Ada-Grad is that the continual decay of learning rates may cause the algorithm to stop too early when training neural networks, as the learning rates may become infinitesimally small. Therefore, AdaDelta [START_REF] Zeiler | Adadelta: an adaptive learning rate method[END_REF] is proposed to address this issue. By using the window of fixed size w of past gradients instead of the full-gradients, AdaDelta is able to continue running even after many iterations have been progressed. Another modified version of AdaGrad is root mean square propagation (RMSProp) Hinton et al. [2012b],

which is designed to perform better in non-convex setting. In RMSProp, the gradient accumulation is replaced by an exponentially weighted moving average, which removes the excessive decay of learning rates of AdaGrad and makes the algorithm converge quickly after finding a convex bowl [START_REF] Goodfellow | Deep learning[END_REF]. One of the most popular optimization algorithm in deep learning is Adam, which incorporate the advantages of AdaGrad and RMSProp [START_REF] Kingma | A method for stochastic optimization[END_REF]. The algorithm adapts learning rates for model parameters from the information of the first and second moments of the gradients [START_REF] Soydaner | A comparison of optimization algorithms for deep learning[END_REF]. Besides computational efficiency and small-memory requirement,

Adam is known to be suitable for non-stationary objectives and systems with highly noisy and sparse gradients [START_REF] Kingma | A method for stochastic optimization[END_REF]. Some well-known modified/improved algorithm based on Adam are AdaMax, Nadam [START_REF] Dozat | Incorporating nesterov momentum into adam[END_REF] and AMSGrad Reddi et al.

[2019], which are proposed to improve the speed of convergence while maintaining the benefits of Adam. While having great advantages in terms of computation efficiency and convergence speed, Adam and its modified algorithms have been found to generalize poorly compared to SGD [START_REF] Keskar | Improving generalization performance by switching from Adam to SGD[END_REF], that suggests using a hybrid strategy that combines an adaptive method and SGD [START_REF] Keskar | Improving generalization performance by switching from Adam to SGD[END_REF]. [START_REF] Roberts | Langevin diffusions and Metropolis-Hastings algorithms[END_REF]:

Generalization and SGD

w k+1 ULA = w k ULA -η∇f (w k ULA ) + 2η β Z k+1 , (3.2)
where k ∈ N + denotes the iterations, η denotes the step-size, and (Z n ) n is a sequence of independent and identically-distributed (i.i.d.) standard Gaussian random variables.

When β = 1, we obtain the classical ULA, which is mainly used for Bayesian posterior sampling. Theoretical properties of the classical ULA have been extensively studied [START_REF] Roberts | Langevin diffusions and Metropolis-Hastings algorithms[END_REF], [START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion: the case of a weakly mean reverting drift[END_REF], Durmus andMoulines [2015, 2016], Dalalyan [2017b].

When β 1, the algorithm is called tempered and becomes more suitable for optimization. Indeed, one can observe that the noise term Z k in (3.2) becomes less dominant, and the overall algorithm can be seen as a 'perturbed' version of the gradient descent (GD) algorithm. The connection between ULA and GD has been recently established in Dalalyan [2017a] for strongly convex f . Moreover, [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF] and [START_REF] Xu | Global convergence of langevin dynamics based algorithms for nonconvex optimization[END_REF] proved non-asymptotic guarantees for this perturbed scheme 1 . Their results

showed that, even in non-convex settings, the algorithm is guaranteed to escape from local minima and converge near the global minimizer. These results were extended in Zhang et al. [2017b] and [START_REF] Tzen | Local optimality and generalization guarantees for the Langevin algorithm via empirical metastability[END_REF], which showed that the iterates converge near a local minimum in polynomial time and stay there for an exponential time. Recently, the guarantees for ULA were further extended to second-order Langevin dynamics Gao et al. [2018b,a].

Fractional Langevin Monte Carlo

Another line of research has extended Langevin Monte Carlo by replacing the Brownian motion with a motion which can incur 'jumps' (i.e. discontinuities), such as the αstable Lévy Motion (see Figure 2.1) S ¸imşekli [2017], [START_REF] Ye | Stochastic fractional Hamiltonian Monte Carlo[END_REF]. Coined under the name of Fractional Langevin Monte Carlo (FLMC) methods, these approaches are motivated by the statistical physics origins of the Langevin equation (3.1). In such a context, the Langevin equation aims to model the position of a small particle that is under the influence of a force, which has a deterministic and a stochastic part. If we assume that the stochastic part of this force is a sum of many i.i.d. random variables with finite variance, then by the central limit theorem (CLT), we can assume that their sum follows a Gaussian distribution, which justifies the Brownian motion in (3.1).

The main idea in FLMC is to relax the finite variance assumption and allow the random pulses to have infinite variance. In such a case, the classical CLT will not hold; however, the extended CLT [START_REF] Lévy | Théorie de l'addition des variables aléatoires[END_REF] will still be valid: the law of the sum of the pulses converges to an α-stable distribution, a family of 'heavy-tailed' distributions that

1 The results given in [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF] are more general in the sense that they are proved for the Stochastic Gradient Langevin Dynamics (SGLD) algorithm [START_REF] Welling | Bayesian learning via stochastic gradient Langevin dynamics[END_REF], which is obtained by replacing the gradients in (3.2) with stochastic gradients.

contains the Gaussian distribution as a special case. Then, by using a similar argument to the previous case, we can replace the Brownian motion with the α-stable Lévy Motion [START_REF] Yanovsky | Lévy anomalous diffusion and fractional Fokker-Planck equation[END_REF], whose increments are α-stable distributed.

Based on an SDE driven by an α-stable Lévy Motion, S ¸imşekli [2017] proposed the following iterative scheme that is referred to as Fractional Langevin Algorithm (FLA):

w k+1 FLA = w k FLA -ηc α ∇f (w k FLA ) + η β 1 α S k+1 , (3.3)
where α ∈ (1, 2] is called the characteristic index, c α is a known constant, and {S k } k∈N+ is a sequence of α-stable distributed random variables. Recently, [START_REF] Ye | Stochastic fractional Hamiltonian Monte Carlo[END_REF] extended FLA to Hamiltonian dynamics. The experimental results in S ¸imşekli [2017] and [START_REF] Ye | Stochastic fractional Hamiltonian Monte Carlo[END_REF] showed that the use of the heavy-tailed increments can provide advantages in multi-modal settings, robustness to algorithm parameters. [START_REF] Ye | Stochastic fractional Hamiltonian Monte Carlo[END_REF] further illustrated that in an optimization context their algorithm achieves better generalization in deep neural networks.

The FLA algorithm is based on a Lévy-driven SDE, that is defined as follows:

dw t = Ψ(w t-, α)dt + (1/β) 1/α dL α t (3.4)
where w t-denotes the left limit of the process at time t and L α t denotes the ddimensional Lévy motion whose components are independent α-stable Lévy motions in R.

We have the following result for FLMC:

Theorem 1 (S ¸imşekli [2017]). Consider the SDE (3.4) in the case d = 1, β = 1, and α ∈ (1, 2], where the drift Ψ is defined as follows:

Ψ(x, α) - D α-2 φ(x) ∂f (x) ∂x φ(x) . (3.5)
where φ(x) exp(-βf (x)) and D denotes the fractional Riesz derivative and is defined as follows for a function u:

D γ u(x) F -1 {|ω| γ û(ω)}.
Here, F denotes the Fourier transform and û F(u).

Let π be a random measure whose density is φ(x) (up to a multiplicative factor).

Then, π is an invariant measure of the Markov process (w t ) t≥0 that is a solution of the SDE given by (3.4).

This theorem states that if the drift (3.5) can be computed, then the sample paths of (3.4) can be considered as samples drawn from π. However, computing (3.5) is in general not tractable, therefore one needs to approximate it for computational purposes.

If we use the alternative definition of the Riesz derivative given by [START_REF] Ortigueira | Riesz potential operators and inverses via fractional centred derivatives[END_REF], we can approximate the drift as follows S ¸imşekli [2017], [START_REF] Ye | Stochastic fractional Hamiltonian Monte Carlo[END_REF]:

- D α-2 φ(x) ∂f (x) ∂x φ(x) ≈ -c α ∂f (x) ∂x ,
where c α Γ(α -1)/Γ(α/2) 2 and Γ denotes the Gamma function. With this choice of approximation, in the d-dimensional case we obtain FLA, as given in (3.3). We can observe that, when α = 2, (3.4) becomes the Langevin equation (3.1) and FLA becomes ULA.

Chapter 4

Heavy-tailed behavior in stochastic gradient descent

The gradient noise in SGD is often considered to be Gaussian for mathematical convenience. This assumption enables SGD to be studied as a stochastic differential equation

Overview

A popular approach for investigating the behavior of SGD is based on considering SGD as a discretization of a continuous-time process [START_REF] Mandt | A variational analysis of stochastic gradient algorithms[END_REF][START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF], Li et al., 2017, Hu et al., 2017[START_REF] Ye | Stochastic fractional Hamiltonian Monte Carlo[END_REF][START_REF] Chaudhari | Stochastic gradient descent performs variational inference, converges to limit cycles for deep networks[END_REF]. This approach models the stochastic gradient noise as a Gaussian distribution, i.e. U k (w)

∇ fk (w) -∇f (w) satisfies U k (w) ∼ N (0, σ 2 I), (4.1)
where N denotes the multivariate (Gaussian) normal distribution and I denotes the identity matrix of appropriate size. 1 The rationale behind this assumption is that, if the size of the minibatch b is large enough, then we can invoke the Central Limit Theorem (CLT) and assume that the distribution of U k is approximately Gaussian. Then, under this assumption, SGD (1.2) can be written as follows:

w k+1 = w k -η∇f (w k ) + √ η ησ 2 Z k , (4.2)
where Z k denotes a standard normal random vector in R d . If we further assume that η is small enough, then the continuous-time analogue of the discrete-time process (4.2) is the following stochastic differential equation (SDE):

dw t = -∇f (w t )dt + ησ 2 dB t , (4.3)
where B t denotes the standard Brownian motion. This SDE is a variant of the wellknown Langevin diffusion and under mild regularity assumptions on f , one can show that the Markov process (w t ) t≥0 is ergodic with its unique invariant measure, whose density is proportional to exp(-f (w)/(ησ 2 )) for any η > 0 [START_REF] Roberts | Langevin diffusions and Metropolis-Hastings algorithms[END_REF].

From this perspective, the SGD recursion in (4.2) can be seen as a first-order Euler-Maruyama discretization of the Langevin dynamics (see also [START_REF] Li | Stochastic modified equations and adaptive stochastic gradient algorithms[END_REF][START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF], Hu et al., 2017]), which is often referred to as the Unadjusted Langevin Algorithm (ULA) [START_REF] Roberts | Langevin diffusions and Metropolis-Hastings algorithms[END_REF][START_REF] Lamberton | Recursive computation of the invariant distribution of a diffusion: the case of a weakly mean reverting drift[END_REF][START_REF] Durmus | Non-asymptotic convergence analysis for the unadjusted Langevin algorithm[END_REF], Durmus et al., 2016]. Based on this observation, [START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF] focused on the relation between this invariant measure and the algorithm parameters, namely the step-size η and minibatch size, as a function of σ 2 . They concluded that the ratio of step-size divided by the batch size is the control parameter that determines the width of the minima found by SGD. Furthermore, they revisit the famous wide minima folklore [START_REF] Hochreiter | Flat minima[END_REF]]: Among the minima found by SGD, the wider it is, the better it performs on the test set. This is visualized in Figure 4.1, where the local minimum on the right lies on a wider valley with width m 2 compared to the local minimum on the left with width m 1 lying in a sharp valley of depth H. However, there are several fundamental issues with this approach, which we will explain below. We first illustrate a typical mismatch between the Gaussianity assumption and the empirical behavior of the stochastic gradient noise in terms of the long term behavior.

In Figure 4.2, we plot the histogram of the norms of the stochastic gradient noise at the first and the last iterations that are computed using a convolutional neural network (AlexNet) in an image classification problem on the CIFAR10 dataset and compare it to the histogram of the norms of Gaussian random vectors. It can be clearly observed that, even though the shape of the histogram corresponding to gradients resembles the one of the Gaussian vectors at the first iteration, throughout training, it drifts apart from the Gaussian and exhibits a heavy-tailed behavior.

In addition to the empirical observations, the Gaussianity assumption also yields some theoretical issues. The first issue with this assumption is that the current SDE analyses of SGD are based on the invariant measure of the SDE, which implicitly assumes that sufficiently many iterations have been taken to converge to that measure.

Recent results on ULA [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF], Xu et al., 2018] have shown that, the required number of iterations to achieve the invariant measure often grows exponentially with the dimension d. This result contradicts with the current practice: considering the large size of the neural networks and limited computational budget, only a limited number of iterations -which is much smaller than exp(O(d)) -can be taken. This conflict becomes clearer in the light of the recent works that studied the local behavior of ULA [START_REF] Tzen | Local optimality and generalization guarantees for the Langevin algorithm via empirical metastability[END_REF], Zhang et al., 2017b]. These studies showed that ULA will get close to the nearest local optimum in polynomial time; however, the required amount of time for escaping from that local optimum increases exponentially with the dimension.

Therefore, the phenomenon that SGD prefers wide minima within a considerably small number of iterations cannot be explained using the asymptotic distribution of the SDE given in (4.3).

The second issue is related to the local behavior of the process and becomes clear when we consider the metastability analysis of Brownian motion-driven SDEs. These studies [START_REF] Freidlin | Random perturbations[END_REF][START_REF] Bovier | Metastability in reversible diffusion processes i: Sharp asymptotics for capacities and exit times[END_REF], Imkeller et al., 2010b] consider the case where w 0 is initialized in a quadratic basin and then analyze the minimum time t such that w t is outside that basin. They show that this so-called first exit time depends exponentially on the height of the basin; however, this dependency is only polynomial with the width of the basin. These theoretical results directly contradict with the wide minima phenomenon: even if the height of a basin is slightly larger, the exit-time from this basin will be dominated by its height, which implies that the process would stay longer in (or in other words, 'prefer') deeper minima as opposed to wider minima. The reason why the exit-time is dominated by the height is due to the continuity of the Brownian motion, which is in fact a direct consequence of the Gaussian noise assumption.

A final remark on the issues of this approach is the observation that landscape is flat at the bottom regardless of the batch size used in SGD [START_REF] Sagun | Empirical analysis of the Hessian of over-parametrized neural networks[END_REF]. In particular, the spectrum of the Hessian at a near critical point with close to zero loss value has many near zero eigenvalues. Therefore, local curvature measures that are used as a proxy for measuring the width of a basin can be misleading. Such measures usually correlate with the magnitudes of large eigenvalues of the Hessian which are few [START_REF] Keskar | On largebatch training for deep learning: Generalization gap and sharp minima[END_REF][START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF]. Besides, during the dynamics of SGD it has been observed that the algorithm does not cross barriers except perhaps at the very initial phase [START_REF] Xing | A walk with SGD[END_REF], Baity-Jesi et al., 2018]. Such dependence of width on an essentially-flat landscape combined with the lack of explicit barrier crossing during the SGD descent forces us to rethink the analysis of basin hopping under a noisy dynamics.

In this study, we aim at addressing these contradictions and come up with an arguably better-suited hypothesis for the stochastic gradient noise that has more pertinent theoretical implications for the phenomena associated with SGD.

Proposed framework

We go back to (1.3) and (4.1) and reconsider the application of CLT. This classical CLT assumes that U k is a sum of many independent and identically distributed (i.i.d.) random vectors, whose covariance matrix exists and is invertible, and then it states that the law of U k converges to a Gaussian distribution, which then paves the way for (4.2). Even though the finite-variance assumption seems natural and intuitive at the first sight, it turns out that in many domains, such as turbulent motions [START_REF] Weeks | Observation of anomalous diffusion and Lévy flights[END_REF]), oceanic fluid flows [START_REF] Woyczyński | Lévy processes in the physical sciences[END_REF]), finance [START_REF] Mandelbrot | Fractals and Scaling in Finance: Discontinuity, Concentration, Risk. Selecta[END_REF]), biological evolution [START_REF] Jourdain | Lévy flights in evolutionary ecology[END_REF]), audio signals [START_REF] Liutkus | Generalized Wiener filtering with fractional power spectrograms[END_REF], S ¸imşekli et al. [2015], [START_REF] Leglaive | Alpha-stable multichannel audio source separation[END_REF], S ¸imşekli et al. [2018]), brain signals [START_REF] Jas | Learning the morphology of brain signals using alpha-stable convolutional sparse coding[END_REF]), the assumption might fail to hold (see [START_REF] Duan | An Introduction to Stochastic Dynamics[END_REF] for more examples). In such cases, the classical CLT along with the Gaussian approximation will no longer hold. While this might seem daunting, fortunately, one can prove a generalized CLT and show that the law of the sum of these i.i.d. variables with infinite variance still converges to a family of heavy-tailed distributions that is called the α-stable distribution [START_REF] Lévy | Théorie de l'addition des variables aléatoires[END_REF].

As we detailed in Section 2.2.2, these distributions are parametrized by their tail-index α ∈ (0, 2] and they coincide with the Gaussian distribution when α = 2.

In this study, we relax the finite-variance assumption on the stochastic gradient noise and by invoking the generalized CLT: we assume that U k follows an α-stable distribution, as hinted in Figure 4.2(d). By following a similar rationale to (4.2) and ( 4.3), we reformulate SGD with this new assumption and consider its continuous-time limit for small step-sizes. Since the noise might not be Gaussian anymore (i.e. when α = 2), the use of the Brownian motion would not be appropriate in this case and we need to replace it with the α-stable Lévy motion, whose increments have an α-stable distribution [START_REF] Yanovsky | Lévy anomalous diffusion and fractional Fokker-Planck equation[END_REF]). Due to the heavy-tailed nature of α-stable distribution, the Lévy motion might incur large discontinuous jumps and therefore exhibits a fundamentally different behavior than the Brownian motion, whose paths are on the contrary almost surely continuous. The discontinuities also reflect in the metastability properties of Lévy-driven SDEs, which indicate that, as soon as α < 2, the first exit time from a basin does not depend on its height; on the contrary, it directly depends on its width and the tail-index α Imkeller and [START_REF] Imkeller | First exit times of SDEs driven by stable Lévy processes[END_REF], Imkeller et al. [2010b,a]. Informally, this implies that the process will escape from narrow minima -no matter how deep they are -and stay longer in wide minima. Besides, as α gets smaller, the probability for the dynamic to jump into a wide basin will increase. Therefore, if the α-stable assumption on the stochastic gradient noise holds, then the existing metastability results automatically provide strong theoretical insights for illuminating the behavior of SGD.

SGD as a Lévy-driven SDEs

Following the above argument, a more general assumption on the stochastic gradient noise (cf. (4.1)) can be given by:

[U k (w)] i ∼ SαS i (σ i (w)), ∀i = 1, . . . , n (4.4)
where [v] i denotes the i'th component of a vector v, and SαS i distributed with α i (w).

Clearly, this assumption is way too general to offer reasonable theoretical treatment.

We will resort to several simplifications: (1) We assume that each coordinate of U k is SαS distributed with the same σ which depends on the state w. Here, this dependency is not crucial since we are mainly interested in the tail-index α, which can be estimated independently from the scale parameter (Section 4.5.1). Therefore, we will simply denote σ(w) as σ for clarity.

(2) We further assume that each coordinate of U k is SαS distributed with the same α independent of the state w. We will demonstrate the state independence at later stages of SGD experimentally in Section 4.6.3, however, imposing the coordinate dependence is a much harder challenge which will be addressed in the section devoted for open problems (Chapter 7).

By using the assumption (4.4), we can rewrite the SGD recursion as follows [S ¸imşekli, 2017]:

w k+1 = w k -η∇f (w k ) + η 1/α η α-1 α σ S k , (4.5)
where

S k ∈ R d is a random vector such that [S k ] i ∼ SαS(1). If the step-size η is small
enough, then we can consider the continuous-time limit of this discrete-time process, which is expressed in the following SDE driven by an α-stable Lévy process:

dw t = -∇f (w t )dt + η (α-1)/α σ dL α t , (4.6)
where L α t denotes the d-dimensional α-stable Lévy motion with independent components. In other words, each component of L α t is an independent α-stable Lévy motion in R. It is easy to check that the noise term in (4.5) is obtained by integrating dL α t from kη to (k + 1)η. When α = 2, L α t coincides with a scaled version of Brownian motion, √ 2B t . SαS and L α t are illustrated in Figure 2.1. Stochastic processes based on Lévy motion such as (4.6) and their mathematical properties have also been studied in the literature, we refer the reader to [START_REF] Tankov | Financial modelling with jump processes[END_REF], [START_REF] Øksendal | Applied stochastic control of jump diffusions[END_REF] for details.

First exit time and metastability properties

Consider a basin in which a particle is initialized and undergoes fluctuations continually, the particle persists in the basin for a long time before exiting it by the influence of fluctuations. This relative instability phenomenon is described by the term 'metastability'. More formally, the metastability studies consider the case where w 0 is initialized in a basin and analyze the minimum time t such that w t exits that basin.

We start by reviewing known metastability properties of the α-stable Lévy process (4.6) from the literature. We will also focus on the first exit time which is, roughly speaking, the average time it takes for the process to exit a neighborhood of a local minima (a quantity we define formally later in (4.12)).

For clarity of the presentation and notational simplicity we focus on the scalar case and consider the SDE (4.6) in R (i.e. d = 1). Multidimensional generalizations of the metastability results presented in this work can be found in [Imkeller et al., 2010a] and will be summarized at the end of this section. We rewrite (4.6) as follows:

dw ε t = -f (w ε t )dt + εdL α t (4.7)
for t ≥ 0, started from the initial point w 0 ∈ R, where L α t is the α-stable Lévy process, ε ≥ 0 is the noise level and f is a non-convex objective with r ≥ 2 local minima.

We denote the derivative of f by f . When ε = 0, we recover the gradient descent dynamics in continuous time: dw 0 t = -f (w 0 t )dt, where the local minima are the stable points of this differential equation. However, as soon as ε > 0, these states become 'metastable', meaning that there is a positive probability for w ε t to transition from one basin to another. However, the time required for transitioning to another basin strongly depends on the characteristics of the injected noise. The two most important cases are α = 2 and α < 2. When α = 2, (i.e. the Gaussianity assumption) the process (w ε t ) t≥0 is continuous, which requires it to 'climb' the basin all the way up, in order to be able to transition to another basin. This fact makes the transition-time depend on the height of the basin. On the contrary, when α < 2, the process can incur discontinuities and does not need to cross the boundaries of the basin in order to transition to another one, since it can directly jump. This property is called the 'transition phenomenon' [START_REF] Duan | An Introduction to Stochastic Dynamics[END_REF] and makes the transition-time mostly depend on the width of the basin. In the rest of the section, we will formalize these explanations.

Gradient-like flows driven by Brownian motion and weak error for their discretization are well studied from a theoretical standpoint (see e.g. [Li et al., 2017, Mertikopoulos and[START_REF] Mertikopoulos | On the convergence of gradient-like flows with noisy gradient input[END_REF]), however their Lévy-driven analogue (4.7) and the discrete-time versions [START_REF] Burghoff | Spectral analysis for a discrete metastable system driven by Lévy flights[END_REF] are relatively less studied. Under some assumptions on the objective f , it is known that the process (4.7) admits a stationary density [START_REF] Samorodnitsky | Tails of solutions of certain nonlinear stochastic differential equations driven by heavy tailed Lévy motions[END_REF]. For a general f , an explicit formula for the equilibrium distribution is not known, however when the noise level ε is small enough, finer characterizations of the structure of the equilibrium density in dimension one is known. We next summarize known results in this area, which show that Lévy-driven dynamics spend more time in 'wide valleys' in the sense of [START_REF] Chaudhari | Entropy-SGD: Biasing gradient descent into wide valleys[END_REF] when ε goes to zero.

Assume that f is smooth with r local minima {m i } r i=1 separated by r -1 local maxima {s i } r-1 i=1 , i.e. -∞ := s 0 < m 1 < s 1 < • • • < s r-1 < m r < s r := ∞.
Furthermore, assume that the local minima and maxima are not degenerate, i.e. f (m i ) > 0 and f (s i ) < 0 for every i. We also assume the objective gradient has a growth condition f (w) > |w| 1+c for some constant c > 0 and when |w| is large enough. Each local minima m i lies in the (interval) valley

S i = (s i-1 , s i ) of (width) length L i = |s i -s i-1 |.
Consider also a δ-neighborhood B i := {|x -m i | ≤ δ} around the local minimum with δ > 0 small enough so that the neighborhood is contained in the valley S i for every i.

We are interested in the first exit time from B i starting from a point w 0 ∈ B i and the transition time T i w0 (ε) := inf{t ≥ 0 : w ε t ∈ ∪ j =i B j } to a neighborhood of another local minimum, we will remove the dependency to w 0 of the transition time in our discussions as it is clear from the context. The following result shows that the transition times are asymptotically exponentially distributed in the limit of small noise and scales like 1/ε α with ε.

Theorem 2 [START_REF] Pavlyukevich | Cooling down lévy flights[END_REF]). For an initial point w 0 ∈ B i , in the limit ε → 0, the following statements hold regarding the transition time:

P w0 (w ε T i (ε) ∈ B j ) → q ij q -1 i if i = j, ε α T i (ε) d -→ E(q i ), E[ε α T i (ε)] → q -1 i .
where E denotes the exponential distribution, d -→ denotes convergence in distribution and

q ij = 1 α 1 |s j-1 -m i | α - 1 |s j -m i | α , i = j, (4.8)
q i = j =i q ij . (4.9)
If the SDE (4.7) would be driven by the Brownian motion instead, then an analogous theorem to Theorem 2 holds saying that the transition times are still exponentially distributed but the scaling ε α needs to be replaced by e 2H/ε 2 where H is the maximal depth of the basins to be traversed between the two local minima [START_REF] Day | On the exponential exit law in the small parameter exit problem[END_REF][START_REF] Bovier | Metastability in reversible diffusion processes ii: Precise asymptotics for small eigenvalues[END_REF]. This means that in the small noise limit, Brownian-motion driven gradient descent dynamics need exponential time to transit to another minimum whereas Lévy-driven gradient descent dynamics need only polynomial time. We also note from Theorem 2 that the mean transition time between valleys for Lévy SDE does not depend on the depth H of the valleys they reside in which is an advantage over Brownian motion driven SDE in the existence of deep valleys. Informally, this difference is due to the fact that Brownian motion driven SDE has to typically climb up a valley to exit it, whereas

Lévy-driven SDE could jump out.

The following theorem says that as ε → 0, up to a normalization in time, the process w ε t behaves like a finite state-space Markov process that has support over the set of local minima {m i } r i=1 admitting a stationary density π = (π i ) r i=1 with an infinitesimal generator Q. The process jumps between the valleys S i , spending time proportional to probability π i amount of time in each valley in the equilibrium where the probabilities π = (π i ) r i=1 are given by the solution to the linear system Qπ = 0.

Theorem 3 [START_REF] Pavlyukevich | Cooling down lévy flights[END_REF]). Let w 0 ∈ S i , for some 1 ≤ i ≤ r. For t ≥ 0, w ε tε -α → Y mi (t), as ε → 0, in the sense of finite-dimensional distributions, where Y = (Y y (t)) t≥0 is a continuous-time Markov chain on a state space {m 1 , m 2 , . . . , m r } with the infinitesimal generator Q = (q ij ) r i,j=1 with

q ij = 1 α 1 |s j-1 -m i | α - 1 |s j -m i | α ,
(4.10)

q ii = - j =i q ij . (4.11)
This process admits a density π satisfying Q T π = 0.

A consequence of this theorem is that equilibrium probabilities π i are typically larger for "wide valleys". To see this consider the special case illustrated in Figure 4.1 with r = 2 local minima m 1 < s 1 = 0 < m 2 separated by a local maximum at s 1 = 0. For this example, m 2 > |m 1 |, and the second local minimum lies in a wider valley. A simple computation reveals

π 1 = |m 1 | α |m 1 | α + m α 2 , π 2 = |m 2 | α |m 1 | α + |m 2 | α .
We see that π 2 > π 1 , that is in the equilibrium the process spends more time on the wider valley. In particular, the ratio π2 π1 = m2 |m1| α grows with an exponent α when the ratio m2 |m1| of the width of the valleys grows. Consequently, if the gradient noise is indeed α-stable distributed, these results directly provide theoretical evidence for the wide-minima behavior of SGD assuming the loss landscape is not degenerate.

In addition to the transition time between the basins of attraction of two local minima, understanding how long it takes for the continuous-time process w t given by (4.6) to exit a neighborhood of a local minimum w (given that it is started in that neighborhood) is also relevant. We formally define the first exit time of the stochastic process (4.6) as follows:

τ a (ε) inf{t ≥ 0 : |w t -w| ∈ [0, a]}.
(4.12)

The following result characterizes the first exit time in dimension one.

Theorem 4 [START_REF] Imkeller | First exit times of SDEs driven by stable Lévy processes[END_REF]). Consider the SDE (4.6) in dimension d = 1 and assume that it has a unique strong solution. Assume further that the objective f has a global minimum at zero, satisfying the conditions f (x)x ≥ 0 for every x ∈ R, f (0) = 0, f (x) = 0 if and only if x = 0, and f (0) > 0. Then, there exist positive constants ε 0 , γ, δ, and C > 0 such that for 0 < ε ≤ ε 0 , the following holds:

e -uε α θ α (1+Cε δ ) (1 -Cε δ ) ≤ P(τ a (ε) > u) ≤ e -uε α θ α (1-Cε δ ) (1 + Cε δ ) (4.13) for all initialization w 0 ∈ [-a + ε γ , a -ε γ ] and u ≥ 0, where θ = 2 a α . Consequently, E[τ a (ε)] = α 2 a α ε α (1 + O(ε δ )), for all w 0 ∈ [-a + ε γ , a -ε γ ]. (4.14)
Extension of Theorem 4 to R d . The exit behavior of the SDE (4.6) from an arbitrary domain in R d has also been studied in the literature. Imkeller et al. [2010a] generalizes Theorem 4 from dimension d = 1 to arbitrary dimensions and shows that in the small noise limit the exit time from a domain is exponentially distributed with a parameter that depends on the tail-index α. In case the components of the Lévy motion in (4.6) is replaced by a process that consists of the sum of finitely many one-dimensional Lévy processes with different tail-indices α i , it is also shown that the first exit time from a domain is determined by the smallest α i when the noise level ε is small enough.

Experimental methodology

In this section, we describe our experimental methodology regarding how we estimate the heavy-tailedness of the stochastic gradients. First, we discuss how we can compute the tail-index α based on a recent estimator proposed in [START_REF] Mohammadi | On estimating the tail index and the spectral measure of multivariate α-stable distributions[END_REF].

Second, we describe the procedure proposed in [START_REF] Brcich | The stability test for symmetric alphastable distributions[END_REF] for testing whether stochastic gradients follow a symmetric α-stable distribution. Our experimental results will be presented in Section 4.6.

Tail index estimation

Estimating the tail-index of an extreme-value distribution is a long-standing topic. Some of the well-known estimators for this task are [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF][START_REF] Pickands | Statistical inference using extreme order statistics[END_REF][START_REF] Dekkers | A moment estimator for the index of an extreme-value distribution[END_REF][START_REF] De Haan | Comparison of tail index estimators[END_REF]]. Despite their popularity, these methods are not specifically developed for α-stable distributions and it has been shown that they might fail for estimating the tail-index for α-stable distributions [START_REF] Mittnik | Tail estimation of the stable index α[END_REF]Rachev, 1996, Paulauskas and[START_REF] Paulauskas | Once more on comparison of tail index estimators[END_REF].

In this section, we use a relatively recent estimator proposed in [START_REF] Mohammadi | On estimating the tail index and the spectral measure of multivariate α-stable distributions[END_REF] for α-stable distributions. It is given in the following theorem.

Theorem 5 [START_REF] Mohammadi | On estimating the tail index and the spectral measure of multivariate α-stable distributions[END_REF]). Let {X i } K i=1 be a collection of random variables with X i ∼ SαS(σ) and

K = K 1 × K 2 . Define Y i K1 j=1 X j+(i-1)K1 for i ∈ 1, K 2 .
Then, the estimator

1 α 1 log K 1 1 K 2 K2 i=1 log |Y i | - 1 K K i=1 log |X i | . (4.15) converges to 1/α almost surely, as K 2 → ∞.
As shown in Theorem 2.3 of [START_REF] Mohammadi | On estimating the tail index and the spectral measure of multivariate α-stable distributions[END_REF], this estimator admits a faster convergence rate and smaller asymptotic variance than all the aforementioned methods.

In order to verify the accuracy of this estimator, we conduct a preliminary experiment, where we first generate K = K 1 × K 2 many SαS(1) distributed random variables with K 1 = 100, K 2 = 1000 for 100 different values of α. Then, we estimate α by using α ( 1 α ) -1 . We repeat this experiment 100 times for each α. As shown in Figure 4.3, the estimator is very accurate for a large range of α. Due to its favorable theoretical properties such as independence of the scale parameter σ, combined with its empirical stability, we choose this estimator in our experiments.

In order to estimate the tail-index α at iteration k, we first partition the set of data points D {1, . . . , n} into many disjoint sets Ω i k ⊂ D of size b, such that the union of these subsets gives all the data points. Formally, for all i, j = 1, . . . , n/b,

|Ω i k | = b, ∪ i Ω i k = D, and Ω i k ∩ Ω j k = ∅ for i = j.
This approach is similar to sampling without replacement. We then compute the full gradient ∇f (w k ) and the stochastic gradients

∇ fΩ i k (w k ) for each minibatch Ω i k .
We finally compute the stochastic gradient noises

U i k (w k ) = ∇ fΩ i k (w k ) -∇f (w k ), vectorize each U i k (w k
) and concatenate them to obtain a single vector, and compute the reciprocal of (4.15). In this case, we have K = dn/b and we set K 1 to the divisor of K that is the closest to √ K.

Stability test

Besides estimating the tail-index of a random process, it is also important to verify whether the process is symmetric α-stable. In this section, we describe a procedure [START_REF] Brcich | The stability test for symmetric alphastable distributions[END_REF]) for obtaining a confidence level for the stability of a random process, based on the following property:

Theorem 6 [START_REF] Brcich | The stability test for symmetric alphastable distributions[END_REF]). A necessary and sufficient condition for a random variable X to have an SαS distribution is

X 1 + X 2 ∼ C 1 X (4.16) X 1 + X 2 + X 3 ∼ C 2 X (4.17)
where C 1 , C 2 > 0 and X 1 , X 2 and X 3 are independent copies of X.

Here we adopt the stability test presented in [START_REF] Brcich | The stability test for symmetric alphastable distributions[END_REF]. To obtain a statistical test from (4.16), we first separate the observations into three equal-size subsets X, X 1 and X 2 , which are considered as independent copies of the observations. We then assign the first subset X to the right side of (4.16) and estimate the tail index α X of this subset using the idea of the previous section. For the left side of (4.16), we sum X 1 and X 2 term by term, and estimate α 12 of the resulting sum. Similarly, by separating the observations into four equal-size subsets X , X 1 , X 2 and X 3 , then repeating these above steps, we get α X from X and α 123 from X 1 + X 2 + X 3 , for a statistical test of (4.17). In the end, the process is considered to be α-stable if the tail indices estimated from the left and the right sides of (4.16) (as well as of (4.17

Numerical results

We investigate the tail behavior of the stochastic gradient noise in a variety of scenarios. We first consider a fully-connected network (FCN) on the MNIST and CIFAR10

datasets. For this model, we vary the depth (i.e. the total number of layers) in the set {2, 3, . . . , 10}, the width (i.e. the number of neurons per hidden layer) in the set {2, 4, 8, . . . , 1024}, and the minibatch size ranging from 1 to full batch.

We then consider a convolutional neural network (CNN) architecture (AlexNet) on the CIFAR10 and CIFAR100 datasets. We scale the number of filters in each convolutional layer in range {2, 4, . . . , 512}. We use the existing random split of the MNIST dataset into train and test parts of sizes 60K and 10K, and CIFAR10 and CIFAR100

datasets into train and test parts of sizes 50K and 10K, respectively. The order of the total number of parameters d range from several thousands to tens of millions.

For both FCN and CNN, we run each configuration with the negative-log-likelihood (i.e. cross entropy) and with the linear hinge loss, and we repeat each experiment with three different random seeds (see [START_REF] Geiger | The jamming transition as a paradigm to understand the loss landscape of deep neural networks[END_REF] for details on the choice of the hinge loss). The training algorithm is SGD with no explicit modification such as momentum or weight decay. The training runs for a fixed number of iterations unless it hits 100% training accuracy first. At every 100th iteration, we log the full training and test accuracies, and the tail estimate of the gradients that are sampled using the corresponding minibatch size. The codebase is implemented in python using pytorch2 .

Below, we present the most relevant and representative results. We have observed that, in all configurations, the three different initializations yielded no significant difference. Therefore, the effects of the randomness in initialization (under a given scheme)

do not appear to affect the gradient noise. Similarly, the choice of the loss function do not yield different behaviors in terms of the tail index. Even though the heavy tailed nature remains the same, the choice of the loss function results in a different way of dependence to the hyperparameters of the system, which we discuss in Section 4.6.4 and leave the investigation to a further study.

Stability test results

We first start by investigating the stability of the stochastic gradient noises under the datasets that we use for our estimation experiments. We will first focus on the later iterations of SGD, where the tail-index becomes stationary. Using an FCN on the MNIST and CIFAR10 datasets, we estimate the condition number c st (as described in Section 4.5.2) at every 50th iteration of the training stage, then take its average over the last 10K iterations to get the final result. Here, we consider c st ≈ 0.05 to be an acceptable level for the test since it is a quite small number with respect to estimated α in our experiments.

The results using the MNIST dataset are illustrated by As can be seen from the figure, the condition c st of the network fails to be smaller than our required criterion in some cases. However, the gap from this number to the criterion is quite small that we can consider that it does not violate the α-stable assumption on the gradient noise of the network. Unlike the MNIST dataset, we observe that for the networks with 256 neurons per layer, even though the overall gradient noise strongly exhibits an α-stable behavior, some of the hidden layers are very far from being α-stable, suggesting that the characteristics of the first layer dominate the overall structure. In contrast, the gradient noise with respect to the parameters of the hidden layers becomes more α-stable with a very high number (512) of neurons per layer.

By these experiments, we observe that the structure of the dataset has a strong impact on the statistical properties of the gradient noise, especially for the layers with smaller number of parameters. When this number of parameters is large (which is usually the case in practice), the gradient noise corresponding to these parameters becomes more α-stable. In short, this means increasing the size of the network (the number of the network parameters) tends to make the gradient noise behave similarly to an α-stable noise.

4.6.2 Effect of varying network size We measure the tail-index for varying the widths and depths for the FCN, and varying widths (i.e. the number of filters) for the CNN. For very small sizes, the networks perform poorly; therefore, we only illustrate sufficiently large network sizes, which yield similar accuracies. For these experiments, we compute the average of the tail-index measurements for the last 10K iterations (i.e. when α becomes stationary) to focus on the late stage dynamics.

Figure 4.5 shows the results for the FCN. The first striking observation is that in all the cases, the estimated tail-index is far from 2, meaning that the distribution of the gradient noise is highly non-Gaussian. For the MNIST dataset, we observe that α systematically decreases for increasing network size, where this behavior becomes more prominent with the depth. This result shows that, for MNIST, increasing the dimension of the network results in a gradient noise with heavier tails and therefore increases the probability of ending up in a wider basin. For the CIFAR10 dataset, we still observe that α is far from 2; however, in this case, increasing the network size does not have a clear effect on α. In all cases, we observe that α is in the range 1.1-1.2. From left to right: depth 3, depth 5, depth 7. Different lines correspond to different widths.

In Figure 4.6, we plot estimated α for each layer of FCNs, using MNIST dataset where the minibatch is of size 100. The resulting α is obtained by averaging α over the last 10K iterations. The layer index '0' corresponds to the estimated α of the whole network. In this experiment, we observe that α becomes smaller (heavier-tailed) for the deeper layers. In addition, the value of the tail-index for the whole network has a strong connection with the first layers: the α for the whole network is closer to that of the first layers than of the last layers. show that, for both CIFAR10 and CIFAR100, the tail-index is extremely low for the under-parametrized regime (e.g. the case when the width is 2, 4, or 8 for CIFAR10). As we increase the size of the network the value of α increases until the network performs reasonably well and stabilizes in the range 1.0-1.1. We also observe that α behaves similarly for both train and test sets3 .

These results show that there is strong interplay between the network architecture, dataset, and the algorithm dynamics: (i) we see that the size of the network can strongly influence α, (ii) for the exact same network architecture, the choice of the dataset has a significant impact on not only the landscape of the problem, but also the noise characteristics, hence on the algorithm dynamics. The results clearly show that there are two distinct phases of SGD (in this configuration before and after iteration 1000). In the first phase, the loss decreases very slowly, the accuracy slightly increases, and more interestingly α rapidly decreases. When α reaches its lowest level, the process possesses a jump, which causes a sudden decrease in the accuracy. After this point the process recovers again and we see a stationary behavior in α and an increasing behavior in the accuracy. We also investigate this behavior for each layer of an FCN with depth 7 and width 512 in Figure 4.9. The estimated tail-index for each layer has a clear phase change at earlier iterations, where we observe that this jump is more prominent in the deeper layers where the tail-index is smaller. On the other hand, unlike the whole network, the tail-index of each layer undergoes a fluctuation period before becoming stationary at the last 2000 iterations. However, this observation might be due to the measurement error since the size of the sample that is used in the estimator (4.15) gets smaller when we make layer-wise measurements.

The fact that the process has a jump when α is at its smallest value provides a strong support to our assumptions and the metastability theory that we discussed in the previous section. Furthermore, these results also strengthen the view that SGD crosses barriers at the very initial phase and continues searching until it reaches a "wide and flat enough" region of a local optimum. On the other hand, our current analysis is not able to determine whether the process jumps in a different basin or a 'better' part of the same basin and we leave it as a future work.

A note on generalization

In this section, we investigate the connection between the tail-index and the generalization performance. In particular, we consider the relation between α and the ratio of the step-size to the batch size η/b which is proportional to the noise scale of SGD when there is no momentum [START_REF] Park | The effect of network width on stochastic gradient descent and generalization: an empirical study[END_REF]. It has been empirically demonstrated that this ratio correlates with performance of the model [START_REF] Jastrzebski | Three factors influencing minima in SGD[END_REF], hence the higher the noise scale, the better the generalization performance until a certain level.

Clearly, when the noise is too high, training may diverge, however, proper level of noise leads to better solutions.

In this section, we will investigate how the tail index of the gradient noise is affected for different noise scales. We reproduce and follow the initialization convention and the hyper-parameter scale that is studied in [Park et al., 2019, Appendix G]: A fully connected model with 3 hidden layers, each hidden layer has 512 nodes. Weights are initialized ∼ N (0, 1), bias terms are set to zero at the initial point. Each layer is then passed through ReLU non-linearity, and multiplied by the inverse of the width of the previous layer. As usual, the network is trained with SGD without momentum; the dataset is the standard MNIST. Minibatch size ranges in the set [24,48,96,192] In both choices of loss functions, hinge and NLL, the behavior of the test error with respect to the noise scale is consistent with previous observations. Similarly, in both cases, the estimated α remains within a narrow band of 1, indicating the heavy tail behavior. However, the trends in estimated α are different depending on the choice of the loss. Therefore, we cannot attribute the improvement in performance to lower α when increasing the noise scale. To better emphasize this point, we plot the correlation of estimated α and test error in Figure 4.11 where the positive and negative correlations are clearly visible depending on the choice of the loss function. This contrasting behavior is another hint that there exists a connection between α and the test performance (since they are correlated in both cases) and suggests us to examine this connection in order to understand when exactly the dynamics falls into basins with better performance.

Chapter 5

Global convergence analysis

In this chapter, we investigate the global convergence property of stochastic gradient descent (SGD) for non-convex optimization by using a stochastic process driven by αstable distribution. Our results show that the weak-error by SGD analyzed by this method increases faster than when analyzed using a stochastic process driven by Gaussian distribution, which suggests using smaller step-sizes.

This chapter is based on the article [Nguyen et al., 2019b].

Summary of the main result

So far, we have illustrated that stochastic processes driven by α-stable noise can be used as a proxy for understanding the dynamics of stochastic gradient descent in deep learning. Therefore, to understand the behaviors of SGD, we will provide an analysis of the non-asymptotic behavior of the following algorithm for non-convex optimization:

w k+1 = w k -ηc α ∇f (w k ) + η β 1 α S k+1 , (5.1)
where α ∈ (1, 2], c α is a known constant, β ∈ R + is called the 'inverse temperature' parameter and {S k } k∈N+ is a sequence of α-stable distributed random variables.

Even though asymptotic convergence properties of (5.1) were established for decreasing step-sizes in S ¸imşekli [2017], [START_REF] Panloup | Recursive computation of the invariant measure of a stochastic differential equation driven by a Lévy process[END_REF], these results do not explain the behavior of SGD algorithm for finite number of iterations. Besides, in practice, using a constant step-size often yields better performance [START_REF] Baker | sgmcmc: An R package for stochastic gradient Markov chain Monte Carlo[END_REF], a situation which cannot be handled by the existing theory.

In particular, we analyze the expected suboptimality E[f (w k ) -f ], where f f (w ) = min w∈R d f (w). As we will describe in detail in Section 5.3, we decompose this suboptimality into four different terms, and we bound each of those terms one by one.

Due to the choice of the α-stable Lévy motion, the standard tools for analyzing SDEs driven by a Brownian motion are not available for our use, and therefore, we cannot use the proof strategies developed for ULA (3.2) as they are (such as [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF], [START_REF] Xu | Global convergence of langevin dynamics based algorithms for nonconvex optimization[END_REF], [START_REF] Erdogdu | Global non-convex optimization with discretized diffusions[END_REF]). Instead, we follow an alternative path, where we first relate the expected discrepancies to Wasserstein distance of fractional orders, and then, inspired by [START_REF] Gairing | Transportation distances and noise sensitivity of multiplicative Lévy SDE with applications[END_REF], we prove a result that expresses the Wasserstein distance (Definition 1) between the laws of two SDEs (driven by α-stable Lévy motion) in terms of their drift functions.

Informally, we show that the expected suboptimality E[f (w k ) -f ] is bounded by a sum of four terms, summarized as follows:

E[f (w k ) -f ] ≤ A 1 + A 2 + A 3 + A 4 ,
where

A 1 = O k 1+max{ 1 q ,γ+ γ q } η 1 q , A 2 = O k 1+max{ 1 q ,γ+ γ q } η 1 q + γ αq d β (q-1)γ αq , A 3 = O β + d exp - λ * kη β , A 4 = O 1 β γ+1 + d β log(β + 1) .
Here γ ∈ (0, 1) is the Hölder exponent of the gradients of f , and q ∈ (1, α), λ * > 0 are some constants. This result has the following implications. For any ε > 0,

1. If 1 q > γ + γ q and k ε -1 and η < ε 2q+1 , then A 1 scales as Cε and A 2 scales as εPoly(β, d).

If 1

q ≤ γ + γ q and k ε -1 and η < ε 2q+γ+γq , then A 1 scales as Cε and A 2 scales as εPoly(β, d).

3. If we choose kη > β λ * log 1 ε , then A 3 scales as εPoly(β, d). where Poly(. . .) denotes a formal polynomial, i.e., an expression containing the realordered exponents of the variables, coefficients, and only the operations of addition, subtraction, and multiplication.

In Section 5.5, we extend our results in two directions: (i) obtaining guarantees for Bayesian posterior sampling and (ii) non-convex optimization where exact gradients are replaced with stochastic gradients. Our results imply that, in the context of global optimization, the error induced by (5.1) has a worse dependency on k and η, as compared to ULA. This suggests that one should use smaller step-sizes in (5.1).

Assumptions and the main result

We start by defining three different stochastic processes x 1 (t), x 2 (t), and x 3 (t), which will be the main constructs in our analysis. We first informally define these processes as follows: x 2 is a continuous-time process that interpolates w k in time and it will let us avoid dealing with the discrete-time process w k directly. x 1 is the limiting process of w k when the step-size goes to zero. Finally, x 3 is a process whose law converges to the Gibbs measure π, whose density is exp(-βf (w)) (up to a multiplicative factor).

In our approach, we will first relate x 2 to its limiting process x 1 . Since it is more challenging to relate x 1 to w , we will then relate x 1 to x 3 , and x 3 to π. By following a similar approach to [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF], we will finally relate π to f . Formally, we decompose the expected suboptimality in the following manner:

Ef (w k ) -f * = Ef (x 2 (kη)) -Ef (x 1 (kη)) + Ef (x 1 (kη)) -Ef (x 3 (kη)) + Ef (x 3 (kη)) -Ef ( ŵ) + Ef ( ŵ) -f * , (5.2)
where x i (kη) with i = 1, 2, 3 denotes the state reached by the three stochastic processes at time kη, and ŵ is a random variable drawn from π. We will now formally define the processes x 1 , x 2 , and x 3 .

The first SDE is the continuous-time limit of the algorithm given in (5.1) and defined as follows for t ≥ 0:

dx 1 (t) = b 1 (x 1 (t-), α)dt + β -1/α dL α t , (5.3)
where the drift function has the following form:

b 1 (x, α) -c α ∇f (x).
The second SDE is a linearly interpolated version of the discrete-time process {w k } k∈N+ , defined as follows:

dx 2 (t) = b 2 (x 2 , α)dt + β -1/α dL α t , (5.4)
where x 2 ≡ {x 2 (t)} t≥0 denotes the whole process and the drift function is chosen as follows:

b 2 (x 2 , α) -c α ∞ k=0 ∇f (x 2 (kη))I [kη,(k+1)η[ (t).
Here, I denotes the indicator function, i.e. for any set [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF].

A ⊂ R d , I A (x) = 1 if x ∈ A and I A (x) = 0 if x / ∈ A. It is easy to verify that x 2 (kη) = w k for all k ∈ N + Dalalyan [2017b],
The last SDE is designed in such a way that its solution has the Gibbs distribution as the invariant distribution and is defined as follows:

dx 3 (t) = b(x 3 (t-), α)dt + β -1/α dL α t , (5.5)
where the drift is a d-dimensional vector whose i-th component, i = 1, . . . , d, has the following form:

(b(x, α)) i - D α-2 i φ(x)∂ i f (x) φ(x) .
(5.6)

Here, D i denotes the Riesz derivative (defined in Theorem 1) along the direction i Ortigueira et al. [2014] and ∂ i denotes the derivative with respect to the i-th component.

With this definition for the drift, we have the following result for the invariant measure of x 3 , which is an extension of Theorem 1 to general d and β.

Lemma 1. The SDE (5.5) with drift b defined by (5.6) admits π as an invariant distribution of its solution (x 3 (t)) t≥0 .

The process {x 3 (t)} t will play an important role in our analysis, since it will enable us to relate x k to the Gibbs measure π, whose samples will be close to the global optimum w with high probability [START_REF] Pavlyukevich | Cooling down lévy flights[END_REF].

We now state our assumptions that will be used to imply the main result of the chapter.

A 1. There exists a constant B ≥ 0 such that

c α ∇f (0) ≤ B.
A 2. The gradient of f is Hölder continuous with constants M > 0, 0 ≤ γ < 1:

c α ∇f (x) -∇f (y) ≤ M x -y γ , ∀x, y ∈ R d . A 3. For some m > 0 and b ≥ 0, f is (m, b, γ)-dissipative: c α x, ∇f (x) ≥ m x 1+γ -b, ∀x ∈ R d .
The assumptions A1-3 are mild and when γ = 1, they become the standard Lipschitz and dissipativity conditions that are often considered in diffusion-based non-convex optimization algorithms [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF], [START_REF] Xu | Global convergence of langevin dynamics based algorithms for nonconvex optimization[END_REF], [START_REF] Erdogdu | Global non-convex optimization with discretized diffusions[END_REF].

However, due to the choice of the α-stable Lévy motion with α ∈ (1, 2), we need to consider a 'fractional' version of those assumptions and exclude the case where γ = 1, which makes A3 weaker and A2 more restrictive than the case where γ = 1. Nevertheless, A2 can be weakened to local Hölder continuity by using the localization techniques given in the proof of Proposition 4.2.2 of [START_REF] Kunze | Stochastic differential equations[END_REF]. This approach requires rewriting all the expressions which employ A2 in our proofs, by using stopping-times in such a way that we can adopt the local Hölder continuity in the same manner of [START_REF] Kunze | Stochastic differential equations[END_REF].

In our analysis, we will make a repeated use of the Hölder and Minkowski inequalities, which require the following condition to hold:

A 4. There exist positive real numbers p, q, p 1 , q 1 such that

1 p + 1 q = 1 p 1 + 1 q 1
= 1, and q < α, γp < 1, γq 1 < 1, (q -1)p 1 < 1.

Even though this assumption looks rather technical, when combined with A2-3, it will in fact impose smoothness constraints on f and restrict γ to be less than 1. We will discuss this observation in more detail in Section 5.4.

Next, we require the drift b (presented in equation (5.5)) to be dissipative for large distances and we assume a bounded moment condition, which will be used for establishing the ergodicity of x 3 .

A 5. 1) For all x, y ∈ R d and for some constants γ ∈ [0, 1], l 0 ≥ 0, K 1 > 0 and K 2 > 0, the following holds:

b(x) -b(y), x -y x -y ≤    K 1 x -y γ , x -y < l 0 , -K 2 x -y , x -y ≥ l 0 .
2) For any t > 0, γ ∈ (0, α), and for any coupling P t of x 3 (t) and ŵ ∼ π, we have:

x 3 (t) -ŵ γ dP t < C * ,
for some constant C * > 0.

Remark The first part of assumption A5 can be satisfied under a set of (rather nontrivial) assumptions (see appendix for more details).

Proposition 1. Under assumptions A1-3 and A5, the distribution of x 3 (t) exponentially converges to its unique invariant distribution π in the Wasserstein metric, i.e., for any λ ≥ 0 such that λ < α, there exist constants C > 0 and C 1 > 0 such that

W λ (µ 3t , π) ≤ Ce -C1t , (5.7)
where µ 3t denotes the probability measure of x 3 (t).

In the rest of the paper, we will assume that the constants C and C 1 behave similarly to the case of the unadjusted Langevin algorithm (α = 2). In particular, we assume that C is proportional to β and C 1 is proportional to β -1 , so that we can rewrite (5.7) as follows:

W λ (µ 3t , π) ≤ Cβe -λ * t/β .
In the unadjusted Langevin algorithm, the constant λ * turns out to be the uniform spectral gap associated with the Gibbs measure π and it has shown to scale exponentially with respect to the dimension d in the worst case [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF]. We believe that a similar property holds in our case as well.

Our next assumption is on the approximation quality of the function b by b 1 .

A 6. There exists a constant L > 0 such that L < m and sup

x∈R d c α ∇f (x) + b(x, α) ≤ L,
where the function b is defined in (5.6).

In Corollary 2 of S ¸imşekli [2017], it has been shown that A6 holds if the tails of π vanish sufficiently quickly (cf. assumption H4 in S ¸imşekli [2017]). On the other hand, the gap between functions b and b 1 can be diminished even more if we consider a more sophisticated numerical approximation scheme, such as the one given in C ¸elik and Duman

[2012] (cf. Theorem 2 of S ¸imşekli [2017]).

In our final condition, we assume that the fractional moments of π is uniformly bounded.

A 7. There exists a constant C > 0 such that

R d x r π(dx) ≤ C b + d/β m for all 0 ≤ r ≤ 2.
Now, we are ready to state the main result of this chapter.

Theorem 7. Under conditions A1-7 and for 0 < η < m M 2 , there exists a positive constant C independent of k and η such that the following bound holds:

E[f (w k )] -f * ≤C k 1+max{ 1 q ,γ+ γ q } η 1 q + k 1+max{ 1 q ,γ+ γ q } η 1 q + γ αq d β (q-1)γ αq + βb + d m exp(- λ * kη β ) + M c -1 α β γ+1 (1 + γ) + 1 β log (2e(b + d β )) d 2 Γ( d 2 + 1)β d (dm) d 2 .
More details on constant C can be found in appendix.

Theorem 7 provides finite-time guarantees for discrete-time dynamics (4.5) in terms of suboptimality with respect to the global minimum as a function of the step-size and the scale parameter. In particular, it is shown that the heavy-tailed system has a worse dependency on both K and η as compared to the Gaussian case. Besides, it is known that if the scale parameter σ gets smaller, the dynamics admits a stationary distribution that will concentrate more and more on the global minimizer although reaching out to stationary would require an exponential number of steps in the dimension in the worst case. Similar to ULA [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF], our bound grows with the number of iterations k. We note that this result sheds light on the explicit dependency of the error with respect to the algorithm parameters (e.g. step-size) for a fixed number of iterations, rather than explaining the asymptotic behavior when k goes to infinity. In the next sections, we will provide an overview of the proof of this theorem along with some remarks and comparisons to ULA.

Proof overview

Our proof strategy consists of bounding each of the four terms in (5.2) separately. Before bounding these terms, we first start by relating the expected discrepancies to the Wasserstein distance between two random processes. The result is formally presented in the following lemma and it extends the 2-Wasserstein continuity result given in [START_REF] Polyanskiy | Wasserstein continuity of entropy and outer bounds for interference channels[END_REF] to Wasserstein distance with fractional orders.

Lemma 2. Let v and w be two random variables on R d which have µ and ν as the probability measures and let g be a function in

C 1 (R d , R). Assume that for some c 1 > 0, c 2 ≥ 0 and 0 ≤ γ < 1, ∇g(x) ≤ c 1 x γ + c 2 , ∀x ∈ R d and max E w γp 1 p , E v γp 1 p < ∞.
Then, the following bound holds:

gdµ -gdν ≤C W q (µ, ν),
for some C > 0.

Lemma 2 lets us upperbound the first three terms of the right hand side of (5.2) by the Wasserstein distance between the appropriate stochastic processes, respectively W q (µ 1t , µ 2t ), W q (µ 1t , µ 3t ), and W q (µ 3t , π), where µ it denotes the law of x i (t).

The term W q (µ 3t , π) is related to the ergodicity of the process (5.5) and it has been shown that this distance diminishes exponentially for a considerably large class of Lévy diffusions [START_REF] Masuda | Ergodicity and exponential β-mixing bounds for multidimensional diffusions with jumps[END_REF], Xie and Zhang [2017]. On the other hand, the term W q (µ 1t , µ 3t ) is related to the numerical approximation of the Riesz derivatives, which is analyzed in S ¸imşekli [2017]. Therefore, in this study, we use the assumptions A5 and A6

for dealing with these terms, and focus on the term W q (µ 1t , µ 2t ), which is related to the so-called 'weak-error' of the Euler scheme for the SDE (5.3). The existing estimates for such weak-errors are typically of order Cη a , where a < 1 and C is a constant that grows exponentially with t Mikulevičius and Zhang [2011]. The exponential growth with t is prohibitive in our case and one of our main technical contributions is that, in the sequel, we will prove a bound that grows polynomially with t, which substantially improves over the one with exponential growth.

We start by bounding W q (µ 1t , µ 2t ) and W q (µ 1t , µ 3t ). In order to do so, we prove the following lemma, which will be the key for our analysis.

Lemma 3. For λ ∈ (1, ∞), i, j ∈ {1, 2, 3} and i = j, we have the following identity:

W λ (µ it , µ jt ) = inf E t 0 λ ∆x ij (s) λ-2 ∆x ij (s), ∆b ij (s-) ds 1/λ
, where the infimum is taken over the couplings whose marginals are µ it and µ jt and

∆x ij (s) x i (s) -x j (s), ∆b ij (s-) b i (x i (s-), α) -b j (x j (s-), α).
This result extends the recent study [START_REF] Gairing | Transportation distances and noise sensitivity of multiplicative Lévy SDE with applications[END_REF] and lets us relate the Wasserstein distance between the distributions of the random processes to their drift functions.

By using Lemma 3, we start by bounding the Wasserstein distance between µ 1t and µ 2t . The result is summarized in the following theorem.

Theorem 8. Assume that the following condition holds: 0 < η ≤ m M 2 . Then, we have

W q q (µ 1t , µ 2t ) ≤ Cq Poly(k, η, β, d),
for some C > 0.

The full statement of the proof and the explicit constants are provided in appendix.

By only considering the leading terms of the bound provided in Theorem 8, we obtain the following corollary.

Corollary 1. Suppose that 0 < η < min 1, m M 2 . Then, the bound for the Wasserstein distance between the laws of x 1 (t) and x 2 (t) can be written as follows:

W q q (µ 1t , µ 2t ) ≤C(k 2 η + k 2 η 1+γ/α β -(q-1)γ/α d).
By combining Corollary 1 with Lemma 2, we obtain the following result, which provides an upperbound for the first term of the right hand side of (5.2).

Corollary 2. For 0 < η < m M 2 , there exists a constant C > 0 such that the following bound holds:

E[f (x 1 (kη))] -E[f (x 2 (kη))] ≤ C k 1+ 1 q η 1 q + k 1+ 1 q η 1 q + γ αq β -(q-1)γ αq d .
Remark For any ε > 0, if we choose k ε -1 Poly(β, d) and η < ε 2q+1 Poly(β, d), then the bound in Corollary 2 scales as εPoly(β, d).

Next, by using a similar approach, we bound the distance between µ 1t and µ 3t . In the next theorem, we show that the error grows polynomially with the parameters.

Theorem 9. We have the following estimate:

W q q (µ 1t , µ 3t ) ≤CqPoly(k, η, β, d)
By considering the leading terms of the bound in Theorem 9 and combining it with Lemma 2, we obtain the following corollaries.

Corollary 3. There exists a constant C ≥ 0 such that the following bound holds:

W q q (µ 1t , µ 3t ) ≤C(k q+γ η + k q+γ η q β -q-1 α d).
Corollary 4. There exists a constant C ≥ 0 such that the following inequality holds:

|E[f (x 1 (kη))] -E[f (x 3 (kη))]| ≤ C k γ+ γ+q q η γ+ 1 q β -γ α d + k γ+ γ+q q η 1 q . Remark For any ε > 0, if we choose k ε -1 Poly(β, d) and η < ε 2q+γq+γ Poly(β, d),
then the bound in Corollary 4 scales as εPoly(β, d).

We now pass to the term Ef (x 3 (kη)) -Ef ( ŵ) of (5.2). Since we have that µ 3t exponentially converges to π in Wasserstein distance (Proposition 1), as a direct application of Lemma 2, we obtain the following result.

Lemma 4. Let ŵ be a random variable drawn from the invariant measure π ∝ exp(-βf ) of (5.5). There exists a constant C ≥ 0 such that the following bound holds:

|E[f (x 3 (t))] -E[f ( ŵ)]| ≤ C bβ + d m exp(-λ * β -1 t).
Remark For any ε > 0, if we take kη > β λ * log 1 ε , then the bound in Lemma 4 can be scaled as εPoly(β, d).

We finally bound the term Ef ( ŵ) -f * , which is the expected suboptimality of a sample from π. By following a similar proof technique presented in Raginsky et al.

[2017], we obtain the following result.

Lemma 5. For β > 0, we have

E[f ( ŵ)] -f ≤β -1 log (2e(b + d β )) d/2 Γ( d 2 + 1)β d (dm) d/2 + β -γ-1 M c -1 α 1 + γ .
Combining Corollary 2, Corollary 4, Lemma 4, and Lemma 5 proves Theorem 7.

Additional remarks

In this part, we first compare our global convergence result with those for Unadjusted Langevin Algorithm presented in [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF]. Then, we discuss the feasibility of the smoothness assumption A4 and provide an explicit condition such that assumption A4 holds.

Comparison with ULA

Let us compare this result with those for ULA presented in [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF], since they use a similar decomposition (as opposed to [START_REF] Xu | Global convergence of langevin dynamics based algorithms for nonconvex optimization[END_REF]). The last two terms of the right hand side of the bound in Theorem 7 have less importance as they can be made arbitrarily small by increasing β. Besides, for β large enough, the first two terms in our bound can be combined in a single term that scales in the order of k 1+max{ 1 q ,γ+ γ q } η 1 q . The corresponding term for ULA is given as follows: kη 5/4 , cf. Section 3.1 of [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF]. This observation shows that (5.1) has a worse dependency both on k and η, which is not surprising and indeed in-line with the existing literature Mikulevičius and Zhang [2011].

Discussion on smoothness assumptions

In this section we will discuss assumption A4 and provide a condition on γ such that A4 holds. Let us recall the four constraints given in A4:

(1/p + 1/q) = (1/p 1 + 1/q 1 ) = 1 γp < 1, γq 1 < 1, (q -1)p 1 < 1.

We will refer to these conditions as the first, second, third, and fourth conditions, respectively. Our aim is to find a condition on γ (more precisely, the maximum value of γ) such that there exist p, q, p 1 , q 1 > 0 satisfying these four conditions. First, suppose that p > q 1 . Then, the maximum value of γ is decided by the second constraint. Since we want γ to be as large as possible, it is natural to choose a smaller p. We observe that, as we decrease p, due to the first and the fourth constraints, the value of q 1 needs to be increased. If we continue decreasing p, then q 1 continues to be increased and soon becomes strictly greater than p. At this moment, the maximum value of γ is decided by the third constraint, not by the second constraint anymore, and from this point on, it is more plausible to decrease q 1 . By this intuition, it is reasonable to choose p to be equal to q 1 , which implies that p 1 = q. Accordingly, the fourth constraint becomes: (q -1)q < 1. By noting that q > 1, solving this constraint gives 1 < q < (1 + √ 5)/2. Then by the first constraint, we have p > (3 + √ 5)/2, and the second constraint gives γ < 1/p < (3 -√ 5)/2. This upper bound for γ is a number between 0.38 and 0.39 and tells us that there exist p, q, p 1 , q 1 satisfying the four constraints if and only if 0 ≤ γ < (3 -√ 5)/2.

Let us take a closer look at Theorem 7. Since γ(q + 1)

< (3 - √ 5)(3 + √ 5)/4 = 1,
we have γ + γ/q = γ(q + 1)/q < 1/q. Hence, 1 + max{1/q, γ + γ/q} = 1 + 1/q.

Let ε 1 and ε 2 be positive numbers such that

1/q -ε 1 = 2/(1 + √ 5) = ( √ 5 -1)/2, γ + ε 2 = (3 - √ 5)/2.
then, if q = p 1 is approximately equal to (1 + √ 5)/2 and γ is approximately equal to (3 -√ 5)/2, we imply that ε 1 and ε 2 become very small and

1/q ≈ ( √ 5 -1)/2, 1/q + γ/(αq) ≈ ( √ 5 -1)/2 + ( √ 5 -2)/α, (q -1)γ/(qα) ≈ (7 -3 √ 5)/(2α).
For example, the values α = 1.65, γ = 0.38, p = q 1 = 2.63, q = p 1 = q/(q -1) ≈ 1.613

satisfy assumption A4. Hence, the bound in Theorem 7 can be expressed as follows:

Corollary 5. Under conditions A1-7, for α = 1.65, γ = 0.38, p = q 1 = 2.63, q = p 1 = q/(q -1) ≈ 1.613 and for 0 < η < m M 2 , there exists a positive constant C independent of k and η such that the following bound holds:

E[f (w k )] -f * ≤C k 1.62 η 0.61 + k 1.62 η 0.75 d β 0.0875 + βb + d m exp(- λ * kη β ) + M c -1 α 1.38β 1.38 + 1 β log (2e(b + d β )) d 2 Γ( d 2 + 1)β d (dm) d 2 .
Proof. The result is a direct consequence of Theorem 7.

As a final remark on this smoothness condition, we note that similar constraints are imposed on Lévy-driven SDEs in other studies as well [START_REF] Panloup | Recursive computation of the invariant measure of a stochastic differential equation driven by a Lévy process[END_REF], S ¸imşekli [2017]. This is due to the fact that such SDEs often require better-behaved drifts in order to be able to compensate the jumps incurred by the Lévy motion.

Extensions

In this section, we discuss the implications of our results in the classical Monte Carlo sampling context. Then, we provide a similar global convergence result for the case where the gradient is replaced by a stochastic gradient.

Guarantees for posterior sampling

If our aim is only to draw samples from the distribution π, then, for a fixed k, we can bound the Wasserstein distance between the law of w k and π. The result is stated as follows:

Corollary 6. For 0 < η ≤ m M 2 , the following bound holds:

W q (µ 2t , π) ≤C k max{2,q+γ} q η 1 q + k max{2,q+γ} q η 1 q + γ qα β -γ(q-1) qα d 1 q + βe -λ * kη β .
As a typical use case, we can consider Bayesian posterior sampling, where we choose β = 1 and

f (X) = -(log P(Y |X) + log P(X)).
Here, Y denotes a dataset, P(Y |X) is the likelihood, P(X) denotes the prior density, and the target distribution π becomes the posterior distribution with density P(X|Y ).

Extension to stochastic gradients

In many machine learning problems, the function f to be minimized has the following form:

f (x) 1 n n i=1 f (i) (x),
where i denotes different data points and n is the total number of data points. In large-scale applications, n can be very large, which renders the gradient computation infeasible. Therefore, at iteration k, we often approximate ∇f by its stochastic version that is defined as follows:

∇f k (x) 1 n s i∈Ω k ∇f (i) (x),
where Ω k is a random subset of {1, . . . , n} with |Ω k | = n s n. The quantity ∇f k (x) is often referred to as the 'stochastic gradient'. If the stochastic gradients satisfy a moment condition, then we have the following results:

Theorem 10. Assume that for each i, the function x → f (i) (x) satisfies the conditions A1-7. Let us replace ∇f by ∇f k in (5.1). If, in addition, there exists δ ∈ [0, 1) for any k, such that

E Ω k c α (∇f (x) -∇f k (x)) q1 ≤δ q1 M q1 x γq1 ,
for x ∈ R d , then we have the following bound:

W q q (µ 1t , µ 2t ) ≤C(1 + δ)(k 2 η + k 2 η 1+γ/α β -γ(q-1)/α d).
60 Similar to our previous bounds, we can use Theorem 10 for obtaining a bound for the expected discrepancy, given as follows:

Corollary 7. Under the same assumptions as in Theorem 10, we have the following bound:

E[f (x 1 (kη))] -E[f (x 2 (kη))] ≤ C(1 + δ) k 1+ 1 q η 1 q + k 1+ 1 q η 1 q + γ αq β -(q-1)γ αq d .
These results show that the guarantees for (5.1) still hold even under the presence of stochastic gradients.

61

Chapter 6

First exit time analysis

While approximating SGD as a continuous-time SDE brings a new perspective for analyzing SGD, this approach might not be accurate for any step-size η, and some theoretical concerns have already been raised for the validity of such approximations [START_REF] Yaida | Fluctuation-dissipation relations for stochastic gradient descent[END_REF].

Intuitively, one can expect that the metastable behavior of SGD would be similar to the behavior of its continuous-time limit only when the discretization step-size is small enough. Even though some theoretical results have been recently established for the discretizations of SDEs driven by Brownian motion [START_REF] Tzen | Local optimality and generalization guarantees for the Langevin algorithm via empirical metastability[END_REF], it is not clear that how the discretized Lévy SDEs behave in terms of metastability. In this chapter, we provide a formal theoretical analysis where we derive explicit conditions for the stepsize such that the metastability behavior of the discrete-time system is guaranteed to be close to its continuous-time limit.

This chapter is based on the article [Nguyen et al., 2019a].

First exit times of continuous-time Lévy stable

SDEs

Due to the discontinuities of the Lévy-driven SDEs, their metastability behaviors also differ significantly from their Brownian counterparts. In this section, we will briefly remind some important theoretical results, that we already mentioned in Section 4.4, about the following SDE:

dw t = -∇f (w t )dt + εdL α t (6.1)
For simplicity, let us consider this SDE in dimension one, i.e. d = 1. In a relatively recent study [START_REF] Imkeller | First exit times of SDEs driven by stable Lévy processes[END_REF], the authors considered this SDE, where the potential function f is required to have a non-degenerate global minimum at the origin. Their result (Theorem 4) indicates that the first exit time of w t needs only polynomial time with respect to the width of the basin and it does not depend on the depth of the basin, whereas Brownian systems need exponential time in the height of the basin in order to exit from the basin [START_REF] Bovier | Metastability in reversible diffusion processes i: Sharp asymptotics for capacities and exit times[END_REF], Imkeller et al. [2010b].

This difference is mainly due to the discontinuities of the Lévy motion, which enables it to 'jump out' of the basin, whereas the Brownian SDEs need to 'climb' the basin due to their continuity. Consequently, given that the gradient noise exhibits similar heavytailed behavior to an α-stable distributed random variable, this result can be considered as a proxy to understand the wide-minima behavior of SGD.

Let us remind that (see Section 4.4) this result has already been extended to R d in Imkeller et al. [2010a]. Extension to state dependent noise has also been obtained in [START_REF] Pavlyukevich | First exit times of solutions of stochastic differential equations driven by multiplicative Lévy noise with heavy tails[END_REF]. We also note that the metastability phenomenon is closely related to the spectral gap of the forward operator corresponding to the SDE dynamics (see e.g. [START_REF] Bovier | Metastability in reversible diffusion processes i: Sharp asymptotics for capacities and exit times[END_REF]) and it is known that this quantity scales like O(ε α ) for ε small which determines the dependency to ε in the first term of the exit time (4.14) due to Kramer's

Law [START_REF] Berglund | Kramers' law: Validity, derivations and generalisations[END_REF], [START_REF] Burghoff | Spectral analysis for a discrete metastable system driven by Lévy flights[END_REF]. [START_REF] Burghoff | Spectral analysis for a discrete metastable system driven by Lévy flights[END_REF] showed that similar scaling in ε for the spectral gap would hold if we were to restrict the SDE dynamics to a discrete grid with a small enough grid size.

The main result

In this chapter, we consider a stochastic differential equation with both a Brownian term and a Lévy term, and its Euler discretization as follows [START_REF] Duan | An Introduction to Stochastic Dynamics[END_REF]:

dw t = -∇f (w t )dt + εσdB t + εdL α t (6.
2) (6.3) with independent and identically distributed (i.i.d.) variables Z k ∼ N (0, I) where I is the identity matrix, the components of S k are i.i.d with SαS(1) distribution, and ε is the amplitude of the noise. Here, we choose σ as a scalar for convenience. We also note that the participation of the Brownian term in (6.2) enables the use of a Girsanov-like change of measures, that is convenient for our analysis.

w k+1 = w k -η∇f (w k ) + εση 1/2 Z k+1 + εη 1/α S k+1 ,
We formally define the first exit times, respectively for w t and w k as follows:

τ ξ,a (ε) inf{t ≥ 0 : w t -w ∈ [0, a + ξ]}, (6.4) τξ,a (ε) inf{k ∈ N : w k -w ∈ [0, a + ξ]}. (6.5)
where w is some local minimum and the processes are initialized at w 0 ≡ w 0 such that

w 0 -w ∈ [0, a].
Understanding the metastability behavior of SGD modeled by these dynamics requires understanding the first exit times for the continuous-time process w t given by (6.2) and its discretization w k (6.3). Similar to (4.12), we will study the first exit times defined by (6.4) and (6.5). Note that in the special case ξ = 0, we recover τ 0,a (ε) = τ a (ε) introduced previously in (4.12).

Our main goal is to obtain an explicit condition on the step-size, such that the first exit time of the continuous-time process τ ξ,a (ε) (6.4) would be similar to the first exit time of its Euler discretization τξ,a (ε) (6.5). In equations (6.2) and ( 6.3), σ is chosen as a scalar for convenience; however, we believe that this analysis can be extended to the case where σ is a function of w t .

Let us now state the assumptions which will imply our result.

A 8. The SDE (6.2) admits a unique strong solution.

A 9. Consider the process d ŵt = g( ŵ)dt + εσdB t + εdL α t , where ŵ ≡ { ŵt } t≥0 denotes the whole process and the drift g is defined as follows1 :

g( ŵ) - ∞ k=0 ∇f ( ŵkη )I [kη,(k+1)η) (t).
Here, I denotes the indicator function, i.e.

I S (x) = 1 if x ∈ S and I S (x) = 0 if x / ∈ S.
Then, the process φ t -g( ŵ)+∇f ( ŵt)

εσ satisfies: E exp 1 2 T 0 φ 2 t dt < ∞.
A 10. The gradient of f is γ-Hölder continuous: There exists a constant

M > 0 such that ∇f (x) -∇f (y) ≤ M x -y γ , ∀x, y ∈ R d .
A 11. The gradient of f satisfies the following assumption: ∇f (0) ≤ B.

A 12. For some m > 0 and b ≥ 0, f is (m, b, γ)-dissipative:

x, ∇f (x) ≥ m x 1+γ -b, ∀x ∈ R d .
We note that, as opposed to the theory of SDEs driven by Brownian motion, the theory of Lévy-driven SDEs is still an active research field where even the existence of solutions with general drift functions is not well-established and the main contributions have appeared in the last decade [START_REF] Priola | Pathwise uniqueness for singular SDEs driven by stable processes[END_REF], [START_REF] Kulik | On weak uniqueness and distributional properties of a solution to an SDE with α-stable noise[END_REF]. For this reason, A8

has been a common assumption in stochastic analysis, e.g. [START_REF] Imkeller | First exit times of SDEs driven by stable Lévy processes[END_REF], Imkeller et al. [2010a], [START_REF] Liang | Gradient estimates and ergodicity for SDEs driven by multiplicative Lévy noises via coupling[END_REF]. Nevertheless, existence and uniqueness results have been very recently established in [START_REF] Kulik | On weak uniqueness and distributional properties of a solution to an SDE with α-stable noise[END_REF] for SDEs with bounded Hölder drifts. Therefore A8 and A9 directly hold for bounded gradients and extending this result to Hölder and dissipative drifts is out of the scope of this study.

On the other hand, the assumptions A10-A12 are standard conditions, which are often considered in non-convex optimization algorithms that are based on discretization of diffusions [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF], [START_REF] Xu | Global convergence of langevin dynamics based algorithms for nonconvex optimization[END_REF], [START_REF] Erdogdu | Global non-convex optimization with discretized diffusions[END_REF], Gao et al. [2018b,a], S ¸imşekli et al. [2018], [START_REF] Liutkus | Sliced-Wasserstein flows: Nonparametric generative modeling via optimal transport and diffusions[END_REF].

The next assumption identifies an explicit condition for the step-size, which is required to make sure that the discrete process well-approximates the continuous one.

A 13. For a given δ > 0, t = Kη, and for some C > 0, the step-size satisfies the following condition:

0 < η ≤ min 1, m M 2 , δ 2 2K 1 t 2 1 γ 2 +2γ-1 , δ 2 2K 2 t 2 1 2γ , δ 2 2K 3 t 2 α 2γ , δ 2 2K 4 t 2 1 γ
, where ε is as in (6.3), the constants m, M, b are defined by A10-A12 and

K 1 = O(dε 2γ 2 -2 ), K 2 = O(ε -2 ), K 3 = O(d 2γ ε 2γ-2 ), K 4 = O(d 2γ ε 2γ-2 ).
More explicit forms of the constants are provided in appendix.

The main result of this chapter is presented in the following theorem, its proof can be found in appendix.

Theorem 11. Under assumptions A8-A13, the following inequality holds:

P[τ -ξ,a (ε) > Kη] -C K,η,ε,d,ξ -δ ≤ P[τ 0,a (ε) > K] ≤ P[τ ξ,a (ε) > Kη] + C K,η,ε,d,ξ + δ, where, C K,η,ε,d,ξ C 1 (Kη(dε + 1) + 1) γ e M η M η ξ + 1 -1 -Cde -ξ 2 e -2M η (εσ) -2 /(16dη) K + 1 -1 -C α d 1+α/2 ηe αM η ε α ξ -α K ,
for some constants C 1 , C α and C that does not depend on η or ε; M is given by A10 and ε is as in (6.2)-(6.3).

Remark Our result (Theorem 11) shows that with sufficiently small discretization step η (as in assumption A13), the probability to exit a given neighborhood of the local optimum at a fixed time t of the discretization process approximates that of the continuous process. This result also provides an explicit condition for the step-size, which explains certain impacts of the other parameters of the problem, such as dimension d, noise amplitude ε, variance of Gaussian noise σ, towards the similarity of the discretization and continuous processes. Theorem 11 enables the use of the metastability results for Lévy-driven SDEs for their discretized counterpart.

Exit time versus problem parameters. In Theorem 11, if we let η go to zero for any δ fixed, the constant C K,η,ε,d,ξ will also go to zero, and since δ can be chosen arbitrarily small, this implies that the probability of the first exit time for the discrete process and the continuous process will approach each other when the step-size gets smaller, as expected. If instead, we decrease d or ε, the quantity C K,η,ε,d,ξ also decreases monotonically, but it does not go to zero due to the first term in the expression of C K,η,ε,d,ξ .

Exit time versus width of local minima. Popular activation functions used in deep learning such as ReLU functions are almost everywhere differentiable and therefore the cost function has a well-defined Hessian almost everywhere (see e.g. Li and Yuan [2017]).

The eigenvalues of the Hessian of the objective near local minima have also been studied in the literature (see e.g. [START_REF] Sagun | Eigenvalues of the Hessian in deep learning: Singularity and beyond[END_REF], [START_REF] Vardan Papyan | The full spectrum of deep net Hessians at scale: Dynamics with sample size[END_REF]). If the Hessian around a local minimum is positive definite, the conditions for the multi-dimensional version of Theorem 4 in Imkeller et al. [2010a]) are satisfied locally around a local minimum. For local minima lying in wider valleys, the parameter a can be taken to be larger; in which case the expected exit time Eτ 0,a (ε) ∼ O(a α ) will be larger by the formula (4.14). In other words, the SDE (4.6) spends more time to exit wider valleys. Theorem 11 shows that SGD modeled by the discretization of this SDE will also inherit a similar behavior if the step-size satisfies the conditions we provide (A13).

Proof overview

Relating the first exit times for w t (6.2) and w k (6.3) often requires obtaining bounds on the distance between w kη and w k . For this purpose, for any given local minimum w of f and a > 0, we define the following set (6.6) which contains the sets of K points in R d , each point at a distance of at most a from the local minimum w. We also define the following set

A (w 1 , . . . , w K ) ∈ R d × . . . × R d : max k≤K w k -w ≤ a ,
N a w ∈ R d : w -w ≤ a . (6.7)
If w kw kη is small with high probability, then we expect that their first exit times from the set N a will be close to each other as well with high probability. For objective functions with bounded gradients, in order to relate τ ξ,a (ε) to τξ,a (ε), one can attempt to use the strong convergence of the Euler scheme (cf. [START_REF] Mikulevičius | On the rate of convergence of strong Euler approximation for SDEs driven by lévy processes[END_REF] Proposition 1): lim η→0 E w kw kη = 0. By using Markov's inequality, this result implies convergence in probability: for any δ > 0 and ε > 0, there exists η such that P( w kw kη > ε) < δ/2. Then, if w kη ∈ N a one of the following events must happen:

1. w k ∈ N a ,
2. w k ∈ N a and w kw kη > (with probability less than δ/2), 3. w k / ∈ N a and distance from w k to N a is at most ε (with probability less than δ/2).

By using this observation, we obtain:

P[w kη ∈ N a ] ≤ P[w k ∈ N a ] + δ.
Even though we could use this result in order to relate τ ξ,a (ε) to τξ,a (ε), this approach would not yield a meaningful condition for η since the bounds for the strong error E w kw kη often grows exponentially in general with k, which means η should be chosen exponentially

The proof technique is similar to the approach of Dalalyan [2017b], Raginsky et al.

[2017]: the idea is to divide the integral in (6.9) into smaller pieces and bounding each piece separately. Once we obtain a bound on KL, by using an optimal coupling argument, the data processing inequality, and Pinsker's inequality, we obtain a bound for the total variation (TV) distance between µ t and μt as follows:

P M [(w η , . . . , w Kη ) = ( ŵη , . . . , ŵKη )] ≤ µ Kη -μKη T V ≤ 1 2 KL(μ Kη , µ Kη ) 1 2 .
where M denotes the optimal coupling between {w s } s∈[0,Kη] and { ŵs } s∈[0,Kη] , i.e., the joint probability measure of {w s } s∈[0,Kη] and { ŵs } s∈[0,Kη] , which satisfies the following identity [START_REF] Lindvall | Lectures on the coupling method[END_REF]:

P M [{w s } s∈[0,Kη] = { ŵs } s∈[0,Kη] ] = µ Kη -μKη T V .
Combined with Theorem 12, this inequality implies the following useful result:

P[(w η , . . . , w Kη ) ∈ A] -δ ≤ P[τ 0,a (ε) > K] ≤ P[(w η , . . . , w Kη ) ∈ A] + δ (6.10)
where we used the fact that the event ( ŵ(η), . . . , ŵ(Kη)) ∈ A is equivalent to the event (τ 0,a (ε) > K). The remaining task is to relate the probability P[(w η , . . . , w Kη ) ∈ A]

to P[τ ξ,a (ε) > Kη]. The event (w η , . . . , w Kη ) ∈ A ensures that the process w t does not leave the set N a when t = η, . . . , Kη; however, it does not indicate that the process remains in N a when t ∈ (kη, (k + 1)η). In order to have a control over the whole process, we introduce the following event:

B max 0≤k≤K-1 sup t∈[kη,(k+1)η] w t -w kη ≤ ξ ,
such that the event [(w η , . . . , w Kη ) ∈ A] ∩ B ensures that the process stays close to N a for the whole time. By using this event, we can obtain the following inequalities:

P[(w η , . . . , w Kη ) ∈ A] ≤P[(w η , . . . , w Kη ) ∈ A ∩ B] + P[(w η , . . . , w Kη ) ∈ B c ] =P[τ ξ,a (ε) > Kη] + P[(w η , . . . , w Kη ) ∈ B c ].
By using the same approach, we can obtain a lower bound on P[(w η , . . . , w Kη ) ∈ A] as well. Hence, our final task reduces to bounding the term P[(w η , . . . , w Kη ) ∈ B c ], which we perform by using the weak reflection principles of Lévy processes Bayraktar et al.

[2015]. This finally yields Theorem 11.

Numerical illustration

To illustrate our results, we perform the experiments on a synthetic setting and the experiments on real data: a multi-layer fully connected neural network with ReLu activations on the MNIST dataset. 5 ×10

-7

α = 1.2 α = 1.4 α = 1.6 α = 1.8 (d)

Synthetic data

We first conduct experiments on a synthetic problem, where the cost function is set to f (x) = 1 2 x 2 . This corresponds to an Ornstein-Uhlenbeck-type process, which is commonly considered in metastability analyses [START_REF] Duan | An Introduction to Stochastic Dynamics[END_REF]. This process locally satisfies the conditions A8-A12.

Since we cannot directly simulate the continuous-time process, we consider the stochastic process sampled from (6.3) with sufficiently small step-size as an approximation of the continuous scheme. Thus, we organize the experiments as follows. We first choose a very small step-size, i.e. η = 10 -10 . Starting from an initial point w 0 satisfying w 0 < a, we iterate (6.3) until we find the first K such that w K > a. We repeat this experiment 100 times, then we take the average Kη as the 'ground-truth' first exit time. We continue the experiments by calculating the first exit times for larger step-sizes (each repeated 100 times), and compute their distances to the ground truth.

The detailed settings of the parameters for this experiment (Figure 6.1) can be found in appendix.

By Theorem 11, the distance between the first exit times of the discretization and the continuous processes depends on two terms C K,η,ε,d,ε and δ, which are used for explaining our experimental results.

We observe from Figure 6.1(a) that the error to the ground-truth first exit time is an increasing function of η, which directly matches our theoretical result. Figure 6.1(b)

shows that, with small noise limit (e.g., in our settings, ε < 1 versus η ≈ 10 -8 ), the error decreases with the parameter ε. By A13, with increased ε, we have the term δ to be reduced. On the other hand, C K,η,ε,d,ε increases with ε. However, at small noise limit, this effect is dominated by the decrease of δ, that makes the error decrease overall. The decreasing speed then decelerates with larger ε, since, the product εη becomes so large that the increase of C K,η,ε,d,ε starts to dominate the decrease of δ. Thus, it suggests that for a large ε, a very small step-size η would be required for reducing the distance between the first exit times of the processes. In Figure 6.1(c), the error decreases when the variance σ increases. The reason for the performance is the same as in (b), and can be explained by considering the expression of δ and C K,η,ε,d,ε in the conclusion of Theorem 11.

In Figure 6.1(d), for small dimension, with the same exit time interval, when we increase d, both processes escape the interval earlier, with smaller exit times. Hence, the distance between their exit times becomes smaller. With larger d, the increasing effect of δ and C K,η,ε,d,ε starts to dominate the above 'early-escape' effect, thus, the decreasing speed of the error diminish. We observe that the error even slightly increases when α = 1.2 and d grows from 70 to 100. In our second set of experiments, we consider the real data setting used in S ¸imşekli et al. [2019]: a multi-layer fully connected neural network with ReLu activations on the MNIST dataset. We adapted the code provided in S ¸imşekli et al. [2019] and we provide our version in https://github.com/umutsimsekli/sgd_first_exit_time. For this model, we followed a similar methodology: we monitored the first exit time by varying the η, the number of layers (depth), and the number of neurons per layer (width). Since a local minimum is not analytically available, we first trained the networks with SGD until a vicinity of a local minimum is reached with at least 90% accuracy, then we measured the first exit times with a = 1 and ε = 0.1. In order to have a prominent level of gradient noise, we set the minibatch size b = 10 and we did not add explicit Gaussian or Lévy noise. The result is given in Figure 6.2. We observe that, even with pure gradient noise, the error in the exit time behaves very similarly to the one that we observed in Figure 6.1(a), hence supporting our theory. We further observe that, the error has a better dependency when the width and depth are relatively small, whereas the slope of the error increases for larger width and depth. This result shows that, to inherit the metastability properties of the continuous-time SDE, we need to use a smaller η as we increase the size of the network. Note that this result does not conflict with Figure 6.1(d), since changing the width and depth does not simply change d, it also changes the landscape of the problem.

Chapter 7

Conclusion and future work

In this thesis, we investigated the tail behavior of the gradient noise in deep neural networks and empirically showed that the gradient noise is highly non-Gaussian. We analyzed the global convergence property of SGD for non-convex optimization via a stochastic process, which can be seen as a perturbed version of the gradient descent algorithm with heavy-tailed α-stable noise, for non-convex optimization and proved finite-time bounds for its expected suboptimality. Our results agreed with the existing related work, and showed that the weak-error of this algorithm increases faster than unadjusted Langevin algorithm (ULA), which suggests using smaller step-sizes.

In terms of metastability, we derived explicit conditions for the step-size such that the discrete-time SGD can inherit the metastability behavior of its continuous-time limit.

These outcomes enabled us to analyze SGD as a stochastic differential equation (SDE) driven by a Lévy motion and establish a bridge between SGD and existing theoretical results, which provides more insights on the behavior of SGD, especially in terms of choosing wide minima.

Future directions

Our study also brings up the following questions:

1. We observe that the tail-index might depend on the current state w k , which suggests analyzing SGD as a 'stable-like process' where the tail-index can depend on time [START_REF] Bass | Uniqueness in law for pure jump Markov processes[END_REF]. However, the metastability behavior of these processes are not clear at the moment and its theory is still in an early phase [START_REF] Kuhwald | Bistable behaviour of a jump-diffusion driven by a periodic stable-like additive process[END_REF].

2. At the initial point, in the over-parametrized regime with large batch sizes, the noise can in fact be of Gaussian nature (Figure 4.2). However, this property is destroyed quickly (see [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF], [START_REF] Der | Beyond Gaussian processes: On the distributions of infinite networks[END_REF], [START_REF] Lee | Deep neural networks as Gaussian processes[END_REF] for a discussion on the infinite width networks, and [START_REF] Panigrahi | Non-Gaussianity of stochastic gradient noise[END_REF] for a discussion on the early phases and large batches). We note that such Gaussianity heavily depends on the structure of the data, initialization scheme, and the size of the network in a sensitive way and may hold in only certain regimes or in specific cases. We think that identifying the crossover between the Gaussian and non-Gaussian regimes depending on the architecture and data is an important open problem.

3. Even though the general heavy-tailed behavior remains unchanged with the choice of the loss function, we still observe different behaviors in terms of relation to generalization (see Figure 4.11). We note that our results are related to the findings of [START_REF] Martin | Traditional and heavy-tailed self regularization in neural network models[END_REF], which modeled the weight matrices as heavy-tailed random matrices and investigated the density of the singular values of those matrices.

Their empirical results on various different types of neural networks show that when the batch size gets smaller, the training process is able to catch finer-scale correlations from the data, leading to more strongly-correlated models between the layers of the network and that the entries of the weight matrices and the density of its singular values have heavier tails. Our results in Section 4.6.4 are partially consistent with the findings of [START_REF] Martin | Traditional and heavy-tailed self regularization in neural network models[END_REF]. Their results combined with ours would shed more light into the heavy-tailedness of the SGD iterates and generalization properties of SGD algorithms. In the future, we would like to investigate the underlying deeper connections between the heavy-tailed behavior and generalization further from both a mathematical and experimental perspectives. 4. We have empirically observed a heavy-tailed behavior in the stochastic gradient noise; however, it is still not (rigorously) clear what the underlying mechanism that drives this heavy-tailed behavior is. One possible idea is to consider a deep network with l layers f (w) w (l) σ l (w (l-1) σ l-1 (...w (2) σ 2 (w (1) x i ) . . .)) where we look for optimal weights w by minimizing the associated loss function using SGD algorithm, and then relate the stochastic gradient noise to a self-similar process which has a close connection with Lévy-driven SDE and α-stable distribution [START_REF] Pipiras | Long-range dependence and self-similarity[END_REF]. Additionally, there have been nice theoretical guarantees established for those self-similar processes [START_REF] Pipiras | Long-range dependence and self-similarity[END_REF] which could be useful for building a theoretical analysis of the heavy-tailed behavior of the stochastic gradient noise. Therefore, investigating this underlying mechanism would be a promising future direction. 5. Another interesting research direction concerns the multifractality, which is one of the important aspects of stochastic processes such as fractional Brownian motion, Lévy motion, etc. These processes are used for modeling phenomena in nature [START_REF] Bobrov | Multifractal analysis: a method to investigate non-stationary properties of geophysical processes[END_REF], [START_REF] Pavlov | Characterizing scaling properties of complex signals with missed data segments using the multifractal analysis[END_REF], [START_REF] Teotia | Role of multifractal analysis in understanding the preparation zone for large size earthquake in the North-Western Himalaya region[END_REF], which are often very irregular.

Multifractal analysis consists of studying the set of the irregular points, characterized by the Hölder exponent Yang et al. [2018], of a function f . Motivated by the first exit time problem from local minima, it is also intriguing to study the connection between FLA (3.3) and the (ir)regularity of function f by using multifractality: The probability that w k in (3.3) resides in a set of low regularity points of f and the role of index α towards this probability. 6. Besides first exit time and metastability, studying the minimal time needed for a stochastic optimization algorithm (e.g. equation (4.2)) to enter a region containing a local minimum is important as well. Such minimal time is called hitting time and has been studied in recent years Zhang et al. [2017b], [START_REF] Chen | On stationary-point hitting time and ergodicity of stochastic gradient Langevin dynamics[END_REF] for Stochastic Gradient Langevin Dynamics, which is obtained by replacing the gradient in ULA (3.2) by stochastic gradient. Thus, extending these results for FLA (3.3) is a natural research direction and would provide a better understanding of the heavytailed behavior in the stochastic gradient noise.

D α i q(x, t).
Using Definition (5.6) of b, we have

∂ t q(x, t) = - d i=1 ∂ i β -1 D α-2 i (-βφ(x)∂ i f (x)) φ(x) q(x, t) -β -1 d i=1 D α i q(x, t) = - d i=1 ∂ i β -1 D α-2 i (-βπ(x)∂ i f (x)) π(x) q(x, t) -β -1 d i=1 D α i q(x, t) = - d i=1 ∂ i β -1 D α-2 i (∂ i π(x)) π(x) q(x, t) -β -1 d i=1 D α i q(x, t).
Here, we used π(x) = φ(x)/ φ(x)dx in the second equality and

-β∂ i f (x) = ∂ i log π(x) = ∂iπ(x)
π(x) in the third equality. Next, by replacing q by π on the right hand side of the above equality, we have: 77

- d i=1 ∂ i β -1 D α-2 i (∂ i π(x)) π(x) π(x, t) -β -1 d i=1 D α i π(x, t) = = - d i=1 ∂ i β -1 D α-2 i (∂ i π(x)) -β -1 d i=1 D α i π(x, t) = - d i=1 ∂ i ∂ i β -1 D α-2 i (π(x)) -β -1 d i=1 D α i π(x, t) = d i=1 D 2 i [β -1 D α-2 i (π(x))] -β -1 d i=1 D α i π(x, t) = d i=1 D α i [β -1 π(x)] -β -1 d i=1 D α i π(x, t) =0.
Here, we used Proposition 1 in S ¸imşekli [2017], D 2 u(x) = -∂ ∂x 2 u(x), and the semigroup property of the Riesz derivation D a D b u(x) = D a+b u(x). This proves that π is an invariant measure of the Markov process (x(t)) t≥0 .

Proof of Proposition 1

Proof. By Corollary 1.3 in [START_REF] Liang | Gradient estimates and ergodicity for SDEs driven by multiplicative Lévy noises via coupling[END_REF], the assumptions imply that there exist constants C > 0 and C1 > 0 such that W 1 (µ 3t , π) ≤ Cβe -C1t . Let P 3t be the coupling of µ 3t and π that such that W 1 (µ 3t , π) =

x 3 (t) -ŵ dP 3t .

For 0 < λ < 1, by Hölder inequality,

W λ λ (µ 3t , π) ≤ x 3 (t) -ŵ λ dP 3t ≤ x 3 (t) -ŵ dP 3t λ =W λ 1 (µ 3t , π) For α > λ > 1, W λ λ (µ 3t , π) ≤ x 3 (t) -ŵ λ dP 3t ≤ x 3 (t) -ŵ δ x 3 (t) -ŵ λ-δ dP 3t ≤ x 3 (t) -ŵ dP 3t δ x 3 (t) -ŵ (λ-δ)/(1-δ) dP 3t 1-δ =W δ 1 (µ 3t , π) x 3 (t) -ŵ (λ-δ)/(1-δ) dP 3t 1-δ ,
where we used Holder's inequality for δ < 1 such that (λ-δ)/(1-δ) < α, and x 3 (t)ŵ (λ-δ)/(1-δ) dP 3t is bounded by a constant, by assumption A5.

Finally, we have

W λ (µ 3t , π) ≤Cβe -C1t ,
for some constants C, C 1 > 0 and for 0 < λ < α. This completes the proof.

Remark Let us consider the case where the dimension d is equal to 1 (the extension for d > 1 is similar). The first part of assumption A5 can be satisfied under the following (rather non-trivial) assumptions. Assume that there exist constants

P, C 1 , C 2 , C 3 , C 4 , C 5 , C 6 > 0 such that: f (z) > 0 if z > P, (7.1) |z|≤P |φ(z)f (z)|dz = C 1 > 0 (7.2) z<-P φ(z)|f (z)||z| 1-α dz = C 2 > 0 (7.3) z>P φ(z)f (z)|z| 1-α dz = C 3 > 0, (7.4) if |z| ≤ P : 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 ≤ C 4 |x -y| ∀x, y ∈ R, (7.5) if z < -P : 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 ≤ C 5 |x -y||z| 1-α ∀x, y ∈ R, (7.6) if z > P : 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 ≤ C 6 |z| 1-α (y -x) ∀x, y ∈ R s.t x > y, (7.7) C 1 C 4 + C 2 C 5 < C 3 C 6 . (7.8)
By definition of Riesz potential, we have:

b(x) -b(y) = R φ(z)f (z) φ(x)|x -z| α-1 dz - R φ(z)f (z) φ(y)|y -z| α-1 dz = R φ(z)f (z) 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz = |z|≤P φ(z)f (z) 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz + z<-P φ(z)f (z) 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz + z>P φ(z)f (z) 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz.
By these assumptions, we estimate the first term on the right hand side in the above expression of b(x) -b(y), for x > y, as follows:

|z|≤P φ(z)f (z) 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz ≤ ≤ |z|≤P |φ(z)f (z)| 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz ≤ |z|≤P |φ(z)f (z)|C 4 |x -y|dz =C 1 C 4 |x -y| =C 1 C 4 (x -y).
For the remaining terms, we have:

z<-P φ(z)f (z) 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz ≤ ≤ z<-P |φ(z)f (z)| 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz ≤ z<-P φ(z)|f (z)|C 5 |z| 1-α |x -y|dz =C 2 C 5 |x -y| =C 2 C 5 (x -y),
and z>P φ(z)f (z) 1 φ(x)|x -z| α-1 - 1 φ(y)|y -z| α-1 dz ≤ ≤ z>P φ(z)f (z)C 6 |z| 1-α (y -x)dz =C 3 C 6 (y -x) = -C 3 C 6 (x -y).
By combining these estimates, we get, for x > y: b Proof of Lemma 2 In this section, we precise the statement of Lemma 2 and provide the proof. Lemma 6. Let v and w be two random variables on R d which have µ and ν as the probability measures and let g be a function in C 1 (R d , R). Assume that for some c 1 > 0, c 2 ≥ 0 and 0 ≤ γ < 1,

(x) -b(y) ≤ (C 1 C 4 + C 2 C 5 -C 3 C 6 )(x -y). Thus, (b(x) -b(y))(x -y) ≤ (C 1 C 4 + C 2 C 5 -C 3 C 6 )(x -y) 2 . Since C 1 C 4 + C 2 C 5 - C 3 C 6 < 0,
∇g(w) ≤ c 1 w γ + c 2 , ∀w ∈ R d
then the following bound holds:

gdµ -gdν ≤ c 1 E w γp 1 p + c 1 E v γp 1 p + c 2 W q (µ, ν).
Proof. We have

g(v) -g(w) = 1 0 w -v, ∇g((1 -t)v + tw) dt ≤ 1 0 w -v ∇g((1 -t)v + tw) dt (by Cauchy-Schwarz) ≤ 1 0 w -v (c 1 ((1 -t) v + t w ) γ + c 2 )dt (by the assumption on ∇g) ≤ w -v c 1 ( v + w ) γ + c 2 ≤ w -v (c 1 v γ + c 1 w γ + c 2 ).
(by Lemma 16)

Now let P be a joint probability distribution of µ and ν that achieves W λ (µ, ν), that is, P = L((w, v)) with µ = L(w) and ν = L(v). We have

gdµ -gdν = E P [g(w) -g(v)] ≤ [E P (c 1 w γ + c 1 v γ + c 2 ) p ] 1 p [E P w -v q ] 1 q ≤ c 1 E P w γp 1 p + c 1 E P v γp 1 p + c 2 W q (µ, ν),
where we have used Holder's inequality and Minkowski's inequality.

Proof of Lemma 3

Proof. We define a real function F λ as follows:

F λ (y) y λ . (7.9) It is clear that F λ is a C 1 function. Let y(t) x 1 (t) -x 2 (t). By the chain rule, dF λ (y(t)) = ∇F λ (y(t)), b 1 (x 1 (t-), α) -b 2 (x 2 (t-), α) dt = λ x 1 (t) -x 2 (t)) λ-2 x 1 (t) -x 2 (t), b 1 (x 1 (t-), α) -b 2 (x 2 (t-), α) dt.
(7.10)

By integrating both sides of (7.10) with respect to t, we arrive at

F λ (y(t)) =F λ (y(0)) + t 0 λ x 1 (t) -x 2 (t)) λ-2 x 1 (t) -x 2 (t), b 1 (x 1 (t-), α) -b 2 (x 2 (t-), α) ds = t 0 λ x 1 (t) -x 2 (t)) λ-2 x 1 (t) -x 2 (t), b 1 (x 1 (t-), α) -b 2 (x 2 (t-), α) ds.
By definition of Wasserstein distance, we have

W λ (µ 1t , µ 2t ) = inf{(E[F λ (y(t))]) 1/λ },
which is the desired result.

Proof of Theorem 8 In this section, we first precise the statement of Theorem 8 and then provide the corresponding proof.

Theorem 13. Let E L α (1) λ l α,λ,d < ∞. We also define the following quantities:

P 1 (η) cη d β 1/α 1 p 1 + (cη) 1 p 1 + (2η(b + m)) (q-1) 2
+ 2

(q-1) 2

(ηB) (q-1)

+ η β (q-1) α l 1 p 1
α,(q-1)p1,d + η q-1 M q-1 (2η(b + m))

(q-1)γ 2

+ 2

(q-1)γ 2

(ηB) (q-1)γ

+ η β (q-1)γ α l 1 p 1 α,(q-1)p1γ,d , P 2 (η) M cη d β 1/α 1 q 1 + (cη) 1 q 1 + (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 q 1 α,γq1,d , Q 1 (η) c 1 p 1 + (E x 2 (0) (q-1)p1 ) 1 p 1 + η q-1 M q-1 (E x 2 (0) (q-1)p1γ ) 1 p 1 + B (q-1) + η β q-1 α l 1 p 1 α,(q-1)p1,d , Q 2 M (E x 2 (0) γq1 ) 1 q 1 + M c 1 q 1 .
Under additional assumption on the step-size: 0 < η ≤ m M 2 , we have

W q q (µ 1t , µ 2t ) ≤ qη k 2 P 1 (η)P 2 (η) + k 1+1/p1 P 1 (η)Q 2 + k 1+1/q1 P 2 (η)Q 1 (η) + kQ 1 (η)Q 2 .
Proof. From Lemma 3, we have

W q q (µ 1t , µ 2t ) = =E t 0 q x 1 (s) -x 2 (s)) q-2 x 1 (s) -x 2 (s), b 1 (x 1 (s-), α) -b 2 (x 2 (s-), α) ds = k-1 j=0 E (j+1)η jη q x 1 (s) -x 2 (s)) q-2 x 1 (s) -x 2 (s), b 1 (x 1 (s-), α) -b 2 (x 2 (s-), α) ds ≤ k-1 j=0 E (j+1)η jη q x 1 (s) -x 2 (s) q-1 c α ∇f (x 1 (s)) -∇f (x 2 (jη)) ds =q k-1 j=0 (j+1)η jη E x 1 (s) -x 2 (s) q-1 c α ∇f (x 1 (s)) -∇f (x 2 (jη)) ds ≤q k-1 j=0 (j+1)η jη E x 1 (s) -x 2 (s) (q-1)p1 1 p 1 E c α (∇f (x 1 (s)) -∇f (x 2 (jη))) q1 1 q 1 ds, 82
where we have used Cauchy-Schwarz inequality in the third line and Holder's inequality in the last line.

Since (q -1)p 1 < 1 by assumption A4, using Lemma 16 twice, we have:

E x 1 (s) -x 2 (s) (q-1)p1 1 p 1 ≤ E x 1 (s) (q-1)p1 + E x 2 (s) (q-1)p1 1 p 1 ≤ E x 1 (s) (q-1)p1 1 p 1 + E x 2 (s) (q-1)p1 1 p 1 .
Then, by applying Lemma 9 and Lemma 12 for s ∈ [jη, (j + 1)η), we obtain:

E x 1 (s)-x 2 (s) (q-1)p1 1 p 1 ≤ ≤ c s d β 1/α + 1 + 1 q-1 + E x 2 (0) (q-1)p1 + j (2η(b + m)) (q-1)p 1 2 + 2
(q-1)p 1 2

(ηB) (q-1)p1 + η β (q-1)p 1 α l α,(q-1)p1,d + (s -jη) (q-1)p1 M (q-1)p1 E x 2 (0) (q-1)p1γ + j (2η(b + m))

(q-1)p 1 γ 2 + 2

(q-1)p 1 γ 2

(ηB) (q-1)p1γ + η β (q-1)p 1 γ α l α,(q-1)p1γ,d + B (q-1)p1 + sjη β (q-1)p 1 α l α,(q-1)p1,d

1 p 1 .
Next, using Lemma 16, the inequalities j < j + 1 and sjη ≤ η for s ∈ [jη, (j + 1)η),

we get E x 1 (s)-x 2 (s) (q-1)p1

1 p 1 ≤ ≤ c s d β 1/α + 1 + 1 q-1 + (E x 2 (0) (q-1)p1 ) 1 p 1 + (j + 1) 1 p 1 (2η(b + m)) (q-1) 2
+ 2

(q-1) 2

(ηB) (q-1) + η β (q-1) α l 1 p 1 α,(q-1)p1,d + η q-1 M q-1 (E x 2 (0) (q-1)p1γ )

1 p 1 + (j + 1) 1 p 1 (2η(b + m)) (q-1)γ 2 + 2 (q-1)γ 2 (ηB) (q-1)γ + η β (q-1)γ α l 1 p 1 α,(q-1)p1γ,d + B (q-1) + η β q-1 α l 1 p 1
α,(q-1)p1,d .

We note that s < (j + 1)η and q -1 < 1 p1 (from the assumptions). Hence,

c s d β 1/α + 1 + 1 q-1 ≤ c (j + 1)η d β 1/α + 1 + 1 1 p 1 ≤(j + 1) 1 p 1 cη d β 1/α + 1 1 p 1 + c 1 p 1 ,
where the last inequality is an application of Lemma 16. By replacing this inequality into the previous one and rearranging the terms, we have

E x 1 (s) -x 2 (s) (q-1)p1 1 p 1 ≤c 1 p 1 + (E x 2 (0) (q-1)p1 ) 1 p 1 + η q-1 M q-1 (E x 2 (0) (q-1)p1γ ) 1 p 1 + B (q-1) + η β q-1 α l 1 p 1
α,(q-1)p1,d + (j + 1)

1 p 1 cη d β 1/α + 1 1 p 1 + (2η(b + m)) (q-1) 2
+ 2

(q-1) 2

(ηB) (q-1) + η β

(q-1) α l 1 p 1
α,(q-1)p1,d + η q-1 M q-1 (2η(b + m))

(q-1)γ 2

+ 2

(q-1)γ 2

(ηB) (q-1)γ + η β

(q-1)γ α l 1 p 1 α,(q-1)p1γ,d ≤c 1 p 1 + (E x 2 (0) (q-1)p1 ) 1 p 1 + η q-1 M q-1 (E x 2 (0) (q-1)p1γ ) 1 p 1 + B (q-1) + η β q-1 α l 1 p 1
α,(q-1)p1,d + (j + 1)

1 p 1 cη d β 1/α 1 p 1 + (cη) 1 p 1 + (2η(b + m)) (q-1) 2 + 2 (q-1) 2 (ηB) (q-1) + η β (q-1) α l 1 p 1
α,(q-1)p1,d + η q-1 M q-1 (2η(b + m))

(q-1)γ 2

+ 2

(q-1)γ 2

(ηB) (q-1)γ + η β

(q-1)γ α l 1 p 1 α,(q-1)p1γ,d
=Q 1 (η) + (j + 1)

1 p 1 P 1 (η).
Here, we have used Lemma 16 in the last inequality. Now, consider the following quantity

E c α (∇f (x 1 (s)) -∇f (x 2 (jη))) q1 1 q 1 ≤ E M x 1 (s) -x 2 (jη) γ q1 1 q 1 ≤ E M x 1 (s) γ + M x 2 (jη) γ q1 1 q 1 ≤ E M q1 x 1 (s) γq1 1 q 1 + E M q1 x 2 (jη) γq1 1 q 1 ,
where we have used assumption A2, Lemma 16 and Minkowski's inequality.

By Lemma 9 and Lemma 12, we have

E c α ∇f (x 1 (s))-c α ∇f (x 2 (jη)) q1 1 q 1 ≤ ≤M c s d β 1/α + 1 + 1 γ + M q1 (E x 2 (0) γq1 ) + M q1 j (2η(b + m)) γq 1 2 + 2 γq 1 2 (ηB) γq1 + η β γq 1 α l α,γq1,d 1 q 1 .
By using Lemma 16 and the inequality j < j + 1, we have

E c α ∇f (x 1 (s)) -c α ∇f (x 2 (jη)) q1 1 q 1 ≤ ≤M c s d β 1/α + 1 + 1 γ + M (E x 2 (0) γq1 ) 1 q 1 + M (j + 1) 1 q 1 (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 q 1 α,γq1,d .

84

We note that s < (j + 1)η and γ < 1 q1 (from the assumptions). Hence,

c s d β 1/α + 1 + 1 γ ≤ c (j + 1)η d β 1/α + 1 + 1 1 q 1 ≤(j + 1) 1 q 1 cη d β 1/α + 1 1 q 1 + c 1 q 1 ,
where the last inequality is an application of Lemma 16. By replacing this inequality into the previous one and rearranging the terms, we have

E c α ∇f (x 1 (s)) -c α ∇f (x 2 (jη)) q1 1 q 1 ≤ ≤M (E x 2 (0) γq1 ) 1 q 1 + M c 1 q 1 + M (j + 1) 1 q 1 cη d β 1/α + 1 1 q 1 + (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 q 1 α,γq1,d ≤M (E x 2 (0) γq1 ) 1 q 1 + M c 1 q 1 + M (j + 1) 1 q 1 cη d β 1/α 1 q 1 + (cη) 1 q 1 + (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 q 1 α,γq1,d
=Q 2 + (j + 1)

1 q 1 P 2 (η).
Here, we used Lemma 16 in the last inequality. By combining the above inequalities,

we get

E t 0 q x 1 (s) -x 2 (s)) q-2 x 1 (s) -x 2 (s), b 1 (x 1 (s-), α) -b 2 (x 2 (s-), α) ds ≤ ≤ k-1 j=0
qη (j + 1)P 1 (η)P 2 (η) + (j + 1)

1 p 1 P 1 (η)Q 2 + (j + 1) 1 q 1 P 2 (η)Q 1 (η) + Q 1 (η)Q 2 ≤ qη k 2 P 1 (η)P 2 (η) + k 1+1/p1 P 1 (η)Q 2 + k 1+1/q1 P 2 (η)Q 1 (η) + kQ 1 (η)Q 2 .
The final conclusion follows from this inequality.

Proof of Corollary 1

Proof. In order to get the results from the bound obtained by Theorem 13, we take the max power of k and the min power of η among the terms containing k and η but not containing β. For the terms containing β, we take the max power of k, min power of η, min power of 1/β and max power of d. We get W q q (µ 1t , µ 2t ) ≤C(k 2 η + k 2 η 1+min{γ,q-1}/α β -(q-1)γ/α d).

Since γ < 1/p = (q -1)/q < q -1, we finally obtain

W q q (µ 1t , µ 2t ) ≤C(k 2 η + k 2 η 1+γ/α β -(q-1)γ/α d).

Proof of Corollary 2

Proof. The proof starts from the bound established in Corollary 8 then, follows the same lines of the proof of Corollary 1.

Proof of Theorem 7

Proof. We have the decomposition:

E[f (w k )] -f * =E[f (x 2 (kη))] -f * =(E[f (x 2 (kη))] -E[f (x 1 (kη))]) + (E[f (x 1 (kη))] -E[f (x 3 (kη))]) + (E[f (x 3 (kη))] -E[f ( ŵ))]) + (E[f ( ŵ))] -f * ).
By Corollary 2, Corollary 4, Lemma 4 and Lemma 5, there exists a constant C independent of k, η and β such that

E[f (w k )] -f * ≤C k 1+ 1 q η 1 q + k 1+ 1 q η 1 q + γ αq β -(q-1)γ αq d + k γ+ γ+q q η γ+ 1 q β -γ α d + k γ+ γ+q q η 1 q + β b + d/β m exp(-λ * β -1 t) + β -γ-1 M c -1 α 1 + γ + β -1 log (2e(b + d/β)) d/2 Γ(d/2 + 1)β d (dm) d/2 .
Here, we note that kη = t. then by taking the largest power of k, smallest powers of η and β -1 among the terms containing all of three parameters k, η and β, there exist a constant C satisfying the following inequality:

E[f (w k )] -f * ≤C k 1+max{ 1 q ,γ+ γ q } η 1 q + k 1+max{ 1 q ,γ+ γ q } η 1 q + γ αq β -(q-1)γ αq d + β b + d/β m exp(-λ * β -1 kη) + β -γ-1 M c -1 α 1 + γ + β -1 log (2e(b + d/β)) d/2 Γ(d/2 + 1)β d (dm) d/2 .
Proof of Theorem 9 In this section, we precise the statement of Theorem 9 and provide the full proof.

Theorem 14. We have the following estimate:

W q q (µ 1t , µ 3t ) ≤qt M (c q-1 + c q-1 b )(c γ + c γ b ) t d β 1/α + 1 + 1 q-1+γ + L(c q-1 + c q-1 b ) t d β 1/α + 1 + 1 q-1
, where c and c b are constants defined in Lemma 9 and Lemma 10.

Proof. From Lemma 3, we have

W q q (µ 1t , µ 3t ) =E t 0 q x 1 (s) -x 3 (s)) q-2 x 1 (s) -x 3 (s), b 1 (x 1 (s-), α) -b(x 3 (s-), α) ds = t 0 q x 1 (s) -x 3 (s)) q-2 x 1 (s) -x 3 (s), b 1 (x 1 (s-), α) -b(x 3 (s-), α) ds ≤E t 0 q x 1 (s) -x 3 (s) q-1 c α ∇f (x 1 (s)) + b(x 3 (s), α) ds =q t 0 E x 1 (s) -x 3 (s) q-1 c α ∇f (x 1 (s)) + b(x 3 (s), α) ds ≤q t 0 E x 1 (s) -x 3 (s) (q-1)p1 1 p 1 E c α ∇f (x 1 (s)) + b(x 3 (s), α) q1 1 q 1 ds,
where we have used Cauchy-Schwarz inequality in the third line and Holder's inequality in the last line.

Since (q -1)p 1 < 1 by assumption A4, using Lemma 16 twice, we have:

E x 1 (s) -x 3 (s) (q-1)p1 1 p 1 ≤ E x 1 (s)| (q-1)p1 + E x 3 (s) (q-1)p1 1 p 1 ≤ E x 1 (s) (q-1)p1 1 p 1 + E x 3 (s) (q-1)p1 1 p 1 .
Then, by applying Lemma 9 and Lemma 10 we obtain:

E x 1 (s)-x 3 (s) (q-1)p1 1 p 1 ≤ c s d β 1/α + 1 + 1 q-1 + c b s d β 1/α + 1 + 1 q-1
. Now, consider the following quantity

E c α ∇f (x 1 (s)) + b(x 3 (s), α) q1 1 q 1 ≤ ≤ E c α ∇f (x 1 (s)) -c α ∇f (x 3 (s)) + c α ∇f (x 3 (s)) + b(x 3 (s), α) q1 1 q 1 ≤ E M x 1 (s) -x 3 (s) γ + L q1 1 q 1 ≤ E M x 1 (s) γ + M x 3 (s) γ + L q1 1 q 1 ≤ E M q1 x 1 (s) γq1 1 q 1 + E M q1 x 3 (s) γq1 1 q 1 + L,
where we have used assumption A2, assumption A6, Lemma 16 and Minkowski's inequality. By Lemma 9 and Lemma 10, we have

E c α ∇f (x 1 (s)) + b(x 3 (s), α) q1 1 q 1 ≤ ≤M c s d β 1/α + 1 + 1 γ + M c b s d β 1/α + 1 + 1 γ + L.
By combining the above inequalities, we get

E t 0 q x 1 (s) -x 3 (s)) q-2 x 1 (s) -x 3 (s), b 1 (x 1 (s-), α) -b(x 3 (s-), α) ds ≤ ≤q t 0 c s d β 1/α + 1 + 1 q-1 + c b s d β 1/α + 1 + 1 q-1 M c s d β 1/α + + 1 + 1 γ + M c b s d β 1/α + 1 + 1 γ + L ds =q t 0 M (c q-1 + c q-1 b )(c γ + c γ b ) s d β 1/α + 1 + 1 q-1+γ + L(c q-1 + c q-1 b ) s d β 1/α + + 1 + 1 q-1 ds ≤qt M (c q-1 + c q-1 b )(c γ + c γ b ) t d β 1/α + 1 + 1 q-1+γ + L(c q-1 + c q-1 b ) t d β 1/α + 1 + + 1 q-1 .
The final conclusion follows from this inequality.

Proof of Corollary 3

Proof. First, we replace t by kη. Then, by following the same lines of the proof of Corollary 1, we get W q q (µ 1t , µ 3t ) ≤C(k q+γ η + k q+γ η q β -q-1 α d q-1+γ ).

By assumption A4, q-1 < 1/p 1 and γ < 1/q 1 . It implies that d q-1+γ < d 1/p1+1/q1 = d. Hence, we have W q q (µ 1t , µ 3t ) ≤C(k q+γ η + k q+γ η q β -q-1 α d).

Proof of Corollary 4

Proof. By Lemma 2, Lemma 9 and Lemma 10, we have

c α |E[f (x 1 (t))] -E[f (x 3 (t))]| ≤ ≤ M E x 1 (t) γp 1 p + M E x 3 (t) γp 1 p + B W q (µ 1t , µ 3t ) ≤ M c t d β 1/α + 1 + 1 γ + M c b t d β 1/α + 1 + 1 γ + B W q (µ 1t , µ 3t ).
Then by Theorem 9, we have

c α |E[f (x 1 (t))] -E[f (x 3 (t))]| ≤ M c t d β 1/α + 1 + 1 γ + M c b t d β 1/α + 1 + 1 γ + + B qt M (c q-1 + c q-1 b )(c γ + c γ b ) t d β 1/α + 1 + 1 q-1+γ + + L(c q-1 + c q-1 b ) t d β 1/α + 1 + 1 q-1 1 q .
Applying Lemma 16 twice, we get

c α |E[f (x 1 (t))] -E[f (x 3 (t))]| ≤ ≤ M (c γ + c γ b ) t γ d γ β γ/α + t γ + 1 + B (qt) 1/q M 1/q (c q-1 + c q-1 b ) 1/q (c γ + c γ b ) 1/q td β 1/α + t + 1 (q-1+γ)/q + L 1/q (c q-1 + c q-1 b ) 1/q td β 1/α + t + 1 (q-1)/q ≤ M (c γ + c γ b ) t γ d γ β γ/α + t γ + 1 + B (qt) 1/q M 1/q (c q-1 + c q-1 b ) 1/q (c γ + c γ b ) 1/q
(td) (q-1+γ)/q β (q-1+γ)/(qα) + t (q-1+γ)/q + 1 + L 1/q (c q-1 + c q-1 b ) 1/q (td) (q-1)/q β (q-1)/(qα) + t (q-1)/q + 1 . Now, by replacing t = kη we find that, among the terms containing β, the largest power of d, the largest power of k and the smallest power of η are γ + q-1+γ q , γ + γ+q q and γ + 1 q , respectively. For the smallest power of β -1 , we need to compare the following quantities: γ/α, (q -1 + γ)/(qα) and (q -1)/(qα).

It is obvious that (q -1 + γ)/(qα) > (q -1)/(qα). Next, from the relation γ < 1/p = (q -1)/q, we have γ/α < (q -1)/(qα). Thus, the smallest power of β -1 is γ/α. Hence, we have the following bound:

c α |E[f (x 1 (t))] -E[f (x 3 (t))]| ≤C k γ+ γ+q q η γ+ 1 q β -γ α d γ+ q-1+γ q + k γ+ γ+q q η 1 q
, for some constant C > 0. For the power of d, using that γ < 1/p, q -1 < 1/p 1 and γ < 1/q 1 we have

γ + q -1 + γ q ≤1/p + 1/p 1 + 1/q 1 q = 1/p + 1/q = 1.
Finally, we have

c α |E[f (x 1 (t))] -E[f (x 3 (t))]| ≤C k γ+ γ+q q η γ+ 1 q β -γ α d + k γ+ γ+q q η 1 q .
Proof of Lemma 4

Proof. By Lemma 2, we have

c α |E[f (x 3 (t))] -E[f ( ŵ)]| ≤ M E x 3 (t) γp 1 p + M E ŵ γp 1 p + B W q (µ 3t , π).
Assumption A7 says that E ŵ γp is bounded by a constant depending on b, m and β. In addition, by Proposition 1, lim t→∞ W γp (µ 3t , π) = 0, and by Theorem 7.12 in [START_REF] Villani | Topics in optimal transportation[END_REF], it follows that

lim t→∞ E x 3 (t) γp = E ŵ γp .
Thus, E x 3 (t) γp is bounded by a constant independent of t. Finally, since q < α, by propostion 1 again, W q (µ 3t , π) ≤ Cβe -λ * t/β . Hence, using the bound in assumption A7, there exists constant C such that

|E[f (x 3 (t))] -E[f ( ŵ)]| ≤ Cβ b + d/β m exp(-λ * β -1 t).

Proof of Lemma 5

Proof. The proof is adapted from [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF], Section 3.5. First, we have the decomposition:

E[f ( ŵ)] = R d f (w) exp(-βf (w)) R d exp(-βf (v))dv dw = 1 β - R d exp(-βf (w)) R d exp(-βf (v))dv log exp(-βf (w)) R d exp(-βf (v))dv dw -log R d exp(-βf (v))dv .
The first term in the parentheses is the differential entropy of the probability density of ŵ, which has a finite second moment (due to assumption A7). Hence, it is upperbounded by the differential entropy of a Gaussian density with the same second moment:

-

R d exp(-βf (w)) R d exp(-βf (v))dv log exp(-βf (w)) R d exp(-βf (v))dv dw ≤ d 2 log 2πe(b + d/β) dm .
Here, π denotes the Archimedes' constant. By Lemma 8, we have

-log R d exp(-βf (w))dw ≤ βf (w * ) + β -γ M c -1 α 1 + γ -log π d/2 β -d Γ(d/2 + 1) .
Then, it implies that

E[f ( ŵ)] ≤ dβ -1 2 log 2πe(b + d/β) dm + f (w * ) + β -γ-1 M c -1 α 1 + γ -β -1 log π d/2 β -d Γ(d/2 + 1) =f (w * ) + β -γ-1 M c -1 α 1 + γ + β -1 log (2e(b + d/β)) d/2 Γ(d/2 + 1)β d (dm) d/2
, which leads to desired result.

Proof of Corollary 6

Proof. By triangular inequality, we have W q (µ 2t , π) ≤ W q (µ 2t , µ 1t ) + W q (µ 1t , µ 3t ) + W q (µ 3t , π).

Then, using Corollary 1, Corollary 3 and Proposition 1, we get W q (µ 2t , π) ≤ ≤C (k 2 η + k 2 η 1+γ/α β -γ(q-1)/α d) 1/q + (k q+γ η + k q+γ η q β -(q-1)/α d) 1/q + βe -λ * kη/β ≤C k 2/q η 1/q + k 2/q η 1/q+γ/(qα) β -γ(q-1)/(qα) d 1/q + k 1+γ/q η 1/q + k 1+γ/q ηβ -(q-1)/(qα) d 1/q + βe -λ * kη/β , where, we used Lemma 16 for the second inequality. Then, similar to the proof of Corollary 1, we obtain W q (µ 2t ,π) ≤ ≤C k max{2,q+γ}/q η 1/q + k max{2,q+γ}/q η 1/q+γ/(qα) β -γ(q-1)/(qα) d 1/q + βe -λ * kη/β .

Proof of Theorem 10

Proof. Since each function x → f (i) (x) satisfies assumptions A1-7, it is easy to check that f k also satisfies these assumptions (with the same constants and the same parameters) for all k. Then by repeating exactly the same lines as in the proof of Lemma 12, we obtain the same estimates for the moments of x 2 . Now by following the same steps as in the proof of Theorem 13, we first have

W q q (µ 1t , µ 2t ) ≤ ≤q k-1 j=0 (j+1)η jη E x 1 (s) -x 2 (s) (q-1)p1 1 p 1 E c α (∇f (x 1 (s)) -∇f k (x 2 (jη))) q1 1 q 1 ds, then E x 1 (s) -x 2 (s) (q-1)p1 1 p 1 ≤ ≤c 1 p 1 + (E x 2 (0) (q-1)p1 ) 1 p 1 + η q-1 M q-1 (E x 2 (0) (q-1)p1γ ) 1 p 1 + B (q-1) + η β q-1 α l 1 p 1
α,(q-1)p1,d + (j + 1)

1 p 1 cη d β 1/α 1 p 1 + (cη) 1 p 1 + (2η(b + m)) (q-1) 2 + 2 (q-1) 2 (ηB) (q-1) + η β (q-1) α l 1 p 1 α,(q-1)p1,d + η q-1 M q-1 (2η(b + m)) (q-1)γ 2 + 2 (q-1)γ 2 (ηB) (q-1)γ + η β (q-1)γ α l 1 p 1 α,(q-1)p1γ,d =Q 1 (η) + (j + 1) 1 p 1 P 1 (η), 91
where P 1 (η) and Q 1 (η) are defined in Theorem 13. Now, by Minkowski's inequality, we have

E c α (∇f (x 1 (s))-∇f k (x 2 (jη))) q1 1 q 1 = = E c α (∇f (x 1 (s)) -∇f (x 2 (jη)) + ∇f (x 2 (jη)) -∇f k (x 2 (jη))) q1 1 q 1 ≤ E c α (∇f (x 1 (s)) -∇f (x 2 (jη))) q1 1 q 1 + E c α (∇f (x 2 (jη)) -∇f k (x 2 (jη))) q1 1 q 1 .
As in the proof of Theorem 13, the following inequality holds:

E c α ∇f (x 1 (s))-c α ∇f (x 2 (jη)) q1 1 q 1 ≤ ≤M (E x 2 (0) γq1 ) 1 q 1 + M c 1 q 1 + M (j + 1) 1 q 1 cη d β 1/α 1 q 1 + (cη) 1 q 1 + (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 q 1 α,γq1,d =Q 2 + (j + 1) 1 q 1 P 2 (η),
where P 2 (η) and Q 2 are defined in theorem 13. Using the additional assumption, lemma 12, and lemma 16, we get

E c α (∇f (x 2 (jη))-∇f k (x 2 (jη))) q1 1 q 1 ≤ ≤δ E M q1 x 2 (jη) γq1 1 q 1 ≤δ M q1 (E x 2 (0) γq1 ) + M q1 j (2η(b + m)) γq 1 2 + 2 γq 1 2 (ηB) γq1 + η β γq 1 α l α,γq1,d 1 q 1 ≤δM (E x 2 (0) γq1 ) 1 q 1 + δM (j + 1) 1 q 1 (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 q 1 α,γq1,d .
By combining the two above inequalities, we obtain

E c α ∇f (x 1 (s))-c α ∇f (x 2 (jη)) q1 1 q 1 ≤ ≤(1 + δ)M (E x 2 (0) γq1 ) 1 q 1 + M c 1 q 1 + M (j + 1) 1 q 1 cη d β 1/α 1 q 1 + (cη) 1 q 1 + (1 + δ)(2η(b + m)) γ 2 + (1 + δ)2 γ 2 (ηB) γ + (1 + δ) η β γ α l 1 q 1 α,γq1,d =Q 2 + (j + 1) 1 q 1 P 2 (η).
92 Finally, we have

W q q (µ 1t , µ 2t ) ≤ qη k 2 P 1 (η)P 2 (η) + k 1+1/p1 P 1 (η)Q 2 + k 1+1/q1 P 2 (η)Q 1 (η) + kQ 1 (η)Q 2 .
By considering the additional term δ, we arrive at the following bound:

W q q (µ 1t , µ 2t ) ≤ C(1 + δ)(k 2 η + k 2 η 1+γ/α β -γ(q-1)/α d).
Proof of corollary 7

Proof. By lemma 2,

c α E[f (x 1 (kη))]-E[f (x 2 (kη))] ≤ ≤ M E x 1 (kη) γp 1 p + M E x 2 (kη) γp 1 p + B W q (µ 1t , µ 2t ).
Then, by following the same proof as in corollary 8, corollary 1 and using theorem 10, we get

c α E[f (x 1 (kη))] -E[f (x 2 (kη))] ≤C(1 + δ) k 1+ 1 q η 1 q + k 1+ 1 q η 1 q + γ αq β -(q-1)γ αq d .

Technical results

Corollary 8. Along with P 1 (η), P 2 (η), Q 1 (η), Q 2 in Lemma 13, we define, in addition, the following quantities:

P 3 (η) M cη d β 1/α 1 p + (cη) 1 p + (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 p α,γp,d Q 3 M (E X 2 (0) γp ) 1 p + M c 1 p + B.
For 0 < η < m M 2 , we have the following bound:

c α E[f (x 1 (kη))] -E[f (x 2 (kη))] ≤ ≤(qη) 1 q k 1+ 1 q (P 1 (η)P 2 (η)) 1 q P 3 (η) + k 1+ 1 qp 1 (P 1 (η)Q 2 ) 1 q P 3 (η) + k 1+ 1 qq 1 (P 2 (η)Q 1 (η)) 1 q P 3 (η) + k(Q 1 (η)Q 2 ) 1 q P 3 (η) + k 2 q (P 1 (η)P 2 (η)) 1 q Q 3 + k 1 q + 1 qp 1 (P 1 (η)Q 2 ) 1 q Q 3 + k 1 q + 1 qq 1 (P 2 (η)Q 1 (η)) 1 q Q 3 + k 1 q (Q 1 (η)Q 2 ) 1 q Q 3 . Proof. By Lemma 2, c α E[f (x 1 (kη))] -E[f (x 2 (kη))] ≤ ≤ M E x 1 (kη) γp 1 p + M E x 2 (kη) γp 1 p + B W q (µ 1t , µ 2t ).
Using Lemma 9 and Lemma 13, we have

M E x 1 (kη) γp 1 p +M E x 2 (kη) γp 1 p + B ≤ ≤M c kη d β 1/α + 1 + 1 γ + M (E x 2 (0) γp ) + k (2η(b + m)) γp 2 + 2 γp 2 (ηB) γp + η β γp α l α,γp,d 1 p + B.
By using Lemma 16, we obtain

M E x 1 (kη) γp 1 p +M E x 2 (kη) γp 1 p + B ≤ ≤M c kη d β 1/α + 1 + 1 γ + M (E x 2 (0) γp ) 1 p + M k 1 p (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 p α,γp,d + B.
We note that γ < 1 p . Hence,

c kη d β 1/α + 1 + 1 γ ≤ c kη d β 1/α + 1 + 1 1 p ≤k 1 p cη d β 1/α + 1 1 p + c 1 p ,
where the last inequality is an application of Lemma 16. By replacing this inequality into the previous one and rearranging the terms, we have

M E x 1 (kη) γp 1 p + M E x 2 (kη) γp 1 p + B ≤ ≤M (E x 2 (0) γp ) 1 p + M c 1 p + B + M k 1 p cη d β 1/α + 1 1 p + (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 p α,γp,d ≤M (E x 2 (0) γp ) 1 p + M c 1 p + B + M k 1 p cη d β 1/α 1 p + (cη) 1 p + (2η(b + m)) γ 2 + 2 γ 2 (ηB) γ + η β γ α l 1 p α,γp,d =Q 3 + k 1 p P 3 (η).
Here, we have used lemma 16 in the last inequality. Next, by Lemma 13 and Lemma 16,

W q (µ 1t , µ 2t ) ≤(qη) 1 q k 2 P 1 (η)P 2 (η) + k 1+1/p1 P 1 (η)Q 2 + k 1+1/q1 P 2 (η)Q 1 (η) + kQ 1 (η)Q 2 1 q ≤(qη) 1 q k 2 q (P 1 (η)P 2 (η)) 1 q + k 1 q + 1 qp 1 (P 1 (η)Q 2 ) 1 q + k 1 q + 1 qq 1 (P 2 (η)Q 1 (η)) 1 q + + k 1 q (Q 1 (η)Q 2 ) 1 q .
Proof. We start by writing: log

R d exp(-βf (w))dw = -βf (w * ) + log R d exp -β(f (w) -f (w * )) dw ≥ -βf (w * ) + log R d exp - βM c -1 α 1 + γ w -w * 1+γ dw.
Here, we used Lemma 7, with ∇f (w * ) = 0. For the second term on the right hand side, we have

R d exp - βM c -1 α 1 + γ w -w * 1+γ dw = w ≤β -1 exp - βM c -1 α 1 + γ w 1+γ dw + w ≥β -1 exp - βM c -1 α 1 + γ w 1+γ dw ≥ w ≤β -1 exp - βM c -1 α 1 + γ β -1-γ dw + 0 = exp - β -γ M c -1 α 1 + γ w ≤β -1 1dw = exp - β -γ M c -1 α 1 + γ π d/2 β -d Γ(d/2 + 1) ,
where, Γ denotes the Gamma function and π denotes Archimedes' constant (here, it is not the invariant distribution). Hence, log

R d exp - βM c -1 α 1 + γ w -w * 1+γ dw ≥ - β -γ M c -1 α 1 + γ + log π d/2 β -d Γ(d/2 + 1)
.

By combining the above inequalities, we have the desired result.

Lemma 9. For λ ∈ (0, 1), there exists a constant c depending on m, b, α, such that

E x 1 (t) λ 1 λ ≤ c t(dβ -1/α + 1) + 1 , ∀t > 0, β ≥ 1, 1 < α < 2.
Proof. We follow exactly the same proof as Lemma 7.1 in Xie and Zhang [2017], with some modifications. Let h(x) (1 + x 2 ) 1/2 . By Itô's formula, we have (7.11) where M (t) is a local martingale. Noticing that ∂ i h(x) = x i (1 + x 2 ) -1/2 /2 and using assumption A3, we have

dh(x 1 (t)) = b 1 (x 1 (t)), ∇h(x 1 (t)) + R d h(x 1 (t) + β -1/α x) -h(x 1 (t)) -I x <1 β -1/α x, ∇h(x 1 (t)) ν(dx) dt + dM (t),
b 1 (x), ∇h(x) = b 1 (x), x (1 + x 2 ) -1/2 /2 ≤(-m x 1+γ + b)(1 + x 2 ) -1/2 /2 =(-m( x 1+γ + 1) + m + b)(1 + x 2 ) -1/2 /2. Since ( x 2 + 1) (1+γ)/2 ≤ ( x 1+γ + 1) by Lemma 16, it follows that b 1 (x), ∇h(x) ≤(-m( x 2 + 1) (1+γ)/2 + m + b)(1 + x 2 ) -1/2 /2 =(-m( x 2 + 1) γ/2 + (m + b)(1 + x 2 ) -1/2 )/2 ≤(-m( x 2 + 1) γ/2 + m + b)/2 =(-mh(x) γ + m + b)/2.
On the other hand, observing that

|h(x + y) -h(x)| ≤ y 1 0 ∇h(x + sy) ds ≤ y /2, and
h(x + y) -h(x) -y, ∇h(x) ≤ y 2 /2, we have R d h(x 1 (t) + x) -h(x 1 (t))-I x <1 x, ∇h(x 1 (t)) ν(dx) ≤ ≤ 1 2β 2/α x <1 x 2 ν(dx) + 1 2β 1/α x ≥1 x ν(dx) ≤C d β 1/α ,
where the last inequality is due to Lemma 15. By integrating (7.11) and combining the above inequalities, we have

h(x 1 (t)) -h(x 1 (0)) ≤ t 0 (-mh(x 1 (s)) γ + m + b)/2 + C d β 1/α ds + M (t) ≤ t 0 (m + b)/2 + C d β 1/α ds + M (t).
By Lemma 3.8 in Xie and Zhang [2017], for λ ∈ (0, 1),

E sup s∈[0,t] h(x 1 (s)) λ ≤ c λ Eh(x 1 (0)) + ((m + b)/2 + C d β 1/α )t λ .
This leads to the conclusion since h(x) ≥ x .

Lemma 10. For λ ∈ (0, 1), there exists a constant c b depending on L, m, b, α, such that

E x 3 (t) λ 1 λ ≤ c b t(dβ -1/α + 1) + 1 , ∀t > 0, β ≥ 1, 1 < α < 2.
Proof. The proof is similar to the proof of Lemma 9.

Lemma 11. Let X be a scalar symmetric α-stable distribution with α < 2, i.e. X ∼ SαS(1) (see Definition 3), then, for -1 < λ < α,

E(|X| λ ) = 2 λ Γ((1 + λ)/2)Γ(1 -λ/α) Γ(1/2)Γ(1 -λ/2) .
Proof. The proof follows from Theorem 3 in [START_REF] Shanbhag | On certain self-decomposable distributions[END_REF] (see also equation ( 13) in [START_REF] Matsui | Fractional absolute moments of heavy tailed distributions[END_REF]).

Corollary 9. The quantity l α,λ,d E L α (1) λ is finite for 0 ≤ λ < α. For details, we have

(a) If 1 < λ < α, then E L α (1) λ ≤ d λ 2 λ Γ((1 + λ)/2)Γ(1 -λ/α) Γ(1/2)Γ(1 -λ/2) . (b) If 0 ≤ λ ≤ 1, then E L α (1) λ ≤ d 2 λ Γ((1 + λ)/2)Γ(1 -λ/α) Γ(1/2)Γ(1 -λ/2) .
Proof. Since L α (1), by definition, is a d-dimensional vector whose components are i.i.d symmetric α-stable distributions L α i (1) for i ∈ {1, . . . , d}, we have

L α (1) ≤ d i=1 |L α i (1)| (a) 1 < λ < α.
By using Minkowski's inequality and Lemma 11,

(E L α (1) λ ) 1/λ ≤ E d i=1 |L α i (1)| λ 1/λ ≤ d i=1 (E|L α i (1)| λ ) 1/λ =d 2 λ Γ((1 + λ)/2)Γ(1 -λ/α) Γ(1/2)Γ(1 -λ/2) 1/λ
. Thus, we have

E L α (1) λ ≤ d λ 2 λ Γ((1 + λ)/2)Γ(1 -λ/α) Γ(1/2)Γ(1 -λ/2) .
(b) 0 ≤ λ ≤ 1. By using Lemma 16 and Lemma 11 ,

E L α (1) λ ≤E d i=1 |L α i (1)| λ ≤ d i=1 E|L α i (1)| λ =d 2 λ Γ((1 + λ)/2)Γ(1 -λ/α) Γ(1/2)Γ(1 -λ/2) .
Lemma 12. Let us denote the value E L α (1) λ by l α,λ,d < ∞. For 0 < η ≤ m M 2 and s ∈ [jη, (j + 1)η), we have the following estimates:

(a) If 1 < λ < α and 1 < γλ < α then E x 2 (jη) λ ≤ B j,λ E x 2 (0) λ 1 λ + j (2η(b + m)) 1 2 + 2 1 2 ηB + η β 1 α l 1 λ α,λ,d λ , E x 2 (s) λ ≤ B 1 λ j,λ + (s -jη) M B 1 λ j,γλ + B + s -jη β 1 α l 1 λ α,λ,d λ . (b) If 0 ≤ λ ≤ 1 then E x 2 (jη) λ ≤ Bj,λ E x 2 (0) λ + j (2η(b + m)) λ 2 + 2 λ 2 (ηB) λ + η β λ α l α,λ,d , E x 2 (s) λ ≤ Bj,λ + (s -jη) λ M λ Bj,γλ + B λ + s -jη β λ α l α,λ,d . (c) If 1 < λ < α and 0 ≤ γλ ≤ 1 then E x 2 (jη) λ ≤ B j,λ , E x 2 (s) λ ≤ B 1 λ j,λ + (s -jη) M B 1 λ j,γλ + B + s -jη β 1 α l 1 λ α,λ,d λ . Proof. Starting from x 2 ((j + 1)η) = x 2 (jη) -ηc α ∇f (x 2 (jη)) + η β 1 α L α (1),
we have either (by Minkowski, if λ > 1) (7.12) or (by Lemma 16, if 0 ≤ λ ≤ 1)

E x 2 ((j + 1)η) λ 1 λ ≤ E x 2 (jη) -ηc α ∇f (x 2 (jη)) λ 1 λ + η β 1 α E L α (1) λ 1 λ ,
E x 2 ((j + 1)η) λ ≤ E x 2 (jη) -ηc α ∇f (x 2 (jη)) λ + η β λ α E L α (1) λ . (7.13)
We have

x 2 (jη) -ηc α ∇f (x 2 (jη)) λ = = x 2 (jη) -ηc α ∇f (x 2 (jη)) 2× λ 2 (7.14) = x 2 (jη) 2 -2ηc α x 2 (jη), ∇f (x 2 (jη) + η 2 c α ∇f (x 2 (jη) 2 λ 2 ≤ x 2 (jη) 2 -2η(m x 2 (jη) 1+γ -b) + η 2 (2M 2 x 2 (jη) 2γ + 2B 2 ) λ 2 , (7.15)
where we used assumption A3 and Lemma 13. For 0

< η ≤ m M 2 , 2ηm( X 2 (jη) 1+γ + 1) ≥ 2η 2 M 2 X 2 (jη) 2γ . (since 1 + γ > 2γ and ηm > η 2 M 2 )
Using this inequality we have

x 2 (jη) -ηc α ∇f (x 2 (jη)) λ ≤ x 2 (jη) 2 + 2η(b + m) + 2η 2 B 2 λ 2 ≤ x 2 (jη) λ + (2η(b + m)) λ 2 + 2 λ 2 (ηB) λ . (by Lemma 16) (7.16)
Consider the case where λ > 1. By (7.12) and (7.16),

E x 2 ((j + 1)η) λ 1 λ ≤ ≤ E x 2 (jη) λ + (2η(b + m)) λ 2 + 2 λ 2 (ηB) λ 1 λ + η β 1 α E L α (1) λ 1 λ ≤ E x 2 (jη) λ 1 λ + (2η(b + m)) 1 2 + 2 1 2 ηB + η β 1 α l 1 λ α,λ,d (by Lemma 16) ≤ E x 2 (0) λ 1 λ + (j + 1) (2η(b + m)) 1 2 + 2 1 2 ηB + η β 1 α l 1 λ α,λ,d .
For the case where 0 ≤ λ ≤ 1, by (7.13) and (7.16),

E x 2 ((j + 1)η) λ ≤ E x 2 (jη) λ + (2η(b + m)) λ 2 + 2 λ 2 (ηB) λ + η β λ α l α,λ,d ≤ E x 2 (0) λ + (j + 1) (2η(b + m)) λ 2 + 2 λ 2 (ηB) λ + η β λ α l α,λ,d .
Now, from the identification, for s ∈ [jη, (j + 1)η),

x 2 (s) = x 2 (jη) + (s -jη)c α ∇f (x 2 (jη)) + s -jη β 1 α L α (1), we have x 2 (s) ≤ x 2 (jη) + (s -jη)c α ∇f (x 2 (jη)) + s -jη β 1 α L α (1) ≤ x 2 (jη) + (s -jη)(M x 2 (jη) γ + B) + s -jη β 1 α L α (1) . 100 For λ > 1, E x 2 (s) λ 1 λ ≤ E x 2 (jη) λ 1 λ + (s -jη) M E x 2 (jη) γλ 1 λ + B + s -jη β 1 α l 1 λ α,λ,d . For λ ≤ 1, E x 2 (s) λ ≤ E x 2 (jη) λ + (s -jη) λ M λ E x 2 (jη) γλ + B λ + s -jη β λ α l α,λ,d .
By replacing the estimate of E x 2 (jη) λ , we obtain the desired result.

Lemma 13. Under assumptions A1 and A2 we have

c α ∇f (w) ≤ M w γ + B, ∀w ∈ R d .
Proof. By assumption A2 we have

c α ∇f (w) -∇f (0) ≤ M w -0 γ .
Since c α ∇f (0) ≤ B by assumption A1, the conclusion follows. Here, we have used the inequality w ≤ w 1+γ + 1.

Lemma 15. Let ν be the Lévy measure of a d-dimensional Lévy process L α whose components are independent scalar symmetric α-stable Lévy processes L α 1 , . . . , L α d . Then there exists a constant C > 0 such that the following inequality holds with β ≥ 1 and

2 > α > 1: 1 β 2/α x <1 x 2 ν(dx) + 1 β 1/α x ≥1 x ν(dx) ≤ C d β 1/α .
Proof. Using Lemma 4.1 in [START_REF] Kallsen | Characterization of dependence of multidimensional Lévy processes using Lévy copulas[END_REF], we have

x <1 x 2 ν(dx) = d i=1 |xi|<1 |x i | 2 1 |x i | 1+α dx i = d i=1 2 2 -α = 2d 2 -α .
Similarly, we have

x ≥1

x ν(dx)

= d i=1 |xi|≥1 |x i | 1 |x i | 1+α dx i = d i=1 2 α -1 = 2d α -1 .
Combining these two equalities, we have the desired conclusion.

Lemma 16. For a, b ≥ 0 and 0 ≤ γ ≤ 1, we have the following inequality:

(a + b) γ ≤ a γ + b γ .
Proof. If a = b = 0, the inequality is trivial. Hence, let us assume that a > b ≥ 0. We have

1 + b a γ ≤ 1 + γ b a (by Bernoulli's inequality) ≤ 1 + b a (since 0 ≤ γ ≤ 1 and b a ≥ 0) ≤ 1 + b a γ . (since 0 ≤ γ ≤ 1 and 0 ≤ b a ≤ 1)
By multiplying both sides by a γ > 0, we have the conclusion.

Proof. We have for t ∈ [kη, (k + 1)η],

w t -w kη ≤ t kη ∇f (w s ) ds + εσ B(t) -B(kη) + ε L α (t) -L α (kη) ≤ t kη ∇f (w s ) -∇f (w kη ) ds + η ∇f (w kη ) + εσ B(t) -B(kη) + ε L α (t) -L α (kη) ≤ t kη M w s -w kη γ ds + η(M w kη γ + B) + εσ B(t) -B(kη) + ε L α (t) -L α (kη) .
For γ < 1, using that w s -w kη γ ≤ w s -w kη + 1, we get:

w t -w kη ≤ t kη M w s -w kη ds + η(M w kη γ + B + M ) + εσ sup t∈[kη,(k+1)η] B(t) -B(kη) + ε sup t∈[kη,(k+1)η] L α (t) -L α (kη) .
Then the Gronwall Lemma gives: L α (t) -L α (kη) .

sup t∈[kη,(k+1)η] w t -w kη ≤e M η η(M w kη γ + B + M ) + εσ sup t∈[kη,(k+1)η] B(t) -B(kη) + ε sup t∈[kη,(k+1)η] L α (t) -L α (kη) .
By Lemma 7.1 in Xie and Zhang [2017], Lemma 9 and Markov's inequality, for any u > 0, we have:

P[ max 0≤k≤K-1 w kη γ ≥ u] ≤ E[max 0≤k≤K-1 w kη γ ] u ≤ C 1 (Kη(dε + 1) + 1) γ u ,
where C 1 is a constant independent of K, η, ε and d. By Lemma 19, we have: L α (t) -L α (kη) ≥ ε -1 ξ/4] ≤ C 1 (Kη(dε + 1) + 1) γ e M η M η ξ + 1 -1 -Cde -ξ 2 e -2M η (εσ) -2 /(16dη) K + 1 -1 -C α d 1+α/2 ηe αM η ε α ξ -α K .

Now we prove the following lemma. L α (t) -L α (kη) ≥ u] ≤ 1 -1 -C α d 1+α/2 ηu -α K .

Proof. We have =1 -1 -Cde -u 2 /(dη) K .

Similarly, we have L α (t) -L α (kη) ≥ u] ≤1 -1 -C α d 1+α/2 ηu -α K .

w satisfies dw t = b(w)dt+εσdB φ (t)+εdL α (t). Since this equation has a unique solution (constructed explicitly with the Euler scheme), we conclude that w has the same law under Q as ŵ under P.

We thus have: KL(μ t , µ t ) = KL(P t , Q t ) = E P log dP dQ Ft = 1 2ε 2 σ 2 E P t 0 b( ŵ) + ∇f ( ŵs ) 2 ds (7.18) By using the same steps of the proof of Raginsky et al. [2017][Lemma 3.6

+ (εη 1/2 ) 2γ d 2γ 2 γ Γ 2γ+1 2 √ π ≤K 1 k 2 η 1+2γ+γ 2 + K 2 kη 1+2γ + K 3 kη 1+ 2γ α + K 4 kη 1+γ ,
where

K 1 CM 2+2γ 3 γ ε 2 σ 2 max (2(b + m)) γ 2 , 2 γ 2 B 2γ 2 , ε 2γ 2 d 2 2γ 2 Γ((1 + 2γ 2 )/2)Γ(1 -2γ 2 /α) Γ(1/2)Γ(1 -γ 2 ) , ε 2γ 2 d 2 γ 2 Γ 2γ 2 +1 2 √ π , K 2 CM 2+2γ 3 γ 2ε 2 σ 2 E ŵ0 2γ 2 + B 2 M 2 , K 3 M 2 3 γ ε 2γ-2 d 2γ 2σ 2 2 2γ Γ((1 + 2γ)/2)Γ(1 -2γ/α) Γ(1/2)Γ(1 -γ) , K 4 M 2 3 γ ε 2γ-2 d 2γ 2σ 2 2 γ Γ 2γ+1 2 √ π .
Proof. Let us consider the term ŵsŵjη , for s ∈ [jη, (j + 1)η]:

ŵsŵjη = -(s -jη)∇f ( ŵjη ) + ε(L s -L jη ) + ε(B s -B jη ) (7.21)

T 1 + T 2 + T 3 (7.22)
Using this equation and (7.20), we obtain:

KL(μ t , µ t ) ≤ M 2 2ε 2 σ 2 k-1 j=0 (
≤ M 2 3 γ 2ε 2 σ 2 k-1 j=0 (j+1)η jη CM 2γ η 2γ E ŵ0 2γ 2 + j (2η(b + m)) γ 2 + 2 γ 2 (ηB) 2γ 2 + ε 2γ 2 η 2γ 2 α d 2 2γ 2 Γ((1 + 2γ 2 )/2)Γ(1 -2γ 2 /α) Γ(1/2)Γ(1 -γ 2 ) + ε 2γ 2 η γ 2 d 2 γ 2 Γ 2γ 2 +1 2 √ π + B 2 M 2 + (εη 1/α ) 2γ d 2γ 2 2γ Γ((1 + 2γ)/2)Γ(1 -2γ/α) Γ(1/2)Γ(1 -γ) + (εη 1/2 ) 2γ d 2γ 2 γ Γ 2γ+1 2 √ π ds = M 2 3 γ 2ε 2 σ 2 kη CM 2γ η 2γ E ŵ0 2γ 2 + k -1 2 (2η(b + m)) γ 2 + 2 γ 2 (ηB) 2γ 2 + ε 2γ 2 η 2γ 2 α d 2 2γ 2 Γ((1 + 2γ 2 )/2)Γ(1 -2γ 2 /α) Γ(1/2)Γ(1 -γ 2 ) + ε 2γ 2 η γ 2 d 2 γ 2 Γ 2γ 2 +1 2 √ π + B 2 M 2 + (εη 1/α ) 2γ d 2γ 2 2γ Γ((1 + 2γ)/2)Γ(1 -2γ/α) Γ(1/2)Γ(1 -γ) + (εη 1/2 ) 2γ d 2γ 2 γ Γ 2γ+1 2 √ π .
By defining the constants K 1 , K 2 , K 3 and K 4 as in the statement of the theorem, we directly have the conclusion.

The proof of Theorem 12 is given bellow.

Proof. By Theorem 15, we have KL(μ t , µ t ) ≤K 1 k 2 η 1+2γ+γ 2 + K 2 kη 1+2γ + K 3 kη 1+ 2γ α + K 4 kη 1+γ .

We can easily check that, for example, if 0 < η ≤ δ 2 2K1t 2 1 γ 2 +2γ-1 , then K 1 k 2 η 1+2γ+γ 2 ≤ δ 2 2 . By the same arguments, we finally have

KL(μ t , µ t ) ≤ δ 2 2 + δ 2 2 + δ 2 2 + δ 2 2 =2δ 2 .
This finalizes the proof.

Technical results

Lemma 21. Under assumptions A10 and A11 we have ∇f (w) ≤ M w γ + B, ∀w ∈ R d .

Proof. By assumption A10 we have ∇f (w) -∇f (0) ≤ M w -0 γ .

Since ∇f (0) ≤ B by assumption A11, the conclusion follows.

For the moments of the noise B(1), we have the following corollary.

Lemma 22. Let X be a scalar standard Gaussian random variable. Then, for λ > -1, we have

E(|X| λ ) = 2 λ/2 Γ λ+1 2 √ π ,
where Γ denotes the Gamma function.

Proof. The result is a direct consequence of equation ( 17 Details of the simulations in Section 6.4

The detailed settings of the parameters for the synthetic experiment (Figure 6.1) are as follows. In this thesis, we are concerned with the Stochastic Gradient Descent (SGD) algorithm. Specifically, we perform theoretical and empirical analysis of the behavior of the stochastic gradient noise (GN), which is defined as the difference between the true gradient and the stochastic gradient, in deep neural networks. Based on these results, we bring an alternative perspective to the existing approaches for investigating SGD. The GN in SGD is often considered to be Gaussian for mathematical convenience. This assumption enables SGD to be studied as a stochastic differential equation (SDE) driven by a Brownian motion. We argue that the Gaussianity assumption might fail to hold in deep learning settings and hence render the Brownian motion-based analyses inappropriate. Inspired by non-Gaussian natural phenomena, we consider the GN in a more general context that suggests that the GN is better approximated by a heavy-tailed α-stable random vector. Accordingly, we propose to analyze SGD as a discretization of an SDE driven by a L évy motion. Firstly, to justify the α-stable assumption, we conduct experiments on common deep learning scenarios and show that in all settings, the GN is highly non-Gaussian and exhibits heavy-tails. Secondly, under the heavy-tailed GN assumption, we provide a non-asymptotic analysis for the discrete-time dynamics SGD to converge to the global minimum in terms of suboptimality. Finally, we investigate the metastability nature of the SDE driven by L évy motion that can then be exploited for clarifying the behavior of SGD, especially in terms of 'preferring wide minima'. More precisely, we provide formal theoretical analysis where we derive explicit conditions for the step-size such that the metastability behavior of SGD, viewed as a discrete-time SDE, is similar to its continuous-time limit. We show that the behaviors of the two systems are indeed similar for small stepsizes and we describe how the error depends on the algorithm and problem parameters. We illustrate our metastability results with simulations on a synthetic model and neural networks. Our results open up a different perspective and shed more light on the view that SGD prefers wide minima.

Institut Polytechnique de Paris 91120 Palaiseau, France

Finally

  , I would like to thank the French National Research Agency (ANR) and the industrial chair Data Science & Artificial Intelligence from Télécom Paris. This thesis is supported by the ANR as a part of the FBIMATRIX (ANR-16-CE23-0014) project, and by the industrial chair Data Science & Artificial Intelligence from Télécom Paris.

  Gürbüzbalaban, M., Nguyen, T. H., Richard, G., & Sagun, L. On the Heavy-Tailed Theory of Stochastic Gradient Descent for Deep Neural Networks. Submitted to Journal of Machine Learning Research, 2019. Under revision. This work is mainly presented in Chapter 4 of this document. • Nguyen, T. H., S ¸imşekli, U., & Richard, G. Non-Asymptotic Analysis of Fractional Langevin Monte Carlo for Non-Convex Optimization. International Conference on Machine Learning, 2019. This work is mainly presented in Chapter 5 of this document. • Nguyen, T. H., S ¸imşekli, U., Gürbüzbalaban, M., & Richard, G. First Exit Time Analysis of Stochastic Gradient Descent Under Heavy-Tailed Gradient Noise. Neural Information Processing Systems, 2019.

•

  Nguyen, T. H., S ¸imşekli, U., Richard, G., & Cemgil, A. T. Efficient Bayesian Model Selection in PARAFAC via Stochastic Thermodynamic Integration. IEEE Signal Processing Letters, 2018. In this work, we develop a novel parallel and distributed Bayesian model selection technique for rank estimation in a large-scale tensor factorization model, called the Parallel factor analysis (PARAFAC). The proposed approach integrates ideas from stochastic gradient MCMC, statistical physics, and distributed stochastic optimization. Our method has a clear mathematical interpretation, and has significantly lower computational requirements, thanks to data sub-sampling and parallelization. • S ¸imşekli, U., Yıldız, C., Nguyen, T. H., Richard, G., & Cemgil, A. T. Asynchronous Stochastic Quasi-Newton MCMC for Non-Convex Optimization. International Conference on Machine Learning, 2018.

Figure 2 . 1 :

 21 Figure 2.1: Illustration of SαS (a), L α t (b). As α gets smaller, SαS becomes heaviertailed and consequently, L α (t) incurs larger jumps.

(

  SDE) driven by a Brownian motion. However, the Gaussianity assumption might fail to hold in deep learning settings and hence render the Brownian motion-based analyses inappropriate. Inspired by non-Gaussian natural phenomena, we consider a better-suited hypothesis for the stochastic gradient noise that has more pertinent theoretical implications for the phenomena associated with SGD. Accordingly, we propose to analyze SGD as a discretization of an SDE driven by a Lévy motion.This chapter is based on the article[S ¸imşekli et al., 2019].

Figure 4

 4 Figure 4.1: A function with a narrow minimum of width m 1 and a wide minimum of width m 2 .

Figure 4

 4 Figure 4.2: A mismatch between the Gaussianity assumption and the empirical behavior of the stochastic gradient noise. (a) and (b) The histogram of the norms of the gradient noises after the first and the last iterations, respectively (computed with AlexNet on CIFAR10). (c) and (d) the histograms of the norms of (scaled) Gaussian and α-stable random vectors, respectively.

Figure 4

 4 Figure 4.3: Illustration of the tail-index estimator α.

  )) are relatively close to each other, i.e. if the 'condition number' c st max{|α X -α 12 |, |α X -α 123 |} is smaller than some threshold.

Figure 4

 4 Figure 4.4: Stability confidence results for each layer as well as for whole network (indicated by layer index at 0). From left to right: depth 3, depth 5, depth 7.

  Figure 4.4(a), in which layer index at 0 corresponds to the whole network while the indices 1, 2, . . . , 7 represent the hidden layers of the network. Our experiments show that the condition number c st for the whole network are always smaller than the threshold 0.05, which means the gradient noise of the network satisfies our required stability criterion, even when we change the number of layers (depths) and the number of neurons per layer (widths). The same conclusion on the stability test is true when we investigate each of the hidden layers of the network.

Figure 4

 4 Figure 4.4(b) shows the results of the stability test for CIFAR10. As can be seen from

Figure 4

 4 Figure 4.5: Estimation of α for varying widths and depths in FCN. The curves in the left figures correspond to different depths, and the ones on the right figures correspond to widths.

Figure 4 . 6 :

 46 Figure 4.6: Estimation of α for varying widths and depths in FCN, dataset MNIST.

Figure 4

 4 Figure 4.7: The accuracy and α of the CNN for varying widths.

Figure 4 .

 4 Figure 4.7 shows the results for the CNN. In this figure, we also depict the train and test accuracy, as well as the tail-index that is estimated on the test set. These results

Figure 4 . 8 :

 48 Figure 4.8: The iteration-wise behavior of of α for the FCN.

Figure 4 . 9 :

 49 Figure 4.9: Estimation of α with an FCN on MNIST.

Figure 4

 4 Figure 4.10: Test error and tail-index in accordance with the change of η/b ratio.

Figure 4 .

 4 Figure 4.10(a) and Figure 4.10(b) visualize the results. The estimated α and the test error are averaged among the candidates with the same η/b, at the last iteration of the training process. We ignore the particular values of test error but rather focus on the way certain choices affect the trends in the system and the behavior of SGD dynamics.

Figure 4 .

 4 Figure 4.11: Test error, estimated α in accordance with the change of η/b ratio.

Figure 6 . 1 :

 61 Figure 6.1: Synthetic experiments. The vertical axis is the quantity |τ 0,a -τ 0,a |.

Figure 6

 6 Figure 6.2: Results of the neural network experiments.

  this inequality for drift b makes the first part of assumption A5 hold.

Lemma 14 .

 14 For the function b defined in Lemma 1, we have, for w ∈ R d , b(w) ≤ M w γ + (B + L), w, b(w) ≤ (L -m) w 1+γ + (b + L). Proof. From assumption A6, it implies that b(w) ≤ c α ∇f (w) + L. Then, by Lemma 13, b(w) ≤ M w γ + (B + L). Next, by Cauchy-Schwarz inequality and assumption A6, we have w, b(w) + c α ∇f (w) ≤ w L. Then, by assumption A3, w, b(w) ≤ -c α w, ∇f (w) + w L ≤ -m w 1+γ + b + w L ≤ -m w 1+γ + b + ( w 1+γ + 1)L =(L -m) w 1+γ + (b + L).

  -B(kη) ≥ u] ≤ 1 -1 -Cde -u 2 /(dη) K and P[ max k∈[0,...,K-1] sup t∈[kη,(k+1)η] L α (t) -L α (kη) ≥ u] ≤ 1 -1 -C α d 1+α/2 ηu -α K .Finally, we get:P[(w η , . . . , w Kη ) ∈B c ] /4] + P[e M η η(B + M ) ≥ ξ/4]

Lemma 18 .

 18 For any u > 0, η > 0 and K ∈ N * , there exist constants C and C α -B(kη) ≥ u] ≤ Cde -u 2 /(dη) . max k∈[0,...,K-1]P[ sup t∈[kη,(k+1)η] L α (t) -L α (kη) ≥ u] ≤ C α d 1+α/2 ηu -α .Proof. To prove the results, we begin with the known results for Brownian motion and α-stable Lévy motion:P[|[B(1)] i | ≥ u] ≤ Ce -u 2 , P[|[L α (1)] i | ≥ u] ≤ C α u -α ,where C and C α are positive constants, [B(1)] i and [L α (1)] i denote the i-th component of the motions respectively, for i from 1 to d. By reflection principle for Brownian motion and α-stable Lévy motion, we haveP[ sup t∈[kη,(k+1)η] |[B(t) -B(kη)] i | ≥ u] ≤ 2P[|[B(η)] i | ≥ u] = 2P[|[B(1)] i | ≥ u/η 1/2 ], P[ sup t∈[kη,(k+1)η] |[L α (t) -L α (kη)] i | ≥ u] ≤ 2P[|[L α (η)] i | ≥ u] = 2P[|[L α (1)] i | ≥ u/η 1/α ]. Since sup t∈[kη,(k+1)η] B(t) -B(kη) 2 ≤ d i=1 sup t∈[kη,(k+1)η] |[B(t) -B(kη)] i | 2 , we have P[ sup t∈[kη,(k+1)η] B(t) -B(kη) ≥ u] =P[ sup t∈[kη,(k+1)η] B(t) -B(kη) 2 ≥ u 2 ] t) -B(kη)] i | 2 ≥ u 2 /d] ] i | ≥ u/(dη) 1/2 ] ≤Cde -u 2 /(dη) .Similarly, we haveP[ sup t∈[kη,(k+1)η] L α (t) -L α (kη) ≥ u] ≤ d i=1 P[|[L α (1)] i | ≥ u/(d 1/2 η 1/α )] ≤C α d 1+α/2 ηu -α .The constants C and C α do not depend on k, hence we have the conclusion.Lemma 19. The following estimates hold: -B(kη) ≥ u] ≤ 1 -1 -Cde -u 2 /(dη) K ,

  k∈[0,...,K-1] sup t∈[kη,(k+1)η]

  25) is obtained from (a + b) γ ≤ a γ + b γ since γ ∈ (0, 1) and a, b ≥ 0. Since 2γ > 1, we have by Corollary 9 E T 2 2γ =E ε(sjη) 1/α L α (1)) 2γ ≤(εη 1/α ) 2γ E L α (1) 2γ ≤(εη 1/α ) 2γ d 2γ 2 2γ Γ((1 + 2γ)/2)Γ(1 -2γ/α) Γ(1/2)Γ(1 -γ) ,and by Corollary 10,E T 3 2γ =E ε(sjη) 1/2 B(1)) 2γ≤(εη 1/2 ) 2γ E B(1) 2γ

  ) in[START_REF] Winkelbauer | Moments and absolute moments of the normal distribution[END_REF].Corollary 10. Let B(1) be a d-dimensional vector whose components are i.i.d standard Gaussian random variable. The quantity E B(1) λ is finite for λ > -1. For details, we have(a) If 1 < λ < α, then E B(1) λ ≤ d λ 2 λ/2 Γ λ+1 Since B(1), by definition, is a d-dimensional vector whose components are i.i.d standard Gaussian random variable B i (1) for i ∈ {1, . . . , d}, i (1)| λ = d 2 λ/2 Γ λ+1 2 √ π .Lemma 23. For 0 < η ≤ m M 2 and s ∈ [jη, (j + 1)η), we have the following estimates:(a) If 1 < λ < α then If 0 ≤ λ ≤ 1 then E ŵjη λ ≤E ŵ0 λ + j (2η(b + m)) Let us denote the value E L α (1) λ by l α,λ,d < ∞ and the value E B(1) λ by b λ,d < ∞. Starting from ŵ(j+1)η = ŵjη -η∇f ( ŵjη ) + εη λ ≤ E ŵjη -η∇f ( ŵjη ) λ + ε λ η λ α E L α (1) λ + ε λ η λ 2 E B(1) λ .(7.37)Consider the first term on the right side:ŵjη -η∇f ( ŵjη ) λ = ŵjη -η∇f ( ŵjη ) 2× λ 2 = ŵjη 2 -2η ŵjη , ∇f ( ŵjη + η 2 ∇f ( ŵjη assumption A12 and Lemma 21. For 0 < η ≤ m M 2 , 2ηm( ŵjη 1+γ + 1) ≥ 2η 2 M 2 ŵjη 2γ . (since 1 + γ > 2γ and ηm > η 2 M 2 )Using this inequality we haveŵjη -η∇f ( ŵjη ) λ ≤ ŵjη 2 + 2η(b + m) + 2η 2 B 2For the case where 0 ≤ λ ≤ 1, by (7.37) and (7.39),E ŵ(j+1)η λ ≤ E ŵjη λ + (2η(b + m)) ) λ + ε λ η λ α l α,λ,d + ε λ η λ 2 b λ,d ≤ E ŵ0 λ + (j + 1) (2η(b + m)) ) λ + ε λ η λ α l α,λ,d + ε λ η λ 2 b λ,d .By using Corollary 9 and Corollary 10, we have the desired results.

Figure 6

 6 Figure 6.1(a) d = 10, α ∈ {1.2, 1.4, 1.6, 1.8}, ε = 0.1, σ = 1, a = 4 × 10 -4 .

Figure 6

 6 Figure 6.1(b) d = 10, α ∈ {1.2, 1.4, 1.6, 1.8}, ε ∈ {10 -3 , 10 -2 , 10 -1 , 10}, σ = 1, a = 4 × 10 -6 .

Figure 6

 6 Figure 6.1(c) d = 10, α ∈ {1.2, 1.4, 1.6, 1.8}, ε = 0.1, σ ∈ {10 -2 , 10 -1 , 1, 10}, a = 4 × 10 -5 .

Figure 6

 6 Figure 6.1(d) d ∈ {10, 40, 70, 100}, α ∈ {1.2, 1.4, 1.6, 1.8}, ε = 0.1, σ = 1, a = 4 × 10 -4 .

  Theorem 15. Under assumptions A10 and A11 we have, for0 < η ≤ min{1, m M 2 }, KL(μ t , µ t ) ≤ ≤ M 2 3 γ 2ε 2 σ 2 kη CM 2γ η 2γ E ŵ0

									], we obtain
			KL(μ t , µ t ) =	1 2ε 2 σ 2	k-1 j=0	(j+1)η jη	E ∇f ( ŵs ) -∇f ( ŵ(jη)) 2 ds	(7.19)
			≤	M 2 2ε 2 σ 2	k-1 j=0	(j+1)η jη	E ŵs -ŵjη	2γ ds.	(7.20)
	Next, we prove the following theorem.
					2γ 2
	+	k -1 2	(2η(b + m)) γ 2	+ 2 γ 2	(ηB) 2γ 2	+ ε 2γ 2 η	2γ 2 α d	2 2γ 2 Γ((1 + 2γ 2 )/2)Γ(1 -2γ 2 /α) Γ(1/2)Γ(1 -γ 2 )
	+ ε 2γ 2 η γ 2 d 2 γ 2 Γ 2γ 2 +1 2 √ π	+	B 2 M 2 + (εη 1/α ) 2γ d 2γ 2 2γ Γ((1 + 2γ)/2)Γ(1 -2γ/α) Γ(1/2)Γ(1 -γ)

  ≤(εη 1/2 ) 2γ d 2γ 2 γ Γ 2γ+1 the equivalence of p -norms and ŵ jη is the concatenation of ŵjη and

	where we used B 1/γ M 1/γ . We then obtain
					E T 1	2γ ≤ CM 2γ η 2γ E ŵ -jη	2γ 2 γ	(7.32)
							≤ CM 2γ η 2γ E ŵ jη	2γ 2 2γ 2	(7.33)
							= CM 2γ η 2γ E ŵjη	2γ 2 2γ 2 +	B 2 M 2	(7.34)
							≤ CM 2γ η 2γ E ŵjη	2γ 2	+	B 2 M 2 .	(7.35)
		By combining the above inequalities and Lemma 23, we obtain
	KL(μ t , µ t ) ≤
	≤	M 2 3 γ 2ε 2 σ 2	k-1 j=0	(j+1)η jη	E T 1	2γ + T 2	2γ + T 3	2γ ds
									2 √ π	,
		By definition, we have
				E T 1	2γ = E (s -jη)∇f ( ŵjη ) 2γ	(7.27)
						≤ η 2γ E ∇f ( ŵjη ) 2γ	(7.28)
						≤ η 2γ E(M ŵjη	γ + B) 2γ	(7.29)
						≤ CM 2γ η 2γ E ŵjη	γ γ +	B 1/γ M 1/γ	γ	2γ	(7.30)
									2γ
						≤ CM 2γ η 2γ E ŵ jη	γ γ	(7.31)
								109

We note that more sophisticated assumptions than (4.1) have been made in terms of the covariance matrix of the Gaussian distribution (e.g. state dependent, anisotropic). However, in all these cases, the resulting distribution is still a Gaussian, therefore the same criticism holds.

The codebase can be found at https://github.com/umutsimsekli/sgd_tail_index.

We observed a similar behavior in under-parametrized FCN; however, did not plot those results to avoid clutter.

Note that this particular scaling is introduced in Jacot-Guillarmod et al.[2018] and it admits slightly larger values of learning rates compared to standard initialization schemes

It is easy to verify that ŵkη = w k for all k ∈ N +[Raginsky et al.,

2017].

Acknowledgements

Firstly, I would like to express my sincere gratitude to my advisors Assoc. Prof.

Umut S ¸imşekli and Prof. Gaël Richard for the continuous support and help during my PhD, and for their patience, motivation, and immense knowledge. Without their guidance, I would not have been able to complete this research and the writing of this thesis.

Besides my advisors, I would like to thank the rest of my thesis committee members: Dr. Olivier Cappé, Assoc. Prof. Zaid Harchaoui, Assoc. Prof. Lenka Zdeborova, Assoc.

Prof. Alain Durmus and Prof. Ali Taylan Cemgil for their insightful comments, valuable suggestions and encouragement.

Bibliography

small for a given k. Therefore, in our strategy, we choose a different path where we do not use the strong convergence of the Euler scheme.

Our proof strategy is inspired by the recent study [START_REF] Tzen | Local optimality and generalization guarantees for the Langevin algorithm via empirical metastability[END_REF], where the authors analyze the empirical metastability of the Langevin equation which is driven by a Brownian motion. However, unlike the Brownian case that [START_REF] Tzen | Local optimality and generalization guarantees for the Langevin algorithm via empirical metastability[END_REF] was based on, some of the tools for analyzing Brownian SDEs do not exist for Lévy-driven SDEs, which increases the difficulty of our task.

We first define a linearly interpolated version of the discrete-time process {w k } k∈N+ , which will be useful in our analysis, given as follows: (6.8) where ŵ ≡ { ŵt } t≥0 denotes the whole process and the drift function b( ŵ) is chosen as follows:

Here, I denotes the indicator function, i.e.

It is easy to verify that ŵkη = w k for all k ∈ N + Dalalyan [2017b], [START_REF] Raginsky | Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis[END_REF].

In our approach, we start by developing a Girsanov-like change of measures [START_REF] Tankov | Financial modelling with jump processes[END_REF] to express the Kullback-Leibler (KL) divergence between µ t and μt , which is defined as follows:

where µ t denotes the law of {w s } s∈[0,t] , μt denotes the law of { ŵs } s∈[0,t] , and dµ t /dμ t is the Radon-Nikodym derivative of µ t with respect to μt . Here, we require A9 for the existence of a Girsanov transform between μt and µ t and for establishing an explicit formula for the transform. In appendix, we show that the KL divergence between µ t and μt can be written as:

(6.9)

While this result has been known for SDEs driven by Brownian motion [START_REF] Øksendal | Applied stochastic control of jump diffusions[END_REF], none of the references we are aware of expressed the KL divergence as in (6.9). We also note that one of the key reasons that allows us to obtain (6.9) is the presence of the Brownian motion in (6.2), i.e. σ > 0. For σ = 0 such a measure transformation cannot be performed [START_REF] Debussche | Existence of densities for stable-like driven SDE's with Hölder continuous coefficients[END_REF].

In the next result, we show that if the step-size is chosen sufficiently small, the KL divergence between µ t and μt is bounded.

Theorem 12. Assume that the conditions A8-A13 hold. Then the following inequality holds:

Appendix

Supplementary materials for Chapter 5

Proof of Lemma 1

Proof. Let q(x, t) be the probability density of x(t). By Proposition 1 in Schertzer et al.

[2001] (see also Section 7 of the same study), the fractional Fokker-Planck equation associated with (5.5) is given as follows:

By combining the above two inequalities, we get

The following lemma is an extension of Lemma 1.2.3 in [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF] to functions with Hölder continuous gradients.

Lemma 7. Under assumption A2, the following inequality holds for any x, y ∈ R d :

By Cauchy-Schwarz inequality and assumption A2, we have

Lemma 8. The normalized factor of π is bounded below, i.e., log

where π denotes the Archimedes' constant.

Supplementary materials for Chapter 6

More details on assumption A 13

In this section, we provide the precise expressions of the constants given in Assumption A13. For a given δ > 0, t = Kη, and for some C > 0, the step-size satisfies the following condition:

, where ε is as in (6.3), the constants m, M, b are defined by A10-A12 and

Proof of Theorem 11

Proof. Note that (w 1 , . . . , w K ) ∈ A is equivalent to τ0,a (ε) > K. Hence, from Lemma 20, the remaining task is to upper-bound P[(w η , . . . , w Kη ) ∈ A]:

and to lower-bound it:

By Lemma 17, the final result follows.

Lemma 17. There exist constants C, C 1 and C α such that:

Lemma 20. Suppose that assumptions A10 and A11 hold. Then, for any δ > 0, we have:

By optimal coupling between two probability measure [START_REF] Lindvall | Lectures on the coupling method[END_REF], Theorem 5.2), there exists a coupling M of (w s ) 0≤s≤Kη and ( ŵs ) 0≤s≤Kη such that

where T V denotes the total variation distance. By Pinsker's inequality, we also have

Then, .

By Theorem 12, we have the desired inequalities.

Proof of Theorem 12

First, we derive a Girsanov-type change of measure [START_REF] Øksendal | Applied stochastic control of jump diffusions[END_REF],

Tankov [2003] for the SDE considered in (6.2). Let P denote the law of w t and Q be an equivalent measure defined by (7.17) where F T denotes the filtration up to time T . Then the process B φ defined by B φ (t) = B(t) -t 0 φ s ds is a Q-Brownian motion. With the choice of φ t given in A9, we see that