Mean-Square Analysis of Discretized It\^o Diffusions for Heavy-tailed Sampling

TL;DR

This paper investigates the complexity of sampling from heavy-tailed distributions using discretized Itô diffusions, providing mean-square analysis and explicit iteration bounds under minimal variance assumptions.

Contribution

It introduces a mean-square analysis framework for heavy-tailed sampling via discretized Itô diffusions, with explicit iteration complexity bounds and gradient estimation methods.

Findings

Established iteration complexity bounds for heavy-tailed distributions

Derived explicit moment bounds for heavy-tailed targets

Extended results to gradient estimation via Gaussian smoothing

Abstract

We analyze the complexity of sampling from a class of heavy-tailed distributions by discretizing a natural class of It\^o diffusions associated with weighted Poincar\'e inequalities. Based on a mean-square analysis, we establish the iteration complexity for obtaining a sample whose distribution is $\epsilon$ close to the target distribution in the Wasserstein-2 metric. In this paper, our results take the mean-square analysis to its limits, i.e., we invariably only require that the target density has finite variance, the minimal requirement for a mean-square analysis. To obtain explicit estimates, we compute upper bounds on certain moments associated with heavy-tailed targets under various assumptions. We also provide similar iteration complexity results for the case where only function evaluations of the unnormalized target density are available by estimating the gradients using a Gaussian smoothing technique. We provide illustrative examples based on the multivariate $t$-distribution.

Equations

$$d{X}_t=\nabla \log \pi (X_t)dt+\sqrt{2}d{B}_t,$$

$$(\beta +1)Va{r}_{{\pi }_{\beta }}(g)\le {\int }_{{\mathbb{R}}^d}\frac{⟨({\nabla }^{2}V{)}^{-1}\nabla G,\nabla G⟩}{V}d{\pi }_{\beta }+\frac{d}{\beta -d}\left({\int }_{{\mathbb{R}}^d}gd{\pi }_{\beta }\right).$$

$$d{X}_t=-(\beta -1)\nabla V(X_t)dt+\sqrt{2V(X_t)}d{B}_t,$$

$$W_2(\mu ,\nu ):=\underset{\zeta \in C(\mu ,\nu )}{\inf }{\left({\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}|x-y{|}^2\zeta (dx,dy)\right)}^{\frac{1}{2}}.$$

$${\chi }^2(\nu |\mu ):={\int }_{{\mathbb{R}}^d}{\left(\frac{\nu (dx)}{\mu (dx)}-1\right)}^2\mu (dx).$$

$$\Gamma (z):={\int }_0^{\infty }{t}^{z-1}{e}^{-t}dt,\forall z>0,\text{and}B(x,y):={\int }_0^1{t}^{x-1}(1-t{)}^{y-1}dt,\forall x,y>0.$$

$$\mathcal{L}\psi =-(\beta -1)⟨\nabla V,\nabla \psi ⟩+V\Delta \psi .$$

$${\partial }_t{\rho }_t=\nabla \cdot \left(\beta {\rho }_t\nabla V+V\nabla {\rho }_t\right)=\nabla \cdot \left({\rho }_{t}V\nabla \log \frac{{\rho }_t}{{\pi }_{\beta }}\right).$$

$$\frac{⟨({\nabla }^{2}V{)}^{-1}(x)\nabla V(x),\nabla V(x)⟩}{V(x)}\le C_V.$$

$$\nabla V(x)=\frac{{\varphi }^{\prime }(|x|)}{|x|}x,\text{and}{\nabla }^{2}V=\left({\varphi }^{\prime \prime }(|x|)-\frac{{\varphi }^{\prime }(|x|)}{|x|}\right)\frac{x\otimes x}{|x{|}^2}+\frac{{\varphi }^{\prime }(|x|)}{|x|}{I}_d,$$

$${\varphi }^{\prime }(r)\le ({\varphi }^{\prime \prime }(r)r)\wedge (C_V\varphi (r)/r),\text{ for all }r\ge 0.$$

$$Va{r}_{{\pi }_{\beta }}(\varphi )$$

$$(\beta +1){\int }_{{\mathbb{R}}^d}g(x{)}^2{\pi }_{\beta }(x)dx\le {\int }_{{\mathbb{R}}^d}\frac{⟨({\nabla }^{2}V{)}^{-1}(x)\nabla (gV)(x),\nabla (gV)(x)⟩}{V(x)}{\pi }_{\beta }(x)dx.$$

$$⟨A(u+v),(u+v)⟩\le r⟨Au,u⟩+\frac{r}{r-1}⟨Av,v⟩,u,v\in {\mathbb{R}}^d,r>1,$$

$$(\beta +1){\int }_{{\mathbb{R}}^d}g(x{)}^2{\pi }_{\beta }(x)dx$$

$$+\frac{r}{r-1}{\int }_{{\mathbb{R}}^d}\frac{⟨({\nabla }^{2}V{)}^{-1}(x)V(x)\nabla g(x),V(x)\nabla g(x)⟩}{V(x)}{\pi }_{\beta }(x)dx.$$

$$(\beta +1){\int }_{{\mathbb{R}}^d}g(x{)}^2{\pi }_{\beta }(x)dx$$

$$+\frac{r}{r-1}{\int }_{{\mathbb{R}}^d}⟨V(x)({\nabla }^{2}V{)}^{-1}(x)\nabla g(x),\nabla g(x)⟩{\pi }_{\beta }(x)dx,$$

$${\int }_{{\mathbb{R}}^d}g(x{)}^2{\pi }_{\beta }(x)dx$$

$$Va{r}_{{\pi }_{\beta }}(\varphi )$$

$$\frac{⟨\nabla V(x),\nabla V(x)⟩}{V(x)}\le \alpha C_V.$$

$${\text{Var}}_{{\pi }_{\beta }}(\varphi )$$

$${\chi }^2({\rho }_t|{\pi }_{\beta })\le \exp \left(-2\alpha {\left(\sqrt{\beta +1}-\sqrt{C_V}\right)}^{2}t\right){\chi }^2({\rho }_0|{\pi }_{\beta }).$$

$$\frac{d}{dt}{\chi }^2({\rho }_t|{\pi }_{\beta })$$

$$=2{\int }_{{\mathbb{R}}^d}{\partial }_t{\rho }_t(x)\left(\frac{{\rho }_t(x)}{{\pi }_{\beta }(x)}-1\right)dx$$

$$=-2{\int }_{{\mathbb{R}}^d}⟨\nabla \left(\frac{{\rho }_t}{{\pi }_{\beta }}\right)(x),\nabla \log \left(\frac{{\rho }_t}{{\pi }_{\beta }}\right)(x)⟩V(x){\rho }_t(x)dx$$

$$=-2{\int }_{{\mathbb{R}}^d}V(x){\left|\nabla \left(\frac{{\rho }_t}{{\pi }_{\beta }}\right)(x)\right|}^2{\pi }_{\beta }(x)dx.$$

$$\frac{d}{dt}{\chi }^2({\rho }_t|{\pi }_{\beta })$$

$$=-2\alpha {\left(\sqrt{\beta +1}-\sqrt{C_V}\right)}^2{\chi }^2({\rho }_t|{\pi }_{\beta }).$$

$$x_{k+1}=x_k-h(\beta -1)\nabla V(x_k)+\sqrt{2hV(x_k)}{\xi }_{k+1},$$

$$\delta :=\frac{\beta -1-\frac{1}{4}{C}_{V}d}{\frac{1}{4}{C}_{V}d}>0.$$

$$h<\min \left(\frac{1}{4(\beta -1)L},\frac{2\delta }{3(1+\delta )\alpha (\beta -1)}\right),$$

Full text

1 Introduction

The problem of sampling from a given target density $\pi:\mathbb{R}^{d}\to\mathbb{R}$ arises in a wide variety of problems in statistics, machine learning, operations research and applied mathematics. Markov chain Monte Carlo (MCMC) algorithms are a popular class of algorithms for sampling (Robert and Casella, 1999, Andrieu et al., 2003, Hairer et al., 2006, Brooks et al., 2011, Meyn and Tweedie, 2012, Leimkuhler and Matthews, 2016, Douc et al., 2018); a widely used approach in this domain is to discretize an Itô diffusion that has the target as its stationary density. A popular choice of diffusion is the overdamped Langevin diffusion,

[TABLE]

where $B_{t}$ is a $d$-dimensional Brownian motion. For example, the Unadjusted Langevin Algorithm (Rossky et al., 1978), the Metropolis Adjusted Langevin Algorithm (Roberts and Tweedie, 1996, Roberts and Rosenthal, 1998) and the proximal sampler (Titsias and Papaspiliopoulos, 2018, Lee et al., 2021, Vono et al., 2022) arise as different discretizations of (1). Under light-tailed assumptions, i.e. when the density $\pi$ has exponentially fast decaying tails, the diffusion $X_{t}$ in (1) converges exponentially fast to $\pi$ as its stationary density, which motivates the use of discretizations of (1) as practical algorithms for sampling. In the last decade, the non-asymptotic iteration complexity of various discretizations have been well-explored, thereby providing a relatively comprehensive story of sampling from light-tailed densities.

Motivated by applications in robust statistics (Kotz and Nadarajah, 2004, Jarner and Roberts, 2007, Kamatani, 2018), multiple comparison procedures (Genz et al., 2004, Genz and Bretz, 2009), Bayesian statistics (Gelman et al., 2008, Ghosh et al., 2018), and statistical machine learning (Balcan and Zhang, 2017, Nguyen et al., 2019, Şimşekli et al., 2020, Diakonikolas et al., 2020), in this work, we are interested in sampling from densities that have heavy-tails, for example, those with tails that are polynomially decaying. When the target density $\pi$ is heavy-tailed, the solution to (1) does not converge exponentially to its stationary density in various metrics of interest. Indeed, Theorem 2.4 by Roberts and Tweedie (1996) shows that if $|\nabla\log\pi(x)|\to 0$ as $|x|\to\infty$, then the solution to (1) is not exponentially ergodic. In the other direction, standard results in the literature, for example Wang (2006), Bakry et al. (2014) show that the solution to (1) converging exponentially fast to its equilibrium density in the $\chi^{2}$ metric, is equivalent to the density $\pi$ satisfying the Poincaré inequality, which in turn requires $\pi$ to have exponentially decaying tails. Furthermore, when $\pi$ has polynomially decaying tails, the convergence is only sub-exponential or polynomial (Wang, 2006, Chapter 4). Consequently, the algorithms obtained as discretizations of the Langevin diffusion in (1) are suited to sampling only from light-tailed exponentially decaying densities, and are rather inefficient for sampling from heavy-tailed densities.

Our approach to heavy-tailed sampling is hence based on discretizing certain natural Itô diffusions that arise in the context of the following Weighted Poincaré inequality (Blanchet et al., 2009, Bobkov and Ledoux, 2009). Such inequalities could be considered generalizations of the Brascamp-Lieb inequality (established for the class of log-concave densities) to a class of heavy-tailed densities.

Theorem 1** (Weighted Poincaré Inequality; Theorem 2.3 in Bobkov and Ledoux (2009)).**

Let the target density be of the form $\pi_{\beta}\propto V^{-\beta}$ with $\beta>d$ and $V\in\mathcal{C}^{2}(\mathbb{R}^{d})$ positive, convex and with $(\nabla^{2}V)^{-1}(x)$ well-defined for all $x\in\mathbb{R}^{d}$. For any smooth and $\pi_{\beta}$-integrable function $g$ on $\mathbb{R}^{d}$ and $G=Vg$,

[TABLE]

A canonical example of a heavy-tailed density that satisfies the conditions in Theorem 1, and hence (2), is the multivariate $t$-distribution. In particular, we consider the following Itô diffusion process

[TABLE]

where $(B_{t})_{t\geq 0}$ is a standard Brownian motion in $\mathbb{R}^{d}$. The Itô diffusion in (3) converges exponentially fast to the target $\pi_{\beta}$ in the $\chi^{2}$-divergence as long as it satisfies the Weighted Poincaré inequality and additional mild assumptions; see Proposition 1 for details. Hence, we study the oracle complexity of the Euler-Maruyama discretization of (3), for sampling from heavy-tailed densities. Our proofs are based on mean-square analysis techniques, a popular technique to analyze numerical discretizations of stochastic differential equations; see, for example, Milstein and Tretyakov (2004) for an overview. Our results in this paper pushes mean-square analysis to its limits; the heavy-tailed densities we consider invariably need to have only finite variance, which is the minimum requirement when using this technique.

1.1 Our Contributions

In this work, we make the following contributions:

• •

In Theorem 2, we provide upper bounds on the number of iterations required by the Euler-Maruyama discretization of (3) to obtain a sample that is $\epsilon$-close in the Wasserstein-2 metric to the target density. The established bounds are in terms of certain (first and second-order) moments of the target density $\pi$. Our proof technique is based on a mean-squared analysis; we demonstrate that for the case of multivariate $t$-distributions, our analysis is non-vacuous as long as the density has finite variance, a necessary condition to carry out the mean-squared analysis.

• •

While the result in Theorem 2 assumes access to the exact gradient of the unnormalized target density function (referred to as the first-order setting), in Theorem 3, we analyze the case when the gradient is estimated based on function evaluations (the zeroth-order setting) based on a Gaussian smoothing technique.

• •

We provide several illustrative examples highlighting the differences between the results in the first and the zeroth-order setting. Specifically, in Section 5 we show that for the multivariate $t$-distribution with smaller degrees of freedom, (and hence the truly heavy-tailed case) the gradient estimation error is dominated by the discretization error. Whereas, in the case with larger degrees of freedom (and hence the comparatively moderately heavy-tailed case), the discretization error is of comparable order to the gradient estimation error. Hence, the zeroth-order algorithm matches the iteration complexity of the first-order algorithm by using mini-batch gradient estimators.

1.2 Related Work

Non-asymptotic iteration complexity of different discretizations of (1) have been analyzed extensively in the last decade. The analysis of the Unadjusted Langevin Algorithm (ULA) under various light-tailed assumptions was carried out, for example, in Dalalyan (2017), Durmus and Moulines (2017), Dalalyan and Karagulyan (2019), Durmus et al. (2019), Lee et al. (2020), Shen and Lee (2019), He et al. (2020), Chen et al. (2020), Durmus et al. (2019), Dalalyan et al. (2019), Li and Erdogdu (2020), Chen et al. (2020), Chewi et al. (2021a) and references therein. In particular, Vempala and Wibisono (2019), Erdogdu and Hosseinzadeh (2021), Chewi et al. (2021a) analyzed the performance of ULA under various functional inequalities suited to light-tailed densities. Furthermore, the recent work of Balasubramanian et al. (2022) analyzed the performance of (averaged) ULA for target densities that are only Hölder continuous, albeit in the weaker Fisher information metric.

Several works, for example, Dwivedi et al. (2019), Chewi et al. (2021b), Wu et al. (2022), analyzed the Metropolis-Adjusted Langevin Algorithm (MALA) in light-tailed settings. The proximal sampler algorithm was analyzed under various light-tailed assumptions in Lee et al. (2021), Chen et al. (2022). The iteration complexity of the widely used Hamiltonian Monte Carlo algorithm and discretizations of underdamped Langevin diffusions were analyzed, for example, in Dalalyan and Riou-Durand (2020), Bou-Rabee et al. (2020), Chen et al. (2020), Ma et al. (2021), Monmarché (2021), Cao et al. (2021), Wang and Wibisono (2022), Chen and Vempala (2022). We also refer interested readers to Lu and Wang (2022), Ding and Li (2021) for non-asymptotic analyses of other MCMC algorithms used in practice in light-tailed settings.

In the context of heavy-tailed sampling, Kamatani (2018) considered the scaling limits of appropriately modified Metropolis random walk in an asymptotic setting. Johnson and Geyer (2012) proposed a variable transformation method in the context of Metropolis Random Walk algorithms. Here, the heavy-tailed density is converted into a light-tailed one based on certain invertible transformations so that one can leverage the rich literature on light-tailed sampling algorithms. Similar ideas were also examined recently in Yang et al. (2022). It is also worth highlighting that Deligiannidis et al. (2019), Durmus et al. (2020) and Bierkens et al. (2019) used the transformation approach for proving asymptotic exponential ergodicity of bouncy particle and zig-zag samplers respectively, in the heavy-tailed setting. We also point out the recent works of Andrieu et al. (2021a) and Andrieu et al. (2021b) that establish similar sub-exponential ergodicity results for other sampling methods such as the piecewise deterministic Markov process Monte Carlo, independent Metropolis-Hastings sampler and pseudo-marginal methods in the polynomially heavy-tailed setting. The works of Şimşekli et al. (2020), Huang et al. (2021) and Zhang and Zhang (2022) established exponential ergodicity results for diffusions driven by $\alpha$-stable processes with heavy-tailed densities as its equilibrium in the continuous-time setting. However, the problem of obtaining convergence results for practical discretizations of these diffusions is still largely open.

The literature on non-asymptotic oracle complexity analysis of heavy-tailed sampling is extremely limited. Chandrasekaran et al. (2009) considered the iteration complexity of Metropolis random walk algorithm for sampling from $s$-concave distributions. He et al. (2022) considered ULA on a class of transformed densities (i.e., the heavy-tailed density is transformed to a light-tailed one with an invertible transformation, similar to Johnson and Geyer (2012)) and established non-asymptotic oracle complexity results. However, they focused mainly on the case of isotropic densities. Li et al. (2019) analyzed a class of discretizations of general Itô diffusions that admit heavy-tailed equilibrium densities. A detailed comparison to Li et al. (2019) is provided in Section 5.

The recent works by Hsieh et al. (2018), Zhang et al. (2020), Chewi et al. (2020), Ahn and Chewi (2021), Jiang (2021), Li et al. (2022) also considered sampling based on discretizations of the Mirror Langevin diffusions. The above-mentioned works mainly focus on sampling from constrained densities. The continuous-time convergence is analyzed typically under the so-called mirror Poincaré inequalities which are generalizations of the Brascamp-Lieb inequalities in a different direction compared to the Weighted Poincaré inequalities. The discretization analysis by Li et al. (2022) is based on mean-squared analysis.

As mentioned previously, our work leverages the literature on weighted functional inequalities, that are satisfied by heavy-tailed densities. The weighted Poincare inequality was introduced in Blanchet et al. (2009) and Bobkov and Ledoux (2009), and using an extension of the Brascamp-Lieb inequality, is shown to hold for the class of $s$-concave densities. We also refer the interested reader to Cattiaux et al. (2010, 2011), Bonnefont et al. (2016), Cordero-Erausquin and Gozlan (2017), Cattiaux et al. (2019) for various extensions and improvements of the works of Blanchet et al. (2009) and Bobkov and Ledoux (2009).

1.3 Notation

We use the following notation throughout the rest of the paper.

• •

$\langle\cdot,\cdot\rangle$ denotes the Euclidean inner product and $|\cdot|$ denotes the Euclidean norm.

• •

For two matrices $A$ and $B$, $A\preceq B$ means that $B-A$ is positive semi-definite. The 2-norm of any $d\times d$ matrix $A$ is denoted as $\left\lVert A\right\rVert_{2}$. $I_{d}$ is the $d\times d$ identity matrix.

• •

$\Delta$ denotes the Laplacian, and $\nabla$ denotes the gradient of a given function.

• •

$\mathcal{C}^{2}(\mathbb{R}^{d})$ refers to the set of all real functions on $\mathbb{R}^{d}$ that are twice continuously differentiable. $\mathcal{C}^{2}_{c}(\mathbb{R}^{d})$ refers to the set of all functions in $\mathcal{C}^{2}(\mathbb{R}^{d})$ with compact support.

• •

The Wasserstein-2 distance between two probability measures on $\mathbb{R}^{d}$, $\mu$ and $\nu$ is given by

[TABLE]

where $C(\mu,\nu)$ is the set of all measures on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ whose marginals are $\mu$ and $\nu$ respectively.

• •

The $\chi^{2}$ divergence from a probability measure $\nu$ to a probability measure $\mu$ is defined as

[TABLE]

• •

The gamma and beta functions are given by:

[TABLE]

• •

For two positive quantities $f(d),g(d)$ depending on $d$, we define $f(d)=O(g(d))$ if there exists a constant $C>0$ such that $f(d)\leq Cg(d)$ for all $d>1$. We define $f(d)=\Theta(g(d))$ if there exist constants $C_{1},C_{2}>0$ such that $C_{1}g(d)\leq f(d)\leq C_{2}g(d)$ for all $d>1$. We use $\tilde{O}$ to hide $\log$ factors in the $O$ notation.

1.4 Organization

In Section 2, we first establish the exponential ergodicity of the Itô diffusion in (3) under certain assumptions that are favorable for the discretization analysis. We next provide our main results on the non-asymptotic oracle complexity of the Euler-Maruyama discretization of (3). In Section 3, we provide moment computations in the heavy-tailed setting that are required to obtain explicit rates from the results in Section 2. In Section 4, we provide an extension of our results to the zeroth-order setting. In Section 5 we provide several illustrative examples. We discuss further implications of our assumptions in Section 6. The proofs are provided in Section 7 and in Appendices A, B and C.

2 Itô Discretizations and Weighted Poincare inequalities

In this section, our goal is to analyze the Itô diffusion in (3) which admits a specific class of heavy-tailed densities as its stationary density. Let $X_{0}$ follow distribution $\rho_{0}$ and denote the distribution of $X_{t}$ by $\rho_{t}$ for all $t\geq 0$. For any function $\psi\in\mathcal{C}^{2}_{c}(\mathbb{R}^{d})$, the infinitesimal generator of (3) is given by

[TABLE]

Hence, the Fokker-Planck equation corresponding to (3) is

[TABLE]

It follows that, under the conditions in Theorem 1, $\pi_{\beta}\propto V^{-\beta}$ is the unique stationary density of (3). We next examine the convergence properties of (3) to its stationary density. To do so, we introduce the following assumption.

Assumption 1**.**

There exists a positive constant $C_{V}$ such that, for all $x\in\mathbb{R}^{d}$,

[TABLE]

When $V$ is radially symmetric, i.e., when $V(x):=\phi(|x|)$ for some $\phi\in\mathcal{C}^{2}(\mathbb{R_{+}})$, the condition in Assumption 1 simplifies as follows. Note that

[TABLE]

where $\otimes$ denotes outer-product. Hence, it follows that it is sufficient for $\phi$ to satisfy

[TABLE]

For example, this property holds with $C_{V}=p$ if $\phi$ is a $p$-order polynomial with $p\geq 2$ and non-negative coefficients.

We next provide the following corollary to Theorem 1, motivated by the discussion in Section 2 of Bobkov and Ledoux (2009).

Corollary 1**.**

Consider the setting of Theorem 1 and suppose further that Assumption 1 holds with $C_{V}\in(0,\beta+1)$, then for any smooth, $\pi_{\beta}$-integrable function, $\phi$ on $\mathbb{R}^{d}$,

[TABLE]

Proof.

We start from (2), assume that $\int_{\mathbb{R}^{d}}gd\pi_{\beta}=0$. Then (2) could be rewritten as

[TABLE]

Now, note that we have the following elementary bound

[TABLE]

for any arbitrary positive definite symmetric matrix $A\in\mathbb{R}^{d\times d}$. Hence, we obtain

[TABLE]

Invoking the condition in Assumption 1, we further obtain

[TABLE]

which then implies that, for any $r\in(1,(\beta+1)/C_{V})$,

[TABLE]

With the choice of $r\coloneqq\sqrt{\frac{\beta+1}{C_{V}}}>1$, we get that for all $g$ such that $\int gd\pi_{\beta}=0$, and

[TABLE]

For all general $\phi$, letting $g=\phi-\int\phi d\pi_{\beta}$, we get

[TABLE]

∎

When $V$ is strongly convex, Assumption 1 holds under the following sufficient condition.

Assumption 2**.**

The function $V:\mathbb{R}^{d}\to(0,\infty)$ is twice continuously differentiable and V satisfies

• (1)

$V$* is $\alpha$-strongly convex, i.e. $\nabla^{2}V(x)\succeq\alpha I_{d}$ for all $x\in\mathbb{R}^{d}$.*

• (2)

There exists a positive constant $C_{V}$ such that, for all $x\in\mathbb{R}^{d}$,

[TABLE]

The following result follows immediately from Assumption 2.

Lemma 1**.**

Let $\beta>d$. If Assumption 2 holds with $C_{V}\in(0,\beta+1)$, then for any smooth, $\pi_{\beta}$ integrable function $\phi$ on $\mathbb{R}^{d}$, we have

[TABLE]

With (7), we can show the exponential decay in $\chi^{2}$-divergence along (3). The proof of the following proposition is standard and we include it here for completeness.

Proposition 1**.**

Under the conditions in Lemma 1, for $(X_{t})$ following diffusion (3) with $\rho_{t}$ being the distribution of $X_{t}$, we have

[TABLE]

Proof of Proposition 1.

First we can calculate the derivative of $\chi^{2}(\rho_{t}|\pi)$ via (5),

[TABLE]

According to (7), we get

[TABLE]

Finally, (8) follows from Gronwall’s inequality. ∎

The above result shows that for the class of $\pi_{\beta}$ satisfying Assumption 2, the Itô diffusion in (3), converges exponentially fast to its stationary density. Hence, time-discretizations of (3) provide a practical way of sampling from that class of densities. The Euler-Maruyama discretization to (3) is given by

[TABLE]

where $h>0$ is the step size and $\{\xi\}_{k=1}^{\infty}$ is a sequence of i.i.d. standard Gaussian random vectors in $\mathbb{R}^{d}$. We now present our main result on the iteration complexity of (9) for sampling from $\pi_{\beta}$. We state our discretization result, based on a mean-square analysis, in the $W_{2}$ metric. In particular, we highlight that Proposition 1 requires that condition that $\beta>d$, in addition to Assumption 2, whereas Theorem 2 below, does not. In Section 6, we revisit these conditions and provide additional insights. Obtaining convergence results in the stronger $\chi^{2}$-divergence is left as future work.

Theorem 2**.**

Let $V$ be gradient-Lipschitz with parameter $L>0$, and satisfying Assumption 2 with

[TABLE]

Let $(x_{k})_{k=0}^{\infty}$ be generated from (9) with $\nu_{k}$ denoting the distribution of $x_{k}$, for all $k\geq 0$. Then with the step-size,

[TABLE]

the decay of Wasserstein-2 distance along the Markov chain $(x_{k})_{k=0}^{\infty}$ can be described by the following equation: For all $k\geq 1$,

[TABLE]

with $A,B$ and $C$ given respectively in (51), (52) and (53).

Remark 1** (Constant $\delta$).**

We now motivate the definition and the condition on the constant $\delta$ based on exponential contractivity arguments.

Definition 1** (Exponential contractivity).**

*Let $X_{t}$, $Y_{t}$ be two different solutions to the same stochastic differential equation (SDE) with initial conditions $x,y$ respectively. We say the SDE is $W_{2}$-exponential contractive if there exists a constant $\kappa>0$, such that *

[TABLE]

where by $L(X)$ we refer to the law of $X$.

Uniform dissipativity is a sufficient condition for exponential contractivity (Gorham et al., 2019, Theorem 10). The uniform dissipativity condition for (3) can be represented as

[TABLE]

or equivalently as

[TABLE]

When $V$ satisfies Assumption 2, a sufficient condition for the above uniform dissipativity condition is given by

[TABLE]

or equivalently,

[TABLE]

The sufficient condition coincides with the condition that $\delta>0$ in Theorem 2, which also motivates the assumption in Theorem 2.

Remark 2** (Iteration complexity).**

With Theorem 2, we can calculate the order of the iteration complexity to reach an $\epsilon$-accuracy in Wasserstein-2 distance. With (51),(52),(53), we have

[TABLE]

The above display implies that

[TABLE]

Hence, we get $\frac{C}{A}+\frac{B}{\sqrt{A(2-A)}}<\epsilon/2$ if the step-size $h$ satisfies

[TABLE]

Defining $K_{\epsilon}=\log\left({2W_{2}(\nu_{0},\pi_{\beta})}/{\epsilon}\right)$, we have $W_{2}(\nu_{k},\pi_{\beta})<\epsilon$ for all $k\geq K$ with

[TABLE]

Recall the definition of $\delta$ in (10). The order of $K$ depends on the order of $\delta$. That is, we have the following two cases:

• •

If $\delta=O(1)$ and $\beta=O(d)$, we have that

[TABLE]

• •

If $\delta=O(1/d)$ and $\beta=O(d)$, we have that

[TABLE]

In order to obtain more explicit iteration complexity bounds from Remark 2, it is required to compute bounds on the following two quantities: $\mathbb{E}_{\pi_{\beta}}\left[|\nabla V(X)|^{2}\right]$ and $\mathbb{E}_{\pi_{\beta}}\left[V(X)\right]$.

3 Moment Bounds

In this section, we compute moment bounds under the conditions in Theorem 2.

3.1 An Example: Multivariate $t$-distribution

We first start with the isotropic case.

Proposition 2**.**

Let $\pi_{\beta}=Z_{\beta}^{-1}V^{-\beta}$ with $\beta>d/2+1$, $V(x)=1+|x|^{2}$ and $Z_{\beta}=\int_{\mathbb{R}^{d}}(1+|x|^{2})^{-\beta}dx$. We have

[TABLE]

Proof.

Let $A_{d}(1)$ denote the surface area of the unit sphere in $d$ dimensions. By a standard calculation, we have that, for all $\beta>\frac{d}{2}$,

[TABLE]

where $B$ is the beta function. In the above calculation, the second identity follows from a change to polar coordinates. The third identity follows from a substitution with $R=r^{2}$ and the fourth identity follows from a substitution $u={R}/{(1+R)}$. Therefore for all $\beta>d/2+1$, we have that

[TABLE]

where the fourth identity follows from the property of Beta function, $B(x,y)=\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$ and the fifth identity follows from the property of $\Gamma$ function, $\Gamma(1+z)=z\Gamma(z)$. For the other expectation, we have

[TABLE]

where we apply the same substitutions and properties of Beta functions and Gamma functions in the above calculation. ∎

Remark 3**.**

If $\pi_{\beta}$ is the class of isotropic multivariate $t$-distributions, with the results in Proposition 2, the order of the two expectations in terms of the dimension parameter $d$ is given as follows,

• •

when $\beta>\frac{d}{2}+1$ and $\beta-1-\frac{d}{2}=O(d)$, we have

[TABLE]

• •

when $\beta>\frac{d}{2}+1$ and $\beta-1-\frac{d}{2}=O(1)$, we have

[TABLE]

For a general class of non-isotropic multivariate t-distribution, we consider $\pi_{\beta}=Z_{\beta}^{-1}V^{-\beta}$ with $V(x)=1+x^{T}\Sigma x$ where $\Sigma$ is a strictly positive-definite $d\times d$ matrix. In Roth (2012), it’s been shown that for any $\beta>\frac{d}{2}$, the normalization constant is

[TABLE]

Therefore for any $\beta>\frac{d}{2}+1$, we have

[TABLE]

and

[TABLE]

The above identity implies

[TABLE]

where the second inequality follows from the fact that $\nabla^{2}V(x)=\Sigma$.

Remark 4**.**

If $\pi_{\beta}$ is in the class of non-isotropic multivariate t-distributions, the order of the two expectations in terms of the dimension parameter $d$ is as follows,

• •

when $\beta>\frac{d}{2}+1$ and $\beta-1-\frac{d}{2}=O(d)$, we have

[TABLE]

• •

when $\beta>\frac{d}{2}+1$ and $\beta-1-\frac{d}{2}=O(1)$, we have

[TABLE]

3.2 Non-isotropic densities with quadratic-like $V$ outside of a ball

In this section, we estimate the expectations for a class of non-isotropic densities in the form of $\pi_{\beta}\propto V^{-\beta}$ with $V$ satisfying the following Lyapunov condition:

[TABLE]

The above Lyapunov condition characterizes the class of $V$ that are ‘quadratic-like’ outside a ball of radius $R$. If we assume that $V$ has Lipschitz gradients, then when $\beta$ is sufficiently large, the above assumption is satisfied if $V$ satisfies the PL inequality $|\nabla V(x)|^{2}\geq a^{2}V(x)$ wherever $|x|\geq R$ with some $a>0$ and it is from this inequality that quadratic growth follows. In particular, if $V$ satisfies the gradient Lipschitz assumption with parameter $L$, we have that for all $\beta\geq 1+a^{-2}(dL+\varepsilon)$,

[TABLE]

thereby leading to the Lyapunov condition in (15).

Proposition 3**.**

If $V\in\mathcal{C}^{2}(\mathbb{R}^{d})$ is positive, $L$-gradient Lipschitz and satisfies (15), then we have

[TABLE]

Proof.

Since $\mathcal{L}$ is ergodic with stationary distribution $\pi_{\beta}$, we have

[TABLE]

with $(X_{t})_{t\geq 0}$ being the solution to (3) with initial condition $X_{0}=x$. We will first bound $\mathbb{E}\left[V(X_{t})\right]$ and then take $t\to\infty$. Let $(P_{t})_{t\geq 0}$ be the Markov semigroup of $\eqref{eq:diffusion strongly convex V}$, then

[TABLE]

With (4), we have

[TABLE]

where the first inequality follows from (15) and the fact that $\Delta V\leq d\left\lVert\nabla^{2}V\right\rVert_{2}$. Therefore we obtain

[TABLE]

and it follows from Gronwall’s inequality that

[TABLE]

We hence have that $\mathbb{E}_{\pi_{\beta}}\left[V(X)\right]\leq\left(dL+\varepsilon\right)\max_{|x|\leq R}V(x)$ by taking $t\to\infty$. For the other expectation, we have

[TABLE]

The above identity implies

[TABLE]

∎

3.3 General Case

Next we discuss the general case where $\pi_{\beta}=Z_{\beta}^{-1}V^{\beta}$ and $V\in\mathcal{C}^{2}(\mathbb{R}^{d})$ is positive such that there exist constants $\alpha,L>0$ and $\alpha I_{d}\preceq\nabla^{2}V(x)\preceq LI_{d}$ for all $x\in\mathbb{R}^{d}$. Since $V$ is strongly convex, there is a unique $x^{*}\in\mathbb{R}^{d}$ such that $V(x)\geq V(x^{*})>0$ for all $x\in\mathbb{R}^{d}$ and $\nabla V(x^{*})=0$. Without loss of generality, we assume $x^{*}=0$.

Proposition 4**.**

Let $\beta>\frac{d}{2}+1$. If $V\in\mathcal{C}^{2}(\mathbb{R}^{d})$ is positive, $\alpha$-strongly convex and $L$-gradient Lipschitz, we have for any $r\in(0,\beta-\frac{d}{2}-1)$,

[TABLE]

Proof.

For any $r\in(0,\beta-\frac{d}{2}-1)$, we have

[TABLE]

where the last inequality follows from Lemma 3. For the other expectation, we have

[TABLE]

The above identity implies

[TABLE]

where the last inequality also follows from Lemma 3. ∎

Remark 5**.**

A ratio between Gamma functions appears in (17) and (18). The ratio can be written explicitly via properties of Gamma functions.

• •

When $d$ is an even number and $d=2k$ for some integer $k$,

[TABLE]

• •

When $d$ is an odd number with $d=2k-1$ for some integer $k$,

[TABLE]

where the first inequality follows from Gautschi’s inequality (Ismail and Muldoon, 1994).

Remark 6**.**

With Theorem 4 and the upper bounds in Remark 5, we can get the estimations for $\mathbb{E}_{\pi_{\beta}}\left[|\nabla V(X)|^{2}\right]$ and $\mathbb{E}_{\pi_{\beta}}\left[V(X)\right]$: for any $r\in(0,\beta-\frac{d}{2}-1)$,

[TABLE]

4 Zeroth-Order Itô Discretization

While previously we consider the case when the gradient of the function $V$ is analytically available to us, we now consider the case when we have access only to the function evaluations. This setting is called the zeroth-order setting and has been recently examined in the context of complexity of sampling in the works of Dwivedi et al. (2019), Lee et al. (2021), Roy et al. (2022). In this setting, we construct an approximation to the gradient via zeroth-order information, i.e., function evaluations. For simplicity, we consider the case of obtaining exact function evaluations. Based on the Gaussian smoothing technique (Nesterov and Spokoiny, 2017, Roy et al., 2022), for any $x\in\mathbb{R}^{d}$, we define the zeroth order gradient estimator $g_{\sigma,,m}(x)$ as

[TABLE]

where $u_{i}\sim\mathcal{N}(0,I_{d})$ are assumed to be independent and identically distributed. The parameter $m$ is called the batch size parameter. Then the zeroth order algorithm to sample $\pi_{\beta}$ is given by

[TABLE]

where $h>0$ is the step size and $\{\xi_{k+1}\}_{k=0}^{\infty}$ is a sequence of independent identically distributed standard Gaussian random vectors in $\mathbb{R}^{d}$. From Balasubramanian and Ghadimi (2022) and Roy et al. (2022), we recall the following property of $g_{\sigma,m}$.

Proposition 5**.**

(Roy et al., 2022, Section 8.1)* Assume $V$ is $L$-gradient Lipschitz. Define $\zeta_{k}=g_{\sigma,m}(x_{k})-\nabla V(x_{k})$ with $g_{\sigma,m}$ defined in (21) and $\{x_{k}\}_{k=0}^{\infty}$ generated by (22). We have for any $k\geq 0$,*

[TABLE]

and

[TABLE]

Theorem 3**.**

Suppose $V$ is gradient-Lipschitz with parameter $L>0$ and satisfies Assumption 2 with $\delta$ in (10). Let $g_{\sigma,m}$ be as defined in (21) and $(x_{k})_{k=0}^{\infty}$ be generated from (22) with $x_{k}\sim\nu_{k}$ for all $k\geq 0$. Then with the time step size

[TABLE]

the decay of Wasserstein-2 distance along the Markov chain $(x_{k})_{k=0}^{\infty}$ can be described by the following equation. For all $k\geq 1$,

[TABLE]

with $A^{\prime},B^{\prime}$ and $C^{\prime}$ given respectively in (60), (61) and (62).

Remark 7**.**

With Theorem 3, we can study the iteration complexity to reach an $\varepsilon$-accuracy in Wasserstein-2 distance. In the following discussion, we focus on the dimension dependence and $\varepsilon$ dependence in the iteration complexity. When $\beta=\Theta(d)$ and $\alpha,L=\Theta(1)$, and when $h$ satisfies (25), we have

[TABLE]

To ensure $W_{2}(\nu_{K},\pi_{\beta})<\varepsilon$, we require that each of

[TABLE]

is smaller than $\varepsilon/3$. Setting $\sigma=\varepsilon\delta/\sqrt{d}$, and

[TABLE]

we hence obtain that the iteration complexity $K$ is of order

[TABLE]

*The number of function evaluations is hence $mK$. *

5 Illustrative Examples

We now provide illustrative examples to highlight the implications of our results.

5.1 Multivariate $t$-distribution: Large Degree of Freedom

We first consider the isotropic multivariate $t$-distribution with the degrees of freedom being $d+2$. We choose $V(x)=1+|x|^{2}$, $\beta=d+1$ and $\pi_{\beta}(x)\propto V(x)^{-\beta}=(1+|x|^{2})^{-(d+1)}$. With this choice of $V$ and $\beta$, $V$ satisfies Assumption 2 with $\alpha=2,\ C_{V}=2$, and $V$ is $L$-Lipschitz gradient with $L=2$. The constant $\delta$ in Theorem 2 becomes $\delta=1$. Furthermore, according to proposition 2, $\mathbb{E}_{\pi_{\beta}}[V(X)]=2$ and $\mathbb{E}_{\pi_{\beta}}[|\nabla V(X)|^{2}]=4$.

5.1.1 First order algorithm

According to Theorem 2 and (2), to obtain $\epsilon$-accuracy in Wasserstein-2 distance, the iteration complexity is of order $\tilde{O}(1/\epsilon^{2})$. With the same choice of $V$ and $\beta$, we check the conditions of Theorem 1 in Li et al. (2019). The diffusion (3) is $\alpha^{\prime}$-uniformly dissipative with $\alpha^{\prime}=d$ and the Euler discretization given in (9) has local deviation with order $(p_{1},p_{2})=(1,{3}/{2})$ and $(\lambda_{1},\lambda_{2})=(\Theta(d^{5}),\Theta(d^{4}))$. The detailed calculation for deriving the constants above is provided in Appendix B. Hence, by Theorem 1 in Li et al. (2019), to reach an $\epsilon$-accuracy in Wasserstein-2 distance, the iteration complexity is of order $\tilde{O}(d^{3}/\epsilon^{2})$. Hence, in comparison with the result in Li et al. (2019), we obtain a dimension-free iteration complexity to ensure an $\epsilon$-accuracy in Wasserstein-2 distance.

5.1.2 Zeroth order algorithm

According to Theorem 3 and (27), to obtain $\varepsilon$-accuracy in Wasserstein-2 distance, the iteration complexity is of order $\tilde{O}\left((1\vee d/m)/\varepsilon^{2}\right)$. When $m=1$, the iteration complexity $K\sim\tilde{O}(d/\varepsilon^{2})$ and the number of functions evaluations $mK$ is also of the same order $\tilde{O}(d/\varepsilon^{2})$. If we choose the batch size $m=d$, we get a dimension independent iteration complexity $K\sim\tilde{O}(1/\varepsilon^{2})$ but the number of function evaluations is of order $\tilde{O}(d/\varepsilon^{2})$. Hence, we notice that in the case of multivariate $t$-distribution distributions with large degrees of freedom, the cost of estimating the gradient has an effect on the sampling complexities.

5.2 Multivariate $t$-distribution: Small Degrees of Freedom

We now consider the isotropic multivariate $t$-distribution with the degrees of freedom being $3$. We denote the corresponding density function by $\pi_{\beta}$. The exact number of 3 is chosen just for convenience; the results of this example apply to all cases where the degrees of freedom is strictly larger than 2 which corresponds to the setting where the variance is finite. We choose $V(x)=1+|x|^{2}$, $\beta={(d+3)}/{2}$ and $\pi_{\beta}(x)\propto V(x)^{-\beta}=(1+|x|^{2})^{-(d+3)/2}$. With the above choice of $V$ and $\beta$, $V$ satisfies Assumption 2 with $\alpha=2,\ C_{V}=2$ and $V$ is $L$-Lipschitz gradient with $L=2$. Hence, the constant $\delta$ in Theorem 2 is given by $\delta=1/d$. According to Proposition 2, $\mathbb{E}_{\pi_{\beta}}[V(X)]=d+1$ and $\mathbb{E}_{\pi_{\beta}}[|\nabla V(X)|^{2}]=4d$.

5.2.1 First order algorithm

According to Theorem 2 and (2), to obtain $\epsilon$-accuracy in Wasserstein-2 distance, the iteration complexity is of order $\tilde{O}(d^{4}/\epsilon^{2})$. With the same choice of $V$ and $\beta$, we check the conditions of Theorem 1 in Li et al. (2019). The diffusion (3) is $\alpha^{\prime}$-uniformly dissipative with $\alpha^{\prime}=1$ and the Euler discretization given in (9) has local deviation with order $(p_{1},p_{2})=(1,{3}/{2})$ and $(\lambda_{1},\lambda_{2})=(\Theta(d^{5}),\Theta(d^{4}))$. The detailed calculation for deriving the constants is provided in Appendix B. Hence, according to Theorem 1 in Li et al. (2019), to reach an $\epsilon$-accuracy in Wasserstein-2 distance, the iteration complexity is of order $\tilde{O}(d^{6}/\epsilon^{2})$. Even in this extremely heavy-tail case (i.e., only the variance exists), to ensure an $\epsilon$-accuracy in Wasserstein-2 distance, we can obtain an iteration complexity with polynomial dimension dependence. Furthermore, in comparison to Li et al. (2019), our analysis helps to decrease the dimension exponent by a factor of $2$.

5.2.2 Zeroth order algorithm

According to Theorem 3 and (27), to obtain $\varepsilon$-accuracy in Wasserstein-2 distance, the iteration complexity is of order $\tilde{O}\left(\max\{d^{4}/\varepsilon^{2},~{}d^{\frac{5}{2}}/\varepsilon,~{}d^{4}/\varepsilon^{2}m\}\right)$. Hence, we have that for any batch size $m$, the iteration complexity $K=\tilde{O}(d^{4}/\varepsilon^{2})$. Picking $m=1$, the number of function evaluations are of the same order, i.e., $mK=\tilde{O}(d^{4}/\varepsilon^{2})$.

Remark 8**.**

The example discussed in Section 5.2.2 highlights the following important observation: Choosing a large batch size does not improve the iteration complexity. To explain this, we understand both (9) and (22) as approximation to the continuous dynamics (3). For the first-order algorithm, the error of the approximation only comes from the Euler-Maruyama discretization. For the zeroth-order algorithm, the error of the approximation comes from both the Euler-Maruyama discretization and the zeroth-order gradient estimate=. When the error from the Euler-Maruyama discretization dominates, the optimal batch size is always $1$ and the oracle complexity of the zeroth order algorithm is the same as the iteration complexity for the first-order algorithm. When the error from the zeroth-order gradient estimate dominates, we need to choose a large batch size depending on $d$ so that the iteration complexity for the zeroth-order algorithm is the same as the iteration complexity for the first-order algorithm while the zeroth-order oracle complexity is of order $m$-times larger.

6 Further Results and Additional Insights on Assumptions

In Section 2, we provide sufficient conditions on $V$ such that when $\beta>d$, $\pi_{\beta}\propto V^{-\beta}$ satisfies the weighted Poincaré inequality with weight $V$. In this section, we relax the conditions in Section 2 by introducing the following assumptions.

Assumption 3**.**

The function $V:\mathbb{R}^{d}\to(0,\infty)$ is twice continuously differentiable and $V$ satisfies

• (1)

$\nabla^{2}V(x)$* is invertible for all $x\in\mathbb{R}^{d}$.*

• (2)

There exists $\gamma\in\left(0,\frac{\beta}{d+2}\right]$, such that

[TABLE]

where $V_{\gamma}:=V^{\gamma}$ and $C_{V}(\gamma)$ is a positive constant depending on $\gamma$.

Lemma 2**.**

Under Assumption 3, for any smooth function $\phi\in L^{2}(\pi_{\beta})$,

[TABLE]

Proof.

First we define $V_{\gamma}:=V^{\gamma}$. Choose $\beta^{\prime}=\beta-2\gamma$. For $\pi_{\beta^{\prime}}\propto V^{-\beta^{\prime}}$, we can write it as $\pi_{\beta^{\prime}}\propto{V_{\gamma}}^{-a}$ with

[TABLE]

where the inequality follows from the fact that $\gamma\in\left(0,\frac{\beta}{d+2}\right]$. Therefore we can apply Theorem 1 to $\pi_{\beta^{\prime}}\propto{V_{\gamma}}^{-a}$ and get for any smooth, $\pi_{\beta^{\prime}}$-square integrable function $g$ with $\mathbb{E}_{\pi_{\beta^{\prime}}}[g(X)]=0$ and $G=V_{\gamma}g$,

[TABLE]

Since $\beta^{\prime}=\beta-2\gamma$, (29) is equivalent to

[TABLE]

Under Assumption 3, we have

[TABLE]

where the last identity follows from the fact that $\beta^{\prime}=\beta-2\gamma$. Along with (30), we get

[TABLE]

Since $G=V^{\gamma}g$, $G$ is smooth, $\pi_{\beta}$-square integrable and $\mathbb{E}_{\pi_{\beta-\gamma}}[G(X)]=0$. For any $\pi_{\beta}$-square integrable $\phi$, let $G=\phi-\mathbb{E}_{\pi_{\beta-\gamma}}[\phi(X)]$ and we get

[TABLE]

Therefore for any smooth, $\pi_{\beta}$-square integrable $\phi$,

[TABLE]

which is equivalent to (28) with $C_{\text{WPI}}=\frac{C_{V}(\gamma)}{a+1}=C_{V}(\gamma)\left(\frac{\beta}{\gamma}-1\right)^{-1}$. ∎

Remark 9**.**

Lemma 2 can be applied to the class of multivariate $t$-distributions with $V(x)=1+|x|^{2}$. When $\beta\in\left(\frac{d+2}{2},d\right]$, with the choice of $\gamma=\frac{\beta}{d+2}$, Assumption 3 holds with

[TABLE]

Hence, Lemma 2 implies that the multivariate $t$-distribution with degree of freedom $\nu\in(2,d]$ satisfies the weighted Poincaré inequality with weight $1+|x|^{2}$ and with

[TABLE]

The detailed calculation for deriving the above mentioned constants is provided in Appendix C.

As an immediate consequence of Lemma 2, we have the following $\chi^{2}$ convergence result for (3).

Proposition 6**.**

Under Assumption 3, with $(X_{t})$ satisfying (3) with $\rho_{t}$ being the distribution of $X_{t}$, we have

[TABLE]

For the case of multivariate $t$-distributions, Proposition 6 allows us to show exponential convergence of (3) in the $\chi^{2}$ divergence with smaller degrees of freedom (and hence heavier tails) compared to Proposition 1.

6.1 Relationship between Lemma 1 and Lemma 2

The result in Lemma 2 complements that in Lemma 1. It can be used to study the WPI for $\pi_{\beta}$ when $\beta\leq d$. In particular, when $\beta\leq d$, if $\pi_{\beta}\propto V^{-\beta}$ and $V$ satisfies Assumption 2 with $C_{V}\in(0,\frac{d+2}{d+2-\beta})$, then V satisfies Assumption 3. Therefore $\pi_{\beta}$ satisfies the WPI. In Proposition 7, this relation is proved formally.

Proposition 7**.**

When $\beta\leq d$, if Assumption 2 holds with $C_{V}\in(0,\frac{d+2}{d+2-\beta})$, then Assumption 3 holds.

Proof.

First $\nabla^{2}V$ is invertible because $\nabla^{2}V\succeq\alpha I_{d}$. Next we show that there exists $\gamma\in(0,\frac{\beta}{d+2}]$ such that $\left\lVert V(x)^{\gamma-1}(\nabla^{2}V_{\gamma})^{-1}(x)\right\rVert_{2}\leq C_{V}(\gamma)$ for all $x\in\mathbb{R}^{d}$. It is equivalent to showing that there exists $\gamma\in(0,\frac{\beta}{d+2}]$ such that $\left\lVert V(x)^{1-\gamma}(\nabla^{2}V_{\gamma})(x)\right\rVert_{2}>0$ for all $x\in\mathbb{R}^{d}$. From the calculations in Section C, we have

[TABLE]

Therefore

[TABLE]

where the inequality follows from Assumption 2. Last we show that there exists $\gamma\in(0,\frac{\beta}{d+2}]$ such that $1-(1-\gamma)C_{V}>0$. Note that

[TABLE]

Since $C_{V}\in\left(0,\frac{d+2}{d+2-\beta}\right)$, we have that

[TABLE]

Therefore there exists a constant $\gamma\in\left(0,\frac{\beta}{d+2}\right]$ such that $\left\lVert V(x)^{1-\gamma}(\nabla^{2}V_{\gamma})(x)\right\rVert_{2}>0$ for all $x\in\mathbb{R}^{d}$. ∎

6.2 Relationship between Theorem 2 and Proposition 6

Proposition 6 studies the convergence of the continuous dynamics (3) while Theorem 2 studies the convergence of the discretization (9). The conditions in Theorem 2 can be shown to imply conditions in proposition 6. In Proposition 6 we only assume Assumption 3. In Theorem 2, we assume (i) Assumption 2, (ii) $\delta=\frac{\beta-1-\frac{1}{4}C_{V}d}{\frac{1}{4}C_{V}d}>0$, and (iii) $V$ is gradient Lipschitz. In the following proposition, we show that these three assumptions together imply Assumption 3.

Proposition 8**.**

If Assumption 2 holds such that $\delta=\frac{\beta-1-\frac{1}{4}C_{V}d}{\frac{1}{4}C_{V}d}>0$ and $V$ is $L$-gradient Lipschitz, then Assumption 3 holds.

Proof of Proposition 6.

Under Assumption 2 and $L$-gradient Lipschitzness assumption, we have that $V$ is ‘essential quadratic’. That is, assuming $V$ attains its global minimum at $x^{*}$, for all $x\in\mathbb{R}^{d}$,

[TABLE]

Therefore for all $x\in\mathbb{R}^{d}$,

[TABLE]

which implies that Assumption 2-(2) is satisfied with $C_{V}=\frac{2L^{2}}{\alpha^{2}}$. Furthermore,

[TABLE]

The condition $\delta=\frac{\beta-1-\frac{1}{4}C_{V}d}{\frac{1}{4}C_{V}d}>0$ is equivalent to the condition $\beta>\frac{L^{2}}{2\alpha^{2}}d+1$. Notice that for all $d\geq 1$, we have

[TABLE]

Therefore for any

[TABLE]

we can choose $\gamma=\frac{\beta}{d+2}$ and obtain

[TABLE]

Therefore Assumption 3-(2) is satisfied with $\gamma=\beta/(d+2)$ and

[TABLE]

The proof is now complete because Assumption 3-(1) is automatically satisfied under Assumption 2.

∎

7 Proofs of the Main Results

7.1 Proofs of Theorem 2 and Theorem 3

In this section, we provide the proof of Theorem 2 and Theorem 3 via mean square analysis. We first start with the following intermediate result.

Proposition 9**.**

Let $(X_{t})_{t\geq 0}$ follow (3) with $X_{t}\sim\rho_{t}$ for all $t\geq 0$. If $V$ is gradient Lipschitz with parameter $L$, then we have

[TABLE]

Proof of Proposition 9.

According to (3),

[TABLE]

where

[TABLE]

and

[TABLE]

With (7.1) and (36), we get

[TABLE]

By Gronwall’s inequality, we hence have

[TABLE]

∎

Based on the above proposition, we now prove Theorem 2 below.

Proof of theorem 2.

We perform mean square analysis to (9). Let $(X_{t})_{t\geq 0}$ follow (3) with $X_{0}\sim\pi_{\beta}$. Since $\pi_{\beta}$ is the unique stationary distribution to (3), $X_{t}\sim\pi_{\beta}$ for all $t\geq 0$. With (9), we can calculate the difference between $X_{h}$ and $x_{1}$,

[TABLE]

where

[TABLE]

Therefore according to triangle inequality,

[TABLE]

Since $U_{1}$ is adapted to $\mathcal{F}_{0}$ and $\mathbb{E}[U_{3}|\mathcal{F}_{0}]=0$, we get

[TABLE]

Since $V$ is $\alpha$-strongly convex and $L$-gradient Lipschitz, it satisfies

[TABLE]

Therefore when $h\leq\frac{2}{(\beta-1)(\alpha+L)}$,

[TABLE]

Meanwhile, for arbitrary $r>0$, we have

[TABLE]

Notice that under Assumption 2, we have

[TABLE]

for all $x\in\mathbb{R}^{d}$. Therefore

[TABLE]

and

[TABLE]

With (42) and (43), we get

[TABLE]

Next we apply Proposition 9 to $\mathbb{E}[|X_{t}-X_{0}|^{2}]$. In particular, when

[TABLE]

we have

[TABLE]

Combining (44) and (7.1), when $h<\frac{1}{4(\beta-1)L}$, we have that

[TABLE]

With (7.1) and (7.1), we get

[TABLE]

Since $C_{V}<\frac{4(\beta-1)}{d}$, denote $\delta=\frac{(\beta-1)-\frac{1}{4}C_{V}d}{\frac{1}{4}C_{V}d}>0$. We have

[TABLE]

By picking $r=\frac{\delta}{3}$, we get for any $h\in\left(0,\frac{2\delta}{3(1+\delta)\alpha(\beta-1)}\right)$ that

[TABLE]

With the choice of $r=\delta/3$, (7.1) could be rewritten as

[TABLE]

Next, with the bound in (7.1), we get when $h<\frac{1}{4(\beta-1)L}$,

[TABLE]

With (7.1) and (7.1), we get when $h<\min\left(\frac{1}{4(\beta-1)L},\frac{2\delta}{3(1+\delta)\alpha(\beta-1)}\right)$,

[TABLE]

with

[TABLE]

The above analysis works for each step, therefore we get for all $k\geq 1$,

[TABLE]

According to (Dalalyan and Karagulyan, 2019, Lemma 9), with $A,B,C$ given in (51),(52),(53), for all $k\geq 1$,

[TABLE]

Choosing $X_{0}$ such that $W_{2}(\nu_{0},\pi_{\beta})=\mathbb{E}\left[|X_{0}-x_{0}|^{2}\right]^{\frac{1}{2}}$, we get (11).

∎

We now prove Theorem 3.

Proof of Theorem 3.

Following the same strategy and notation in the proof of Theorem 2, we have

[TABLE]

where $U_{1},U_{2},U_{3}$ are defined in (38),(39),(40) respectively and $\zeta_{0}=g_{\sigma,m}(x_{0})-\nabla V(x_{0})$. Therefore we have

[TABLE]

From the proof of Theorem 2 and Proposition 5, when

[TABLE]

we have that

[TABLE]

where $A,B,C$ are defined in (51),(52),(53). Using the fact that $V$ is gradient Lipshcitz, we have

[TABLE]

Plugging (7.1) in (7.1), we get

[TABLE]

When we pick the step-size such that

[TABLE]

we have

[TABLE]

Therefore we have

[TABLE]

where

[TABLE]

The rest of the proof is the same as the proof of Theorem 2, and hence we get (26). ∎

Acknowledgements

YH was supported in part by NSF TRIPODS grant CCF-1934568. TF was supported by the Engineering and Physical Sciences Research Council (EP/T5178) and by the DeepMind scholarship. KB was supported in part by NSF grant DMS-2053918. MAE was supported by NSERC Grant [2019-06167], Connaught New Researcher Award, CIFAR AI Chairs program, and CIFAR AI Catalyst grant. Parts of this work was done when YH, KB and MAE visited the Simons Institute for the Theory of Computing as a part of the “Geometric Methods in Optimization and Sampling” program during Fall 2021.

Appendix A Computations for Section 3.3

Lemma 3**.**

Let $\beta>\frac{d}{2}+1$. If $V\in\mathcal{C}^{2}(\mathbb{R}^{d})$ is positive, $\alpha$-strongly convex and $L$-gradient Lipschitz, we have for any $r\in(0,\beta-\frac{d}{2}-1)$,

[TABLE]

Proof.

Since $V(x)\leq V(0)+\frac{L}{2}|x|^{2}$, we know that for any $r\in(0,\beta-\frac{d}{2}-1)$, $Z_{\frac{d}{2}+r}$ is finite and $\pi_{\frac{d}{2}+r}$ is a probability measure. Therefore

[TABLE]

For the integral $\int_{\mathbb{R}^{d}}(V(0)+\frac{L}{2}|x|^{2})^{-\beta}dx$, we can calculate it via change of polar coordinates and substitutions,

[TABLE]

where the second identity follows from a substitution with $R_{L}=LR^{2}/(2V(0))$ and the fourth identity follows from a substitution with $u=\frac{R_{L}}{1+R_{L}}$. For $Z_{\frac{d}{2}+r}$, we have

[TABLE]

Therefore, we can further get

[TABLE]

∎

Appendix B Computations for Sections 5.1 and 5.2

Let $\pi_{\beta}(x)\propto V(x)^{-\beta}=(1+|x|^{2})^{-\beta}$ with $\beta>\frac{d+2}{2}$. The gradient and Hessian of $V$ is

[TABLE]

Therefore $V$ is $\alpha$-strongly convex with $\alpha=2$ and $L$-gradient Lipschitz with $L=2$. (3) reduces to

[TABLE]

with $b(x)=-2(\beta-1)x$ and $\sigma(x)=\sqrt{2}(1+|x|^{2})^{\frac{1}{2}}I_{d}$.

Next we look at the uniform dissipativity condition:

[TABLE]

where the inequality follows from the fact that $x\mapsto(1+|x|^{2})^{\frac{1}{2}}$ is $1$-Lipschitz. Therefore diffusion (64) is $\alpha^{\prime}$-uniform dissipative with $\alpha^{\prime}=2(\beta-1-\frac{d}{2})$. In particular, $\alpha^{\prime}=d$ when $\beta=d+1$ and $\alpha^{\prime}=1$ when $\beta=\frac{d+3}{2}$.

Last we look at the local deviation for the Euler discretization to (64). We use the same notations in Li et al. (2019). According to (Li et al., 2019, lemma 29), $p_{1}=1$ and

[TABLE]

According to (Li et al., 2019, lemma 29), $p_{2}=\frac{3}{2}$ and

[TABLE]

with

[TABLE]

The order of $\lambda_{1}$ and $\lambda_{2}$ in dimension parameter $d$ is given by:

[TABLE]

Therefore, we have that

• •

when $\beta=d+1$, $(\lambda_{1},\lambda_{2})=(\Theta(d^{5}),\Theta(d^{4}))$,

• •

when $\beta=\frac{d+3}{2}$, $(\lambda_{1},\lambda_{2})=(\Theta(d^{5}),\Theta(d^{4}))$.

Appendix C Computations for Remark 9

In the example of Cauchy class distributions, $V(x)=1+|x|^{2}$ and $V_{\gamma}:=V^{\gamma}$. When $\gamma>\frac{1}{2}$,

[TABLE]

Plug in $V(x)=1+|x|^{2}$, we get

[TABLE]

and

[TABLE]

When $\beta\in\left(\frac{d+2}{2},d\right]$, $\gamma=\frac{\beta}{d+2}\in\left(\frac{1}{2},1\right]$,

[TABLE]

Therefore Assumption 3 holds with $C_{V}(\gamma)=\frac{(d+2)^{2}}{2\beta(2\beta-d-2)}$. For the Cauchy distribution $\pi_{\beta}\propto(1+|x|^{2})^{-\beta}=(1+|x|^{2})^{-\frac{d+\nu}{2}}$ with $\beta\in(\frac{d+2}{2},d]$, i.e. $\nu\in(2,d]$, according to lemma 2, $\pi_{\beta}$ satisfies the weighted Poincaré inequality with weight $1+|x|^{2}$ with weighted Poincaré constant

[TABLE]