Stochastic Zeroth-order Discretizations of Langevin Diffusions for Bayesian Inference

TL;DR

This paper introduces stochastic zeroth-order discretizations for Langevin diffusions, enabling Bayesian sampling without gradient evaluations, and provides theoretical analysis of their sample complexity and effectiveness in noisy, high-dimensional scenarios.

Contribution

The work develops and analyzes zeroth-order Langevin discretizations using Gaussian Stein's identities, extending sampling methods to black-box and high-dimensional contexts.

Findings

Provides sample complexity bounds for zeroth-order Langevin algorithms.

Introduces a variable selection technique with theoretical guarantees.

Extends Langevin sampling to noisy black-box and high-dimensional settings.

Abstract

Discretizations of Langevin diffusions provide a powerful method for sampling and Bayesian inference. However, such discretizations require evaluation of the gradient of the potential function. In several real-world scenarios, obtaining gradient evaluations might either be computationally expensive, or simply impossible. In this work, we propose and analyze stochastic zeroth-order sampling algorithms for discretizing overdamped and underdamped Langevin diffusions. Our approach is based on estimating the gradients, based on Gaussian Stein's identities, widely used in the stochastic optimization literature. We provide a comprehensive sample complexity analysis -- number noisy function evaluations to be made to obtain an $\epsilon$-approximate sample in Wasserstein distance -- of stochastic zeroth-order discretizations of both overdamped and underdamped Langevin diffusions, under various noise models. We also propose a variable selection technique based on zeroth-order gradient estimates and establish its theoretical guarantees. Our theoretical contributions extend the practical applicability of sampling algorithms to the noisy black-box and high-dimensional settings.

Equations

$$\pi (\theta )=\frac{e^{-f(\theta )}}{{\int }_{{\mathbb{R}}^d}{e}^{-f(r)}dr}$$

$$d{L}_T=-\nabla f(L_T)dT+\sqrt{2}d{W}_T,$$

$$x_{n+1}=x_n-h_{n+1}\nabla f(x_n)+\sqrt{2h_{n+1}}{\epsilon }_{n+1},$$

$$W_2(p,q):={\left(\underset{\varrho \in \varrho (p,q)}{\inf }{\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}\Vert \theta -{\theta }^{\prime }{\Vert }_2^{2}d\varrho (\theta ,{\theta }^{\prime })\right)}^{1/2},$$

$$\mathbf{E}[\Vert \mathbf{E}({\zeta }_n|x_n){\Vert }^2]\le {\delta }_b^{2}d\text{and}\mathbf{E}[\Vert {\zeta }_n-\mathbf{E}({\zeta }_n|x_n){\Vert }^2]\le {\delta }_v^{2}d.$$

$$W_2({\varpi }_n,\pi )\le (1-mh{)}^n{W}_2({\varpi }_0,\pi )+1.65\frac{M}{m}(hd{)}^{1/2}+\frac{{\delta }_b\sqrt{d}}{m}+\frac{{\delta }_v^2(hd{)}^{1/2}}{1.65M+\sigma \sqrt{m}}.$$

$$\mathbf{E}[\Vert \mathbf{E}({\zeta }_n|x_n){\Vert }^2]\le {\delta }_b^2{d}^{\alpha }\text{and}\mathbf{E}[\Vert {\zeta }_n-\mathbf{E}({\zeta }_n|x_n){\Vert }^2]\le {\delta }_v^2{d}^{\beta },$$

$$W_2({\varpi }_n,\pi )\le$$

$$+\frac{{\delta }_v^{2}h{d}^{\beta }}{1.65M(hd{)}^{1/2}+{\delta }_b{d}^{\alpha /2}+{\delta }_v(mh{)}^{1/2}{d}^{\beta /2}}.$$

$$g_{\nu ,b}(\theta )=\frac{1}{b}\sum _{i=1}^b\frac{F(\theta +\nu u_i,{\xi }_i)-F(\theta ,{\xi }_i)}{\nu }{u}_i$$

$$\mathbf{E}\left[\Vert g_{\nu ,b}(\theta )-\nabla f_{\nu }(\theta ){\Vert }^2\right]\le \frac{2(d+5)(\Vert \nabla f(\theta ){\Vert }^2+{\sigma }^2)}{b}+\frac{{\nu }^2{M}^2(d+3{)}^3}{2b}$$

$$\mathbf{E}\left[\Vert g_{\nu ,b}(\theta )-\nabla f(\theta ){\Vert }^2\right]\le \frac{4(d+5)(\Vert \nabla f(\theta ){\Vert }^2+{\sigma }^2)}{b}+\frac{3{\nu }^2{M}^2(d+3{)}^3}{2}$$

$$x_{n+1}=x_n-h{g}_{\nu ,b}(x_n)+\sqrt{2h}{\epsilon }_{n+1}$$

$$h=\frac{{ϵ}^2}{d^2},b=\max (1,{\sigma }^2)d,\nu =\frac{ϵ}{\sqrt{d}},$$

$$N=\mathcal{O}\left(\frac{d}{{ϵ}^2}\cdot \log \left(\frac{d}{ϵ}\right)\right).$$

$$Nb=\mathcal{O}\left(\frac{\max (1,{\sigma }^2)\cdot d^2}{{ϵ}^2}\cdot \log \left(\frac{d}{ϵ}\right)\right).$$

$$d{V}_T$$

$$d{L}_T$$

$${\tilde{x}}_{n+1}$$

$$x_{n+1}$$

$${\tilde{x}}_{n+1}$$

$$x_{n+1}$$

$$h=\frac{mϵ}{12\gamma M\sqrt{d}},\nu =\frac{ϵ}{\sqrt{d}},b=\frac{d^{1.5}\max (1,{\sigma }^2)}{ϵ},$$

$$N=\tilde{\mathcal{O}}\left(\frac{\sqrt{d}}{ϵ}\right)$$

$$Nb=\tilde{\mathcal{O}}\left(\frac{d^2\max (1,{\sigma }^2)}{{ϵ}^2}\right).$$

$$x_{n+\frac{1}{2}}=x_n+\frac{1-e^{-2{\alpha }_{n}h}}{2}{v}_n-\frac{u}{2}\left({\alpha }_{n}h-\frac{1-e^{-2({\alpha }_{n}h)}}{2}\right)g_{\nu ,b}(x_n)+\sqrt{u}{ϵ}_{n+1}^{(1)}$$

$$x_{n+1}=x_n+\frac{1-e^{-2h}}{2}{v}_n-\frac{uh}{2}(1-e^{-2(h-{\alpha }_{n}h)})g_{\nu ,b}\left(x_{n+\frac{1}{2}}\right)+\sqrt{u}{ϵ}_{n+1}^{(2)}$$

$$v_{n+1}=v_n{e}^{-2h}-uh{e}^{-2(h-{\alpha }_{n}h)}{g}_{\nu ,b}\left(x_{n+\frac{1}{2}}\right)+2\sqrt{u}{ϵ}_{n+1}^{(3)}.$$

$$h=C\min \left(\frac{(ϵ\sqrt{m}{)}^{\frac{1}{3}}}{(d\kappa {)}^{\frac{1}{6}}\log {\left(\frac{1}{ϵ}\right)}^{\frac{1}{6}}},\min \left({\left(\frac{m}{d}\right)}^{\frac{1}{3}},{\left(\frac{Mm}{16{\sigma }^2}\right)}^{\frac{1}{3}},\sqrt{m}\right){ϵ}^{\frac{2}{3}}\log {\left(\frac{1}{ϵ}\right)}^{-\frac{2}{3}}\right)b=\frac{d\kappa }{h^3}\nu =\frac{u{h}^2}{d^{1.5}}$$

$$N=\tilde{O}\left(\max \left(\frac{d^{\frac{1}{6}}{\kappa }^{\frac{7}{6}}}{(ϵ\sqrt{m}{)}^{\frac{1}{3}}},\frac{\kappa \max \left({\left(\frac{d}{m}\right)}^{\frac{1}{3}},{\left(\frac{{\sigma }^2}{Mm}\right)}^{\frac{1}{3}},\frac{1}{\sqrt{m}}\right)}{{ϵ}^{\frac{2}{3}}}\right)\right)$$

$$2Nb=\tilde{O}\left(\max \left(\frac{d^{\frac{5}{3}}{\kappa }^{\frac{8}{3}}}{{ϵ}^{\frac{4}{3}}},\frac{d{\kappa }^2\max {\left({\left(\frac{d}{m}\right)}^{\frac{1}{3}},{\left(\frac{{\sigma }^2}{Mm}\right)}^{\frac{1}{3}},\frac{1}{\sqrt{m}}\right)}^4}{{ϵ}^{\frac{8}{3}}}\right)\right).$$

$${\int }_{{\mathbb{R}}^d}{g}^2(\theta )\log g^2(\theta )\pi (\theta )d\theta -{\int }_{{\mathbb{R}}^d}{g}^2(\theta )\pi (\theta )d\theta \log {\int }_{{\mathbb{R}}^d}{g}^2(\theta )\pi (\theta )d\theta \le \frac{2}{\lambda }{\int }_{{\mathbb{R}}^d}\Vert \nabla g(\theta ){\Vert }^2\pi (\theta )d\theta .$$

$$H_{\pi }(\varpi )={\int }_{{\mathbb{R}}^d}\varpi (\theta )\log \frac{\varpi (\theta )}{\pi (\theta )}d\theta J_{\pi }(\varpi )={\int }_{{\mathbb{R}}^d}\varpi (\theta ){\Vert \nabla \log \frac{\varpi (\theta )}{\pi (\theta )}\Vert }^{2}d\theta .$$

$$H_{\pi }(\varpi )\le \frac{1}{2\lambda }{J}_{\pi }(\varpi ).$$

Full text

1 Introduction

First generation sampling algorithms, for example, Metropolis-Hastings algorithm are oblivious to the geometry of the target density as a result of which they suffer from slower rates of convergence. However, they are efficiently implementable and widely applicable, as they are based only on exact density function evaluations; see, for example, [19, 43, 54, 50, 51, 27, 52], for more details about such algorithms. Motivated by statistical physics principles, various researchers developed second-generation of sampling algorithms, that leverage geometric information regarding the target density [59, 55, 60, 66, 67, 34, 6]. Such algorithms are based on gradient-based discretizations of continuous-time underdamped or overdamped Langevin diffusions. Although such algorithms were developed much earlier, recently strong theoretical guarantees have been established for sampling in the works of [30, 31, 24, 29, 21, 20, 25] and several others. Such algorithms typically perform empirically better and exhibit faster rates of convergence compared to the first generation sampling algorithms mentioned above.

In this work, given a density function $\pi:\mathbb{R}^{d}\to\mathbb{R}$, with potential function $f:\mathbb{R}^{d}\to\mathbb{R}$, of the form

[TABLE]

we consider the problem of sampling when we only have access to noisy evaluations of the potential function $f$. We refer to this problem as stochastic zeroth-order sampling. Our approach is based on discretizing overdamped and underdamped Langevin diffusions using stochastic zeroth-order oracles, which, when queried returns noisy unbiased evaluations, $F(x,\xi)$, of the function value $f(x)$. That is, we have ${\bf E}[F(x,\xi)]=f(x)$, where $\xi$ is the random noise in our function evaluations, which is not necessarily an additive noise. Our motivations for studying such problems are three-fold:

• •

Computationally Complexity of Gradient-evaluation: A majority of existing discretizations of Langevin diffusions require computing the gradient of the potential function $f$ in each iteration. It is well-know that for a wide class of functions which could be expressed based on compositions of elementary differentiable functions, the computational cost of evaluating the gradient is $4$ to $5$ times more than that of evaluating the function; see, for example [38]. Furthermore, in order to compute the gradient, it is necessary to store several intermediate gradients, which increases the memory requirement. Hence, for several potential functions, Langevin-discretization based sampling algorithms might end up spending more time and memory for computing and storing gradients in each iteration. To reduce the wall-clock runtimes of such sampling algorithms, it is of interest to develop discretization of Langevin diffusions based only on function evaluations.

• •

Non-availability of Analytic form of Potential Function: In a variety of scientific problems, the potential function $f$ might not even be available in closed form, either due to the sheer size of the dataset (see, for example [69]), or due to the constraints in the physical process underlying the statistical model (see, for example [8, 39, 42]). In these situations, we do not have access to the analytical from of true potential function, let alone its gradients, which are required for discretizing Langevin diffusions. Hence, it is of great interest to develop discretization of Langevin diffusions based on noisy function evaluations to widen the applicability of Bayesian inference. It is worth mentioning here that, in the case of Metropolis-Hastings algorithms, [5, 69] developed and analyzed the so-called Pseudo-Marginal Metropolis-Hasting algorithms which work with unbiased noisy density evaluations. However, similar algorithms for sampling based on discretizing Langevin diffusions are lacking in the literature, except for the recent work on [1] which considered a pseudo-marginal Hamiltonian Monte Carlo algorithm. Our second motivation for this work is to fill this gap and to develop and analyze a unified framework for stochastic zeroth-order discretization of Langevin diffusions for Bayesian inference.

• •

Automating Bayesian Inference: From a practitioner’s perspective, statistical modeling is an inherently iterative process. The probabilistic model is typically refined during the scientific process based on the fit to the data. In the context of sampling, this process could be understood as changing the potential function $f$ in the modeling process. However, each time the function $f$ is changed, it is also invariably required to re-code the sampling algorithm based on the analytically computed gradient of the function $f$ under consideration. Our third motivation in this work is to automate this process, to help the practitioner with quick experimentation. As we will see later, our proposed methodology allows for sampling from a wide variety of density functions in a unified manner, as long as we have an oracle to obtain (noisy) evaluations of the potential function $f$. It is worth mentioning that recently [58, 61, 45] developed related automated Bayesian inference algorithms based on variational inference.

1.1 Preliminaries

Consider the continuous-time Langevin diffusion process $\{L_{T}:T\in\mathbb{R}_{+}\}$ given by the following stochastic differential equation,

[TABLE]

where $T\in\mathbb{R}_{+}$ and $\{W_{T}:T\in\mathbb{R}_{+}\}$ is a $d$-dimensional Brownian motion and $\nabla f(\theta)\in\mathbb{R}^{d}$ denotes the gradient of $f(\theta)$. The Euler-Maruyama discretization of the process in (2) is given by the following Markov chain:

[TABLE]

for the discrete time index $n=0,1,2\ldots$. Here $\varepsilon_{n}\in\mathbb{R}^{d}$ is a sequence of independent standard Gaussian vectors, $h_{n}$ denotes the step-size and an initial point $x_{0}$ is assumed to be given. The above discretization is called as the Langevin Monte Carlo (LMC) sampling algorithm. The update step of the LMC sampling algorithm shares similarity with the standard gradient descent algorithm from the optimization literature. As a prelude to the rest of the paper, our main idea in this work is to provide a non-asymptotic analysis of using stochastic zeroth-order gradient estimators (described in details in Section 1.2) in place of the true gradient in (3) and related discretizations.

Denoting the distribution of the random vector $x_{n}$ by $\varpi_{n}$, to evaluate the performance of the sampling algorithm, the 2-Wasserstein distance between $\varpi_{n}$ and the target density $\pi(\theta)$ is considered. For measures, $p$ and $q$ defined on $(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$, the 2-Wasserstein distance is defined as:

[TABLE]

where $\varrho(p,q)$ is the set of joint distribution that has $p$ and $q$ as its marginals. The performance of the sampling updates is measured by the above 2-Wasserstein distance between the distribution $\varpi_{n}$ and the target density $\pi$, i.e., $W_{2}(\varpi_{n},\pi)$. Specifically, the iteration complexity of the algorithm is defined as the number of iterations $N$, required to get $W_{2}(\varpi_{N},\pi)\leq\epsilon$. We also define the notion of oracle complexity which is the number of calls to the first-order or stochastic zeroth-order oracle used to obtain $W_{2}(\varpi_{N},\pi)\leq\epsilon$. For the LMC algorithm in (3), as we use only one gradient evaluation in each iteration, the oracle and iteration complexity becomes the same.

In order to obtain theoretical guarantees, a common assumption made in the literature on LMC is that the function $f$ is smooth and strongly convex.

Assumption 1.1

Letting $\|\cdot\|=\|\cdot\|_{2}$ denote the Euclidean norm on $\mathbb{R}^{d}$, the potential function $f$,

1. A1:

is strongly convex i.e., $f(\theta)-f(\theta^{\prime})-\nabla f(\theta^{\prime})^{\top}(\theta-\theta^{\prime})\geq\frac{m}{2}\|\theta-\theta^{\prime}\|^{2}$, for $m>0$. 2. A2:

has Lipschitz continuous gradient, i.e., $\|\nabla f(\theta)-\nabla f(\theta^{\prime})\|\leq M\|\theta-\theta^{\prime}\|$ for $M>0$.

The above assumptions on the potential function in-turn makes the density function $\pi$ strongly log-concave and smooth. Such an assumption is satisfied in several sampling and Bayesian inference problems including sampling from mixture of Gaussian distributions and Bayesian logistic regression. Further assuming access to certain inaccurate gradients, [29] provide theoretical guarantees for sampling under Assumption 1.1. Specifically, instead of the true gradient $\nabla f(x_{n})$ in each step, it is assumed that we observe $g_{n}=g(x_{n})=\nabla f(x_{n})+\zeta_{n}$, for a sequence of random noise vectors $\zeta_{n}$ that satisfies certain bias and variance assumption. Then, the noisy LMC updates corresponds to the case of the updates in Equation 3, with $\nabla f(x_{n})$ replaced by $g_{n}$. For such an update, [29] have the following non-asymptotic result. Before providing the result, we remark that due to the assumptions on the stochastic gradient made in (5), this setting is referred to as stochastic first-order setting.

Theorem 1.2

[29]** Assume that the bias and variance of $\zeta_{n}$ satisfies respectively, for all $n=1,2,\ldots$,

[TABLE]

Let the function $f$ satisfy Assumption 1.1. If $h\leq 2/(m+M)$, the following result holds true.

[TABLE]

Remark 1

More generally, if the bounded bias and variance condition are changed to

[TABLE]

respectively, for some $\alpha,\beta>0$, the conclusion turns into

[TABLE]

Furthermore, in the case that $\beta>\max\{1,\alpha\}$, the last term is dominated by $d^{\beta/2}$.

1.2 The Zeroth-order Methodology

The use of zeroth-order information (i.e., noisy function evaluations) for optimizing a function goes back to the works of [46, 16], that used stochastic version of finite-difference gradient approximation methods for estimating the maximum of a regression function (or equivalently mode of a density function). Since, then zeroth-order optimization has developed into an independent field in itself; see, for example [64, 23, 3, 49] for an more up-to-date account of this field. More recently, the focus has been more on developing a non-asymptotic understanding of stochastic zeroth-order optimization [36, 28, 56, 11]. Despite the fact that stochastic zeroth-order optimization is a well-developed field, to the best of our knowledge, there is no prior work on using related techniques for the closely related problem of zeroth-order discretizations of Langevin diffusions; specifically in terms of non-asymptotic analysis.

We now describe the precise assumption made on the stochastic zeroth-order oracle in the first part of this work.

Assumption 1.3

For any $\theta\in\mathbb{R}^{d}$, the stochastic zeroth-order oracle outputs an estimator $F(\theta,\xi)$ of $f(\theta)$ such that, $\mathbf{E}\left[F(\theta,\xi)\right]=f(\theta)$, $\mathbf{E}\left[\nabla F(\theta,\xi)\right]=\nabla f(\theta)$, and $\mathbf{E}\left[\|\nabla F(\theta,\xi)-\nabla f(\theta)\|^{2}\right]\leq\sigma^{2}$.

The assumption above assumes that we have accesses to a stochastic zeroth-order oracle which provides unbiased function evaluations with bounded variance. It is worth noting that in the above, we do not necessarily assume the noise $\xi$ is additive. Our gradient estimator is then constructed by leverage the Gaussian smoothing technique [56, 36, 11], which is amenable for fine-grained non-asymptotic analysis. Specifically, for a point $\theta\in\mathbb{R}^{d}$, we define an estimate $g_{\nu,b}(\theta)$, of the gradient $\nabla f(\theta)$ as follows:

[TABLE]

where $u_{i}\sim N(0,I_{d})$ and are assumed to be independent and identically distributed. An interpretation of the gradient estimator in (6) as a consequence of Gaussian Stein’s identity, popular in the statistics literature [68], was provided in [11]. Finally, the parameter $b$ is called as the batch-size parameter. It turns out that in the stochastic zeroth-order setting invariably we require $b>1$, which in turn leads to the (zeroth-order) oracle complexity being an order $b$ times that of iteration complexity. In Section 2 and 3.1, we use the above gradient estimation technique in the context of discretizing overdamped and underdamped Langevin diffusion and develop their oracle and iteration complexities. In order to establish the results, we will use the following Lemma due to [11] which provides an upper bound on the variance of $g_{\nu,b}$.

Lemma 1.1

[11]** Let $g_{\nu,b}$ be defined as in (6). Then under Assumption 1.3, and condition A1 of Assumption 1.1, we have,

[TABLE]

where $f_{\nu}(\theta)=\mathbf{E}_{u}[f(\theta+\nu u)]$.

One-point versus two-point evaluation: The gradient estimator in (6) is referred to as the two-point estimator in the literature. The reason is that, for a given random vector $\xi$, it is assumed that the stochastic function in (6) could be evaluated at two points,$F(\theta_{1},\xi)$ and $F(\theta_{2},\xi)$. Such an assumption is satisfied in several statistics, machine learning and simulation based optimization and sampling; see for example in [64, 53, 26, 2, 28, 36, 56]. Yet another estimator is the one-point estimator which assumes that for each $\xi$, we observe only one noisy function evaluation $F(\theta,\xi)$. Admittedly, the one-point setting is more challenging than the two-point setting. Specifically, in the one-point feedback setting, Lemma 1.1 no longer holds. From a theoretical point of view, the use of two-point evaluation based gradient estimator is primarily motivated by the sub-optimality (in terms of oracle complexity) of one-point feedback based stochastic zeroth-order optimization methods either in terms of the approximation accuracy or dimension dependency.

The use of one-point feedback for stochastic zeroth-order gradient estimation could be traced back to [57]. Motivated by this, there has been several works in the machine learning community focusing on leveraging it for zeroth-order convex optimization. Specifically, considering the class of convex functions (without any further smoothness assumptions) and adversarial noise (i.e., roughly speaking, with noise vectors not necessarily assumed to be independent and identically distributed (i.i.d.)), [14] proposed a polynomial-time algorithm and an oracle complexity of $\mathcal{O}(d^{21}/\epsilon^{2})$. This was improved to $\mathcal{O}(d^{5}/\epsilon^{2})$ recently in [47]. Further assuming Lipschitz smooth convex functions, [15] and [35], in the i.i.d noise case, obtained an oracle complexity of $\mathcal{O}(d^{7.5}/\epsilon^{2})$ and $\mathcal{O}(d/\epsilon^{3})$ respectively. The best known lower bound in this case is known to be $\mathcal{O}(d^{2}/\epsilon^{2})$, which was established by [62]. Further assuming $(\beta-1)$ differentiable derivatives, for $\beta>2$, [18] obtained as oracle complexity of $\mathcal{O}(d^{2}/\epsilon^{2\beta/(\beta-1)})$ and $\mathcal{O}(d^{2}/\epsilon^{(\beta+1)/(\beta-1)})$ respectively for the convex and strongly-convex setting, with i.i.d. noise case; see also [4]. In contrast to the above discussion, with two-point feedback it is possible to obtain much improved oracle complexities (i.e., linear in dimension and optimal in $\epsilon$) for stochastic zeroth-order optimization, as illustrated in [56, 36, 28, 2]. Given this subtle differences between the two-point and one-point evaluation settings for stochastic zeroth-order gradient estimation, in Section 4 we consider the effect of one-point gradient estimation technique for stochastic zeroth-order discretization of overdamped and underdamped Langevin diffusions.

1.3 Our Contributions

Under the availability of the stochastic zeroth-order oracles, we make the following contributions to the literature on sampling.

We first consider the case of strongly log-concave and smooth densities and analyze a stochastic zeroth-order version of Euler-discretization of overdamped and underdamped Langevin diffusions, under the two-point feedback setting in Section 2. For both cases, we characterize the oracle and iteration complexities to obtain $\epsilon$-approximate samples in term of $W_{2}$ metric. 2. 2.

We next consider in Section 2.3, a stochastic zeroth-order version of the recently proposed Randomized Midpoint Sampling method of the underdamped Langevin diffusion and characterize the oracle and iteration complexities to obtain $\epsilon$-approximate samples in term of $W_{2}$ metric. We show that for certain range of $\epsilon$, this method achieves improved oracle complexity compared to the above method. 3. 3.

While the above contributions are for strongly log-concave densities, in Section 3.1, we consider the more general class of densities satisfying log-Sobolev inequality and establish the oracle and iteration complexities of stochastic zeroth-order discretizations. 4. 4.

While all of the above contributions use the two-point stochastic zeroth-order feedback setting, in Section 4, we next consider the case of one-point feedback and characterize the corresponding oracle and iteration complexities for all the above discretizations. 5. 5.

Next, in Section 5, we consider variable selection for zeroth-order sampling. We specifically assume the unobserved function $f$ is sparse in the sense that it depends only on $s$ of the $d$ coordinates. We provide a variable selection method based on the estimated zeroth-order gradient, which in conjunction with the above discretizations reduces the oracle and iteration complexities to be only poly-logarithmically dependent on the dimensionality $d$ thereby enabling high-dimensional sampling.

Our contributions provide several theoretical insights on the performance of stochastic zeroth-order sampling algorithms, and widen the applicability of theoretically sound Bayesian inference to various practical situations where we do not know the analytical form of the potential function. All proofs are relegated to the appendix.

2 Oracle Complexity Results under Strong Log-concavity

We now leverage the stochastic zeroth-order gradient estimation methodology introduced in Section 1.2 for discretizing underdamped and overdamped Langevin diffusions. Throughout this section, we assume the target density is strongly log-concave and smooth (recall Assumption 1.1).

2.1 Zeroth-Order Langevin Monte Carlo

Replacing the true gradient in the first-order Langevin Monte Carlo algorithm in (3), with the zeroth-order gradient estimation in (6), we obtain the following Zeroth-Order LMC (ZO-LMC) algorithm:

[TABLE]

for $n=0,1,2,\cdots,N-1$. Apart from the choice of step-size $h$, the ZO-LMC also requires two additional tuning parameters, the smoothing parameter $\nu$ and the batch-size $b$ of the zeroth-order gradient estimator $b$, that need to be set. Although ZO-LMC could be interpreted as a form of LMC with inaccurate gradient as in [29], the corresponding theoretical result from [29] cannot be used directly for obtaining the oracle complexity of ZO-LMC, as the variance of the gradient in (6) is not bounded unless we make restrictive assumptions on the true gradient of $f$. We now state the main result of this section, which describes the oracle complexity of ZO-LMC.

Theorem 2.1

Let the potential function $f$ satisfy Assumption 1.1. Then, for the ZO-LMC algorithm in (9), under Assumption 1.3, by choosing

[TABLE]

we have $W_{2}(\varpi_{N},\pi)\leq\epsilon$ for $0<\epsilon\leq\min\left(d\sqrt{\frac{2}{M+m}},\sqrt{\frac{m(d+5)}{8M^{2}}}\right)$, after

[TABLE]

iterations. Hence, the total number of calls to the stochastic zeroth-order oracle is given by,

[TABLE]

Remark 2

Recall that for the exact gradient based LMC algorithm, to obtain $W_{2}(\varpi_{N},\pi)\leq\epsilon$, we require $N={\cal O}\left(d/\epsilon^{2}\cdot\log(d/\epsilon)\right)$ (see [29]) which matches (11). Thus, ZO-LMC matches the performance of LMC in terms of iteration complexity required to obtain $W_{2}(\varpi_{N},\pi)\leq\epsilon$. However, in each iteration of the LMC algorithm, we only require one gradient evaluation. Hence, the total number calls to the first-order oracle is also given by ${\cal O}\left(d/\epsilon^{2}\cdot\log(d/\epsilon)\right)$. For the ZO-LMC, in contrast we require $b=d$ calls to the stochastic zeroth-order oracle in each iteration. Hence, the oracle complexity is given by (12). By a straight forward modification of the proof of Theorem 2.1, for the ZO-LMC, if we restrict ourself to $b=1$, the iteration complexity increases to $N={\cal O}\left(d^{2}/\epsilon^{2}\cdot\log(d/\epsilon)\right)$, which will then also be the oracle complexity. Thus, the price we pay to match LMC in the absence of true gradient information is $O(d)$.

Remark 3

Recently [25] analyzed the standard Metropolis Random Walk algorithm (MRW), which is a zeroth-order algorithm, in the non-noise setting. Specifically, [25] showed that to achieve samples that are $\epsilon$-close to the target $\pi$ in total variation distance, MRW requires $\mathcal{O}(d^{2}\log(1/\epsilon)$ calls to the non-noisy zeroth-order oracles. Considering the non-noisy setting, the result appears to seemingly have an exponential improvement in terms of $\epsilon$. However, this result was obtained under the so-called warm start condition on the distribution of the initial vector $x_{o}$, which seems to be an opaque condition hiding the true complexity of the problem. For example, it is not clear how to pick such a warm start distribution for a given target $\pi$, in particular in the stochastic zeroth-order setting that we consider in this work. As a way to potentially avoid this opaque warm start condition, [25] suggests to set $x_{0}\sim N(x^{*},\mathbf{I}_{d})$, where $x^{*}$ is the unique minimizer of $f(x)$ and $\mathbf{I}_{d}$ is the $d\times d$ identity covariance matrix. For this choice of initial vector, to obtain a sample which is $\epsilon$-close to the target $\pi$ in total variation distance, [25] showed that MRW requires an oracle complexity of $\mathcal{O}(d^{2}\log(1/\epsilon)$. However, in the zeroth-order setting, the oracle complexity of finding an $\epsilon$-minimizer of a strongly-convex and smooth function $f(x)$, is well-studied problem in stochastic optimization – it is upper and lower bounded by $\mathcal{O}(d/\epsilon)$; see for example [28, 41, 56, 36]. This seems to negate the actual oracle complexity improvements shown in [25], as it really seems to require extremely careful initial distributions (i.e., knowledge of the exact minimzer), even in the non-noisy setting. Notwithstanding the fact that the results in [25] are for the non-noisy setting, they are essentially no better than the oracle complexity results established for ZO-LMC algorithm in Theorem 2.1, which also has the advantage that it does not require any opaque warm start conditions or special initial distributions.

2.2 Zeroth-Order Kinetic Langevin Monte Carlo

In the previous section, we consider the stochastic zeroth-order discretizations of the overdamped Langevin diffusions. It is known that in the first-order setting, discretizations of underdamped Langevin diffusion obtain improved oracle complexities [32, 21]. Under Langevin diffusion process (also called as kinetic Langevin diffusion process) is given by the following stochasic differential equation:

[TABLE]

where $\mathbf{I}_{d}$ is the $d\times d$ identity matrix. We refer the reader to [33, 21, 32] for more details about the above diffusion process and related theoretical results. Specifically, it was shown in [21, 32] that first-order discretizations of the kinetic diffusion process (referred to as KLMC in [21]) in (13) have better rates of convergence compared to similar first-order discretizations of the continous-process in (2). Specifically, recall that for the right choice of tuning parameters, LMC (i.e., first-order discretizations of (2)) requires that $N=\mathcal{O}(d/\epsilon^{2}\cdot\log(d/\epsilon))$ for $W_{2}(\varpi_{N},\pi)\leq\epsilon$. Whereas, it was shown in [21, 32] $N=\mathcal{O}(\sqrt{d}/\epsilon\cdot\log(d/\epsilon))$ suffices ([32] provides a much sharper result compared to [21]). We emphasize that the above result does not immediately imply that KLMC might be the algorithm to use always (in comparison to LMC); indeed when considering also the dependence of the bound on the strong-convexity and smoothness parameters (though the condition number of the sampling density defined as $M/m$), [32] precisely characterize when KLMC might be preferred over the vanilla LMC. The bottom line of their analysis is none of the method is uniformly better over the other method.

The Euler-discretization of the SDE in (13), which is a first-order sampling algorithm is given by the following iterations:

[TABLE]

where $(\tilde{\epsilon}_{n+1},\epsilon_{n+1})\in\mathbb{R}^{2d}$ is a a sequence of i.i.d standard Normal vectors, independent of $(\tilde{x}_{0},x_{0})$ and $\psi_{0}(t)=e^{-\gamma t}$ and $\psi_{n+1}=\int_{0}^{T}\psi_{n}(s)ds$. We refer to this algorithm as KLMC following the terminology of [32]. Based on this, we now consider the ZO-KLMC updates as:

[TABLE]

where $g_{\nu,b}$ is the zeroth-order gradient estimator as in (6). In comparison to the ZO-LMC algorithm, the ZO-KLMC algorithm has an additional tuning parameter $\gamma$ that needs to be set. For the ZO-KLMC algorithm, we have the following complexity result.

Theorem 2.2

Let the potential function $f$ satisfy Assumption 1.1. If the initial point $(\tilde{x}_{0},x_{0})$ is chosen such that $\tilde{x}_{0}\sim N(0,\bm{I}_{d})$, then, ensuring $\gamma\geq\sqrt{m+M}$, for the ZO-KLMC, under Assumption 1.3, by choosing,

[TABLE]

we have $W_{2}(\varpi_{N},\pi)\leq\epsilon$ for $0<\epsilon\leq\frac{12M\gamma^{2}\sqrt{d}}{m^{2}}$ after

[TABLE]

iterations. Here $\tilde{\mathcal{O}}$ hides poly-logarithmic factors in $1/\epsilon$. Hence, the total number of calls to the stochastic zeroth-order oracle is given by

[TABLE]

Remark 4

We note that compared to ZO-LMC, while ZO-KLMC obtains improved iteration complexity, the iteration complexity still remains the same. The improvement in the iteration complexity is indeed a consequence of a similar improvement in the first-order setting as demonstrated in [21, 32].

2.3 Zeroth Order Randomized Midpoint Method

Given that the ZO-KLMC offers no improvement over ZO-LMC in terms of oracle complexity despite its improved iteration complexity, it is worth examining if there are other discretizations that obtain improvements in oracle complexities. Towards that, in this section we analyze the zeroth-order version of the Randomized Mid-Point discretization of the underdamped Langevin diffusion, proposed in [63]. In the first-order setting, [22] recently showed that the Randomized Mid-Point discretization of underdamped Langevin diffusion achieves the information theoretic lower bounds for sampling. See also [40] for additional probabilistic results.

The crux of the randomized midpoint method is based on first representing the kinetic Langevin Monte Carlo in (13) in its integral format, and estimating the integrals based on a randomization technique. We also mention in passing that the randomized midpoint idea shares some similarities to symplectic integration methods [65] from the sampling literature and extragradient method [44] from the optimization literature, with the main difference being the randomized choice of step-size which leads to improved oracle complexities. We now provide the algorithm in the zeroth-order setting and the corresponding theoretical result. Let $\epsilon^{(i)}_{n}\in\mathbb{R}^{d}$, $i=1,2,3$, be a sequence of Gaussian random vectors generated according to the procedure described in Appendix A of [63]. Let $\alpha_{n}$ be a sequence of uniform random variables supported on the interval $[0,1]$. Then the zeroth-order Randomized Mid-Point Method (ZO-RMP) is given by the following updates:

[TABLE]

We remark that we use the same choice of batch-size, $b$, in (19), (20) and (21), as using different batch sizes has no effect on the oracle complexity.

Theorem 2.3

Define $\kappa=M/m$ to be the condition number of the potential $f$ which satisfies Assumption 1.1. Furthermore, let the stochastic zeroth-order oracle satisfy Assumption 1.3. Let $x^{*}$ be the minimizer of $f$, and $x_{0}$ be such that $\mathbf{E}\left[f(x_{0})-f(x^{*})\right]=O(d)$, and $v_{0}=0$. Then, for $0\leq\epsilon\leq 1$, by choosing,

[TABLE]

for the ZO-RMP method described in (19)-(21), with $u=1/M$, we have $W_{2}(\varpi_{N},\pi)\leq\epsilon$ after

[TABLE]

iterations. Hence, the total-number of zeroth-order oracle calls are given by

[TABLE]

Remark 5

The analysis of the randomized midpoint algorithm in [63], for the first-order setting, requires access to exact minimizer $x^{*}$ as the initializer. We relax this requirement to the having a point $x_{0}$ satisfying $\mathbf{E}\left[f(x_{0})-f(x^{*})\right]=O(d)$, which is a milder requirement. It is well-known from the stochastic optimization literature, that under Assumption 1.3, and 1.1, in the zeroth-order setting, using the zeroth-order version of stochastic gradient algorith, the oracle complexity of finding a point $x_{0}$ such that $\mathbf{E}\left[f(x_{0})-f(\bar{x})\right]=O(d)$ where $\bar{x}$ is the minimizer of $f$, is $O(\kappa\log d)$ [28, 56].

Remark 6

Note that even though the iteration complexity of ZO-RMP still matches with RMP (except for the dimension dependence which is unavoidable in the zeroth-order setting), and is better than KLMC for all values of $\epsilon$, the oracle complexity for ZO-RMP is not uniformly better than ZO-KLMC for all $\epsilon$. However, observe that when $h=C\frac{(\epsilon\sqrt{m})^{\frac{1}{3}}}{(d\kappa)^{\frac{1}{6}}\log\left(\frac{1}{\epsilon}\right)^{\frac{1}{6}}}$, i.e., when $\epsilon\geq\max\left(\sqrt{\frac{d}{M}},\frac{16\sigma^{2}}{M^{\frac{3}{2}}\sqrt{d}},\frac{1}{\sqrt{dmM}}\right)$ the oracle complexity of ZO-RMP is $\tilde{O}\left(\frac{d^{\frac{5}{3}}\kappa^{\frac{8}{3}}}{\epsilon^{\frac{4}{3}}}\right)$ which is indeed better compared to $\tilde{O}\left(\frac{d^{2}}{\epsilon^{2}}\right)$ for ZO-KLMC.

We end this section by mentioning that developing lower bounds on the oracle complexity of sampling from strongly log-concave densities, in the stochastic zeroth-order setting that we consider is an interesting open problem.

3 Oracle Complexity Results under Log-Sobolev Inequality

The algorithms and oracle complexity results in the previous sections were stated for smooth and strongly log-concave densities (i.e., under Assumption 1.1), which covers important classes of problems in sampling and Bayesian inference. However, the fundamental idea behind the non-asymptotic convergence results of the discretization based sampling algorithm are essentially based on the following facts: (i) the underlying continuous (underdamped or overdamped) Langevin diffusion converges to its equilibrium state (i.e., to the target distribution $\pi$ in this case) exponentially fast in various metrics, and (ii) consequently, the potential function is smooth enough that the error due to discretization is not extremely large. Roughly speaking, condition A1 and A2 in Assumption 1.1 corresponds respectively to the above facts, respectively. However, it is well-know that the overdamped Langevin diffusion converges to its equilibrium under much weaker conditions that strong log-concavity; indeed as long as the target density satisfies functional inequalities like Poincare or Log-Sobolev inequalities, the overdamped Langevin diffusion converges to its equilibrium exponentially faster in various metrics; see, for example [12]. Motivated by the above fact, recently [72] demonstrated that the LMC algorithm also exhibits rapid convergence to the target density if it has access to the exact gradient evaluations of the potential function $f$. As a consequence, one could sample from densities that are not essentially strongly log-concave, there by extending the applicability of LMC algorithms for a wider class of Bayesian inference problems. In this section, we analyze the performance of stochastic zeroth-order discretization of overdamped Langevin diffusions when the target density satisfies log-Sobolev inequality.

Assumption 3.1

A density $\pi$ is said to satisfy Log-Sobolev Inequality (LSI) with a constant $\lambda>0$ if for all smooth function $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ with finite variance,

[TABLE]

In Section LABEL:sec:lsinumerics, we show that mixture of Gaussian densities with unequal covariance satisfies the above assumption, while it does not satisfy condition A1 of Assumption 1.1, and discuss applications to Bayesian variable selection. The above assumption also leads to the following equivalent formulation. Let $H_{\pi}(\varpi)$, and $J_{\pi}(\varpi)$ be the Kullback-Leibler (KL) divergence of $\varpi$ with respect to $\pi$, and the relative Fisher Information respectively which are defined as follows:

[TABLE]

One can verify that LSI is equivalent to the following condition by plugging $g^{2}=\varpi/\pi$ in (25):

[TABLE]

We also know that when $\pi$ satisfies LSI, Talagrand inequality holds [12], i.e., for all $\varpi$,

[TABLE]

With this background, we provide our oracle complexity result of ZO-LMC algorithm when the density satisfies LSI and is smooth.

Theorem 3.2

Let the target density $\pi$ satisfy Assumption 3.1 and let the potential function $f$ be satisfy condition A2 in Assumption 1.1. Let $x_{0}\sim\varpi_{0}$ which satisfies $H_{\pi}(\varpi_{0})\leq\infty$. Then for the ZO-LMC update as in (9), under Assumption 1.3, by choosing,

[TABLE]

we have $W_{2}(\varpi_{N},\pi)\leq\epsilon$, for all $0\leq\epsilon\leq\frac{\alpha}{4L^{2}}$, after $N$ iterations where

[TABLE]

Hence, the total number of calls to the stochastic zeroth-order oracle is given by

[TABLE]

Remark 7

Note that in comparison to condition A1 of Assumption 1.1, the assumptions required for the above theorem are weaker. Specifically, in place of condition A1 in Assumption 1.1, we have Assumption 3.1. Condition A2 is regarding the smoothness is required to handle error that arises due to discretization of continuous time dynamics. For this wider class of densities, the price to pay is that the dependency on both the dimension $d$ and $\epsilon$ increases in comparison to Theorem 2.1.

Remark 8

Given that ZO-LMC exhibits convergence (albeit with a slightly weaker $\epsilon$ and $d$ dependency, it is natural to ask if ZO-KLMC also exhibits similar convergence. However, even in the first-order setting this question is open. Indeed, kinetic Langevin diffusions are a class of degenerate diffusions which require a different class of function inequalities (called as hypocoercivity [71]) for them to converge to their equilibrium. It is an open question to show that the discretize sampling algorithm (KLMC or appropriate modifications) also convergence under hypocoercivity and appropriate smoothness assumptions on the potential function $f$, either given access to exact first-order oracles or stochastic zeroth-order oracles. We leave this question as future work.

4 One-Point Setting: Independent noise per function evaluation

As discussed in Section 1.2, there are subtle differences between the availability of one and two-point evaluation based stochastic zeroth-order gradients. In this section, we examine this difference in more detail. Recall that while defining the zeroth-order gradient estimator in (6), we assumed that the function can be evaluated at two points, namely, $\theta+\nu u_{i}$, and $\theta$, with the same noise $\xi_{i}$. This implies, when the noise is additive, i.e., $F(\theta,\xi)=f(\theta)+\xi$, the gradient estimator is not affected by the noise. Because in that case we have, $F(\theta+\nu u_{i},\xi_{i})-F(\theta,\xi_{i})=f(\theta+\nu u_{i})-f(\theta)$. We emphasize that this is our main reason for consider general non-additive noise in the previous sections. For example, under multiplicative noise, consider the case where $F(\theta,\xi)=\xi f(\theta)$, $\mathbf{E}\left[\xi\right]=1$, and $f(\theta)$ is $L$-Lipschitz continuous; then Assumption 1.3 holds.

Now we will examine the one-point setting in that the noise in the two function evaluations of the gradient estimator is not the same. Specifically, first we show that allowing the noise $\xi_{i}$, and $\xi_{i}^{\prime}$ in $F(\theta+\nu u_{i},\xi_{i})$, and $F(\theta,\xi_{i}^{\prime})$ to be independent additive noise, deteriorates the iteration and/or oracle complexities of zeroth-order discretizations considered in the previous settings. Formally, we work under the following assumption in the one-point stochastic zeroth-order setting.

Assumption 4.1

The stochastic zeroth-order oracle is such that for each point $x$, the observed function evaluation $F(\theta,\xi)$ is given by $F(\theta,\xi)=f(\theta)+\xi$ where $\mathbf{E}\left[\xi\right]=0$, and $\mathbf{E}\left[\xi^{2}\right]=\sigma^{2}$.

Under Assumption 4.1, the upper bound on the variance of the gradient estimator as stated in Lemma 1.1 no longer holds. Instead, we have the following result.

Lemma 4.1

Let $g_{\nu,b}(\theta)$ in (6), be defined under the one-point setting. Then under Assumption 4.1 and condition A1 of Assumption 1.1, we have

[TABLE]

The main difference in the one-point setting, in terms of the variance of the gradient estimator is the presence of the third term, which is of the order of $1/b\nu^{2}$. This causes the additional difficulties in terms of setting the parameters $b$ and $\nu$ in the zeroth-order gradient estimator, which in turn causes the oracle complexities to deteriorate. Based on the above result on the variance, we provide the oracle complexity results for ZO-LMC, ZO-KLMC and ZO-RMP under Assumption 4.1 on the stochastic zeroth-order oracle, in Theorem 4.2, 4.3 and 4.4 respectively.

Theorem 4.2** (ZO-LMC under Strong Log-concavity)**

Let the potential function $f$ satisfy Assumption 1.1. Then, for ZO-LMC algorithms under Assumption 4.1, by choosing

[TABLE]

we have $W_{2}(\varpi_{N},\pi)\leq\epsilon$ for $0<\epsilon\leq\min\left(d\sqrt{\frac{2}{M+m}},\sqrt{\frac{m(d+5)}{8M^{2}}}\right)$, after $N$ iterations, where

[TABLE]

Hence, the total number of calls to the stochastic zeroth-order oracle is given by,

[TABLE]

Theorem 4.3** (ZO-KLMC under Strong Log-concavity)**

Let the function $f$ satisfy Assumption 1.1. If the initial point $(\tilde{x}_{0},x_{0})$ is chosen such that $\tilde{x}_{0}\sim N(0,\bm{I}_{d})$, then, ensuring $\gamma\geq\sqrt{m+M}$, under Assumption 4.1 for the ZO-KLMC, by choosing,

[TABLE]

we have $W_{2}(\varpi_{N},\pi)\leq\epsilon$ for $0<\epsilon\leq\frac{12M\gamma^{2}\sqrt{d}}{m^{2}}$, after

[TABLE]

iterations. Hence, the total number of oracle calls to the stochastic zeroth-order oracle is given by

[TABLE]

Theorem 4.4** (ZO-RMP under Strong Log-concavity)**

Let the potential function satisfy Assumption 1.1 and let $x^{*}$ be the minimizer of $f$, $x_{0}$ be such that $\mathbf{E}\left[f(x_{0})-f(\bar{x})\right]=O(d)$, and $v_{0}=0$. Let the stochastic zeroth-order oracle satisfy Assumption 4.1. Then, for $0\leq\epsilon\leq 1$, by choosing,

[TABLE]

for the ZO-RMP described in (19)-(21), we have $W_{2}(\varpi_{N},\pi)\leq\epsilon$ after

[TABLE]

iterations. Hence, the total-number of zeroth-order oracle calls are given by

[TABLE]

Remark 9

As before, the oracle complexity of ZO-RMP in this setting is not uniformly better than that of ZO-KLMC. We do observe that when $h=C\frac{(\epsilon\sqrt{m})^{\frac{1}{3}}}{(d\kappa)^{\frac{1}{6}}\log\left(\frac{1}{\epsilon}\right)^{\frac{1}{6}}}$, i.e., when $\epsilon\geq\max\left(\sqrt{\frac{d}{M}},\frac{16\sigma^{2}}{M^{\frac{3}{2}}\sqrt{d}},\frac{1}{\sqrt{dmM}}\right)$ the oracle complexity of ZO-RMP is $\tilde{O}\left(\frac{d^{\frac{16}{3}}\kappa^{\frac{10}{3}}}{\epsilon^{\frac{8}{3}}}\right)$ which is worse compared to $\tilde{O}\left(\frac{d^{2}}{\epsilon^{4}}\right)$ for ZO-KLMC. However, it is better than that of ZO-KLMC in the opposite regime.

We now present the corresponding result when the target density is not strongly log-concave but satisfies LSI.

Theorem 4.5** (ZO-LMC under Log-Sobolev Inequality)**

Let the target density $\pi$ satisfy Assumption 3.1 and let the potential function $f$ satisfy condition A2 of Assumption 1.1. Let $x_{0}\sim\varpi_{0}(x)$ which satisfies $H_{\pi}(\varpi_{0})\leq\infty$. Then for the ZO-LMC update as in (9), under Assumption 4.1, by choosing,

[TABLE]

we have, $W_{2}(\varpi_{N},\pi)\leq\epsilon$, for all $0\leq\epsilon\leq\frac{\lambda}{4L^{2}}$ after

[TABLE]

iterations. Hence, the total number of calls to the zeroth-order oracle is given by

[TABLE]

Remark 10** (One-point setting with non-additive noise )**

Given the above result, it is natural to examine the effect of non-additive noise on the oracle complexities. For this case, we have the following result under an additional smoothness assumption on the stochastic function evaluations $F(x,\xi)$.

Lemma 4.2

Let the function $F(\theta,\xi)$ be Lipschitz continuous in its second argument, i.e., $|F(\theta,\xi)-F(\theta,\xi^{\prime})|\leq L|\xi-\xi^{\prime}|$. Under the above condition, Lemma 4.1 holds. Consequently, all the above complexity results in this section holds.

Remark 11** (Effect of Higher-order smoothness)**

While the oracle complexities under the one-point evaluation setting are worse than that of the two-point setting, they could be made to approach that of the two-point setting when we make the stronger assumption that the potential function is assumed to the $\beta$-times differentiable and the $(\beta-1)$-derivatives are Lipschitz continuous. Similar phenomenon has been observed in the case of highly-smooth convex stochastic zeroth-order optimization; see, for example [18, 4]. As the precise statements and the proofs are similar to that of the above theorems, we omit the details.

5 Variable Selection for High-dimensional Black-box Sampling

In practical black-box settings, due to the non-availability of the analytical form of $f(\theta)$, one might potentially over-parametrize $f(\theta)$, in terms of number of covariates selected for modeling. Hence, the problem of variable selection, in a zeroth-order setting becomes crucial. To address this issue, in this section, we study variable selection under certain sparsity assumptions on the objective function $f$, to facilitate sampling in high-dimensions. Throughout this section, we assume one could observe exact function evaluations, without any noise. We emphasize that we make this assumption purely for technical convenience and to convey the theoretical results insightfully; all results presented in this section extends to the noisy setting in a straightforward manner. Specifically, we make the following assumption on the structure of $f$.

Assumption 5.1

We assume that $f(\theta):\mathbb{R}^{d}\to\mathbb{R}$ is $s$ sparse, i.e., the function $f$ depends only on (the same) $s$ of the $d$ coordinates, for all $\theta$, where $s\ll d$. We denote the true support set as $S^{*}$. This implies that for any $\theta\in\mathbb{R}^{d}$, we have $\|\nabla f(\theta)\|_{0}\leq s$, i.e., the gradient is $s$-sparse. Furthermore, define $\nabla f_{\nu}(\theta)={\bf E}_{u}\left[\nabla f(\theta+\nu u)\right]$ for a standard gaussian random vector $u$. Then the gradient sparsity assumption also implies that $\|\nabla f_{\nu}(\theta)\|_{0}\leq s$ for all $\theta\in\mathbb{R}^{d}$. Furthermore, we assume that the gradient lies in the following set that characterizes the minimal signal strength in the relevant coordinates of the gradient vector:

[TABLE]

As a consequence, we also have that $\nabla f_{\nu}(\theta)\in\mathcal{G}_{a,s}$. The above assumption makes a homogenous sparsity assumption on the sparsity and the minimum signal strength of the gradient. Roughly speaking, $a$ represents the minimum signal strength in the gradient so that efficient estimation of the support $S^{*}$ is possible in the sample setting. The above sparsity model on the function $f$, converts the problem to variable selection in a non-Gaussian sequence model setting:

[TABLE]

Hence, $\zeta_{j}$ are zero-mean random variables as $[g_{\nu,n}]_{j}$ is an unbiased estimator of $[\nabla f_{\nu}(\theta)]_{j}$. We refer the reader to [17] for recent results on variable selection consistency in Gaussian sequence model setting. We also make the following assumption on the query point selected to estimate the gradient.

Assumption 5.2

The query point $\theta\in\mathbb{R}^{d}$ selected is such that $\|\nabla f(\theta)\|_{2}\leq R$.

Our algorithm for high-dimensional black-box sampling with variable selection is as follows:

• •

Pick a point $\theta$ (which is assumed to satisfy Assumption 5.2) and estimate the gradient $g_{\nu,n}$ at that point and compute the estimator $\hat{S}$ of $S^{*}$ as $\hat{S}=\{j:|[g_{\nu,n}]_{j}|\geq\tau\}$.

• •

Run any of the zeroth-order sampling algorithm on the selected set of coordinates $\hat{S}$ of $f(\theta)$.

Here, for the first step, we need to select $n,\tau$ and $\nu$. We separate the set of relevant variables by thresholding $|[g_{\nu,n}]_{j}|$ at $\tau$. We now provide our result on the probability of erroneous selection.

Theorem 5.3

Let $f$ satisfy Assumption 1.1 and the query point selected satisfy Assumption 5.2. Set $\tau=(a-M\nu\sqrt{s})/2$ and assume that $\nu\leq\min\left(\frac{a}{2M\sqrt{s}},\frac{R}{MC_{2}\sqrt{s}}\right)$ and

[TABLE]

where $C,C_{2}$ are constants. Then we have $\Pr\{\hat{S}\neq S^{*}\}\leq\epsilon$.

Remark 12

The number of queries $n$ to the function $f$ depends only logarithmically on the dimension $d$ and is a (low-degree) polynomial in the sparsity level $s$. Combining this fact with the result in Theorem 2.1 we see that the total number of queries to the function $f$ (for the sampling error measured in 2-Wasserstein distance) is only poly-logarithmic in the true dimension $d$ and is a low-degree polynomial in the sparsity level $s$. Thus when $s\ll d$, we see the advantage of variable selection in black-box sampling using the two-step approach. The above results assumes that the sparsity level $s$ and signal strength is known. It would be interesting to construct adaptive estimators similar to those for Gaussian sequence model in [17]. Furthermore, exploring appropriately defined notions of non-homogenous sparsity assumptions is also challenging.

6 Discussion

In this work, we proposed and analyzed zeroth-order discretizations of overdamped or underdamped Langevin diffusions. We provide a through analysis of the oracle complexity of such sampling algorithms under various noise models on the zeroth-order oracle and provide simulation results corroborating the theory. Recall that our zeroth-order gradient estimators used in this work were based on Gaussian Stein’s identity and could be used for the case when $f$ is defined on the entire Euclidean space $\mathbb{R}^{d}$.In several situation, for example, in sampling from densities with compact support [7, 9] and in computing volume of convex body [13], one needs to compute the gradient of the function (and density) supported on $\mathcal{M}\subset\mathbb{R}^{d}$. For these situations, one can use a version of Stein’s identity based on score functions to compute the gradient and Hessian. To explain more, we first recall some definitions. The score function $S_{p}\colon\mathcal{M}\rightarrow\mathbb{R}^{d}$ associated to density $p(u)$ defined over $\mathcal{M}$ is defined as

[TABLE]

In the above definition, the derivative is taken with respect to the argument $u$ and not the parameters of the density $p(u)$. Based on the above definition, we have the following versions of Stein’s identity; see, for example, [37].

Proposition 6.1

Let $U$ be a $\mathcal{M}$-valued random vector with density $p(u)$. Assume that $p\colon\mathcal{M}\rightarrow\mathbb{R}$ is differentiable. In addition, let $g:\mathcal{M}\to\mathbb{R}$ be a continuous function such that ${\bf E}_{U}[\nabla g(U)]$ exists and the following is true: $\int_{u\in\mathcal{M}}\nabla_{u}\left(g(u)p(u)\right)du=0$. Then it holds that

[TABLE]

where $S(u)=-\nabla p(u)/p(u)$ is the score function of $p(u)$.

In order to leverage the above identities to estimate the gradient of a given function $f(\theta):\mathcal{M}\to\mathbb{R}$, consider $g(U)=f(\theta+U)$ where $U\sim p(u)$ is a $\mathcal{M}$-valued random variable and appeal to the above Stein’s identity above, as done in Section 2 for with Gaussian random variables. A special case of the above idea, when the space $\mathcal{M}$ is a Riemannian sub-manifold embedded in an Euclidean space was considered in [48] in context of stochastic zeroth-order Riemannian optimization. We postpone a rigorous analysis of the estimation and approximation rates in the general setting, and their applications to black-box sampling on non-Euclidean spaces for future work.

7 Notations

We use $a\wedge b$ and $a\vee b$ to denote the minimum and maximum of $a$ and $b$ respectively. The $L_{2}$ norm of a random vector $X:\Omega\to\mathbb{R}^{d}$ is defined to be $\|X\|_{L_{2}}={\bf E}[\|X\|_{2}^{2}]^{1/2}$. The $L_{p}$ norms of a random matrix $\bm{M}:\Omega\to\mathbb{R}^{d\times d}$ are defined as follows.

[TABLE]

where $\|\cdot\|_{2}$ is the spectral norm, and $\|\cdot\|_{F}$ is the Frobenius norm. For simplicity, we write $\|\cdot\|=\|\cdot\|_{2}$ and $\|\cdot\|_{L_{p}}=\|\cdot\|_{L_{p},\bm{\bullet}}$ when there is no ambiguity. Furthermore, we omit the subscript $h$ in $x_{t,h}$ in places where is no confusion for simplicity.

8 Proofs for Section 2

8.1 Proofs for Oracle Complexity of ZO-LMC

Proof. [of Theorem 2.1] The proof follows by first calculating the bias and variance of the gradient estimator in our zeroth-order setting, where the error term $\zeta_{n}=g_{\nu,b}(x_{n})-\nabla f(x_{n})$. First, by Stein’s identity, ${\bf E}[g_{\nu,1}(x,u)]={\bf E}[\nabla f(x+\nu u)]=\nabla f_{\nu}(x)$, where we denote $f_{\nu}(x)={\bf E}[f(x+\nu u)]$. Under Assumption 1.1 on smoothness of $f$, in the case where $b=1$, we have the following calculation for the bias.

[TABLE]

Next, for $b\geq 1$ in general, $g_{\nu,b}(x)=\frac{1}{b}\sum_{k=1}^{b}g_{\nu,1}(x,u_{k})$, the bias and variance could be calculated as follows. Specifically, for the bias, we have

[TABLE]

From Lemma 2.1 of [11], we have,

[TABLE]

Next, we follow a similar framework to the proof of Theorem 4 in [29], but with modifications to adapt to the variance that is not uniformly bounded. Recall that $\Delta_{n}=L_{0}-x_{n},\;\Delta_{t+1}=L_{h}-x_{t+1}$, where $L_{n}=L_{0}-\int_{0}^{T}\nabla f(L_{s})ds+\sqrt{2}W_{n}$ follows the Langevin diffusion with stationary distribution $\pi$. Moreover, $\|\Delta_{n}-hU\|=\|\Delta_{n}-h[\nabla f(x_{n}+\Delta_{n})-\nabla f(x_{n})]\|\leq(1-mh)\|\Delta_{n}\|,\;\|V\|=\|\int_{0}^{h}[\nabla f(L_{s})-\nabla f(L_{0})]ds\|\leq 1.65M(h^{3}d)^{1/2}$. Thus,

[TABLE]

Here we use the fact that $\sqrt{a^{2}+b+c}\leq\sqrt{a^{2}+b}+\frac{c}{2a}$, ${\bf E}[\|\nabla f(L)\|^{2}]\leq Md$, and that we choose $h$, and $b$ such that $h/\left(b(1-mh)\right)\leq m/(8M^{2}(d+5))$. By Lemma 9 in [29], the above inequality leads to

[TABLE]

Therefore, using the fact $W_{2}(\varpi_{n+1},\pi)\leq\|\Delta_{n+1}\|_{L_{2}}$, and $W_{2}(\varpi_{0},\pi)=\|\Delta_{0}\|_{L_{2}}$, we obtain the bound in Wasserstein distance.

[TABLE]

Choosing $h$, $b$, $\nu$, and $N$ as in (10), and (11) we have $W_{2}(\varpi_{N},\pi)\leq\epsilon$.

8.2 Proofs for Oracle Complexity of ZO-KMLC

Proof. [of Theorem 2.2] Let $(V_{n,n},L_{n,t}),\;t\in[0,h]$ be a stationary kinetic Langevin process for each $n\in\mathbb{N}$, i.e.,

[TABLE]

starting from $V_{0,0}\sim N(0,\bm{I}_{d}),\;L_{0,0}\sim\pi$, and satisfying $V_{n,h}=V_{n+1,0},L_{n,h}=L_{n+1,h}$. Define $(\tilde{V}_{n,t},\tilde{L}_{n,t})$ by the following discretized version of kinetic Langevin diffusion,

[TABLE]

or equivalently,

[TABLE]

Define a different kinetic Langevin process $(\hat{V}_{n,t},\hat{L}_{n,t})$ with initial condition $\hat{V}_{n,0}=\tilde{V}_{n,0},\hat{L}_{n,0}=\tilde{L}_{n,0}$, i.e.,

[TABLE]

Assume that $(\tilde{V}_{0,0},\tilde{L}_{0,0})$ is chosen such that $\tilde{V}_{0,0}=V_{0,0}$ and $W_{2}(\varpi_{0},\pi)=\|\tilde{L}_{0,0}-L_{0,0}\|_{L_{2}}$. By definition of Wasserstein distance, we have $W_{2}(\varpi_{n},\pi)\leq\|\tilde{L}_{n,0}-L_{n,0}\|_{L_{2}}$.

Now we denote $e_{n}=\left\|\bm{P}^{-1}\begin{bmatrix}\tilde{V}_{n,0}-V_{n,0}\\ \tilde{L}_{n,0}-L_{n,0}\end{bmatrix}\right\|_{L_{2}}$, where $\bm{P}^{-1}=\begin{bmatrix}\bm{I}_{d}&\gamma\bm{I}_{d}\\ -\bm{I}_{d}&\bm{0}\end{bmatrix},\;\bm{P}=\gamma^{-1}\begin{bmatrix}0&-\gamma\bm{I}_{d}\\ \bm{I}_{d}&\bm{I}_{d}\end{bmatrix}$ corresponds to the contraction to the kinetic Langevin process. See [32]. Note that $\|\tilde{L}_{n,0}-L_{n,0}\|_{L_{2}}\leq\sqrt{2}\gamma^{-1}e_{n}$ and $\|\tilde{V}_{n,0}-V_{n,0}\|_{L_{2}}\leq e_{n}$. Observe that,

[TABLE]

Similarly,

[TABLE]

Combining (46), and (47), we have

[TABLE]

Now we will upper bound the above three terms. Observe that,

[TABLE]

[TABLE]

So, combining (49), and (50), and using the fact $\|\hat{V}_{n,u}\|_{L_{2}}\leq\|V_{n,u}\|_{L_{2}}+\|\hat{V}_{n,u}-V_{n,u}\|_{L_{2}}\leq\sqrt{d}+e_{n}$, and choosing $h\leq\sqrt{2}/(10\gamma)$, we obtain

[TABLE]

Using (44) we have

[TABLE]

and

[TABLE]

Note that,

[TABLE]

Using Lemma 1.1, we also have,

[TABLE]

Using the fact that $\|\hat{L}_{n,0}-L_{n,0}\|_{L_{2}}^{2}\leq 2\gamma^{-2}e_{n}^{2}$, and $\|\nabla f(L_{n,0})\|_{L_{2}}^{2}\leq Md$, we hence obtain

[TABLE]

where $A_{4}=\frac{2(d+5)(2Md+\sigma^{2})}{b}+\frac{\nu^{2}M^{2}(d+3)^{3}}{2b}$. Similarly, we have

[TABLE]

So, using (54), and (55), we have

[TABLE]

Now using (53), (56), and using the facts that, $\mathbf{E}\left[A_{2}|\hat{L}_{n,0}\right]=0$, $\mathbf{E}\left[B_{2}|\hat{L}_{n,0}\right]=0$, $\hat{V}_{n,h}-V_{n,h}$ is independent of $A_{2},B_{2}$ given $\hat{L}_{n,0}$, and $\hat{L}_{n,h}-L_{n,h}$ is independent of $A_{2},B_{2}$ given $\hat{L}_{n,0}$, we get

[TABLE]

The second inequality follows as $h\leq\sqrt{2}/(10\gamma)$, the third inequality follows if we choose $h\leq\min(\gamma/m,\sqrt{2}/(10\gamma))$, and the last inequality follows if we choose $b\geq\frac{512M^{2}(d+5)}{3m^{2}}$. Combining, (48), (51), (52), and (57), we get

[TABLE]

Using Lemma 9 of [29], and choosing $h\leq\min(\gamma/m,m/(12\gamma M))$ we have $mh/(4\gamma)-3Mh^{2}/2\geq mh/(8\gamma)$, and thus

[TABLE]

Then we obtain

[TABLE]

Now, choosing $h$, $\nu$, $b$, and $N$ as in (16), we get (17), and (18).

8.3 Proofs for Oracle Complexity of ZO-RMP

Before proceeding, we also recall that $\left(x_{n}^{*}(t),v_{n}^{*}(t)\right)$ when $t\in[0,h]$ is the true solution to the underdamped Langevin diffusion with the initial point $(x_{n},v_{n})$ coupled with $x_{n+\frac{1}{2}}$ through a shared Brownian motion defined as follows:

[TABLE]

We also recall some preliminary results from [63].

Lemma 8.1** (Lemma 6[63])**

Let $\{x(t)\}_{t\in[0,h]}$, and $\{v(t)\}_{t\in[0,h]}$ be the true solution to the underdamped Langevin diffusion (58), and (59) on ${t\in[0,h]}$. Then for $h\leq 1/20$, and $u=1/M$, we have

[TABLE]

Lemma 8.2

Let $\alpha_{n}$ be sampled uniformly randomly from $[0,1]$ at iteration $n$. Let $x_{n+\frac{1}{2}}$ be the intermediate value at step $n$. Let $\{x_{n}^{*}(t)\}_{t\in[0,h]}$ be the true solution to (58), and (59) with the initial point $x_{n}^{*}(0)=x_{n}$ coupled to $x_{n+\frac{1}{2}}$ through a shared Brownian motion. Then, under Assumption 1.3 and Assumption 1.1, for $h\leq 1/20$, we have

[TABLE]

Proof. [of Lemma 8.2] For notational simplicity, we drop the subscript $n$ from $\alpha_{n}$ below. First, note that we have

[TABLE]

The first and the sixth inequality follows from the first condition of Assumption 1.1, the second inequality follows from (19), and (58), the third inequality follows from Cauchy-Schwarz inequality, the fourth inequality follows from choosing $u=1/M$ and the fact that $1-e^{-2(\alpha h-s)}\leq 2h$, the fifth inequality follows from Young’s inequality, the seventh inequality follows from Lemma 1.1 and Lemma 8.1. Choosing $b$, and $\nu$ as in (22), we have,

[TABLE]

Lemma 8.3

Let $g_{\nu,b}(x_{n})$ be defined as in (6). Then under the conditions of Lemma 8.2, we have

[TABLE]

Proof. [of Lemma 8.3] Using Lemma 1.1 and Young’s inequality, we have

[TABLE]

Furthermore, using Lemma 8.2, and (61), and the fact that $h$ is small, we get

[TABLE]

Hence, we have

[TABLE]

Lemma 8.4

Let $\mathbf{E}_{\alpha}$ denote the expectation with respect to $\alpha$ at each iteration $n$. Let $\mathbf{E}\left[\cdot\right]$ be the expectation with respect to other randomness present in iteration $n$. Let $\{x_{n}^{*}(t)\}_{t\in[0,h]}$ be the true solution to (58), and (59) with the initial point $x_{n}^{*}(0)=x_{n}$ coupled to $x_{n+\frac{1}{2}}$, $v_{n}$, and $x_{n+1}$ through a shared Brownian motion. Then, under Assumption 1.3–1.1, for $h\leq 1/20$, and $u=1/M$, we have

[TABLE]

Proof. [of Lemma 8.4]

1. a)

Using Lemma 8.2, and 8.3, we have

[TABLE]

The second inequality follows from Young’s inequality, the third and fifth inequality uses the fact $1-e^{-2(\alpha-\alpha h)}\leq 2h$, the fifth inequality follows from the fact $\frac{uh}{2}\mathbf{E}_{\alpha}(1-e^{-2(h-\alpha h)})\nabla f(x_{n}^{*}(\alpha h))-\frac{u}{2}\int_{0}^{h}(1-e^{-2(h-s)})\nabla f(x_{n}^{*}(s))ds=0$. 2. b)

Next, note that

[TABLE]

The first inequality follows from using, $\mathbf{E}_{\alpha}\nabla f(x_{n}^{*}(\alpha h)))-u\int_{0}^{h}e^{-2(h-s)}\nabla f(x_{n}^{*}(s))ds=0$, and $e^{-2(h-\alpha h)}\leq 1$, and the second inequality follows from Lemma 8.2, and 8.3. 3. c)

For the next part, note that we have

[TABLE]

The third inequality follows from Young’s inequality, the fourth inequality follows from Lemma 8.2, and 8.3, and the fact $1-e^{-2(\alpha-\alpha h)}\leq 2h$, and the fifth inequality follows from (60),and (61). 4. d)

Finally, note that we have

[TABLE]

The first inequality follows from Lemma 2 of [63], and the second inequality follows from Lemma 8.3.

Lemma 8.5

Under conditions of Lemma 8.4,

[TABLE]

Proof. [of Lemma 8.5] Note that, we have

[TABLE]

Lemma 8.6

At iteration $n$, with the initial point $(x_{n},v_{n})$, for the updates (19), (20), and (21), we have

[TABLE]

Proof. [of Lemma 8.6] From Lemma 11 of [63], we have

[TABLE]

From Lemma 11 we also have,

[TABLE]

The second inequality above follows from (62). From Lemma 8.5, we have

[TABLE]

Now, from (68), (69) and Lemma 8.5, we get

[TABLE]

Choosing $h$ such that, $\frac{1}{3}hM\geq Mh^{2}$, i.e., $h\leq\frac{1}{3}$, we get

[TABLE]

Summing both sides from $n=0$ to $N-1$, we get

[TABLE]

Since, $\|v_{0}\|=0$, and $\mathbf{E}\left[f(x_{0})\right]-f(x^{*})\leq O(d)$, and consequently, $\mathbf{E}\left[f(x_{0})-f(x_{N})\right]\leq O(d)$, we have

[TABLE]

Lemma 8.7

At iteration $n$, with the initial point $(x_{n},v_{n})$, for the updates (19), (20), and (21), we have

[TABLE]

Proof. [of Lemma 8.7] From (15) in Lemma 12 of [63] we have,

[TABLE]

From Lemma 12 of [63] we also have,

[TABLE]

Now from (65c), (65d), and Lemma 8.1 we have,

[TABLE]

Combining (70), and (71), we get

[TABLE]

Summing both sides from $n=0$ to $N-1$, and using Lemma 8.6, we get

[TABLE]

Now choosing $\frac{1}{24}uh\geq uh^{2}$, and $v_{0}=0$, we have,

[TABLE]

Using Lemma 8.6, we have,

[TABLE]

Proof. [of Theorem 2.3] From Theorem 3 in [63], we have

[TABLE]

where $q_{N}=\mathbf{E}\left[\|x_{N}-y_{N}\|^{2}+\|x_{N}+v_{N}-y_{N}-w_{N}\|^{2}\right]$. We also have,

[TABLE]

From Lemma 8.4,

[TABLE]

[TABLE]

Combining (74), and (75), we have

[TABLE]

From the proof of Theorem 3 of [63] we have,

[TABLE]

Then we have,

[TABLE]

and

[TABLE]

From (76), (77), and (78), and setting $N=\frac{2\kappa}{h}\log\left(\frac{20}{\epsilon^{2}}\right)$ we have

[TABLE]

From (72), and (73),

[TABLE]

Using, $(h^{3}+\kappa h^{7})\leq 1/2$, we have,

[TABLE]

Choosing, $h=C\min\left(\frac{(\epsilon\sqrt{m})^{\frac{1}{3}}}{(d\kappa)^{\frac{1}{6}}\log\left(\frac{1}{\epsilon}\right)^{\frac{1}{6}}},\min\left(\left(\frac{m}{d}\right)^{\frac{1}{3}},\left(\frac{Mm}{16\sigma^{2}}\right)^{\frac{1}{3}},\sqrt{m}\right)\epsilon^{\frac{2}{3}}\log\left(\frac{1}{\epsilon}\right)^{-\frac{2}{3}}\right)$, we get,

[TABLE]

So, the iteration complexity is given by,

[TABLE]

The total number of zeroth-order oracle calls are given by,

[TABLE]

9 Proofs for Section 3.1

Proof. [of Theorem 3.2] Let us define the following continuous time SDE with the initial point ${\hat{x}}_{0}$:

[TABLE]

Observe that $\hat{x}_{h}$ has the same distribution as $x_{n+1}$ when $\hat{x}_{0}=x_{n}$. Let $z$ denote $(\{u_{i}\}_{i=1}^{b},\{\xi_{i}\}_{i=1}^{b}\})$. To show the dependence of $g_{\nu,b}(\hat{\hat{x}}_{0})$ on $z$ we will use $g_{\nu,b}(\hat{x}_{0},z)$ to denote $g_{\nu,b}(\hat{x}_{0})$ just for this proof. Let $\rho_{t0z}({\hat{x}}_{t},{\hat{x}}_{0},z)$ be the joint distribution of ${\hat{x}}_{t}$, ${\hat{x}}_{0}$, and $z$. Observe that conditioned on ${\hat{x}}_{0}$, and $z$, $g_{\nu,b}(\hat{x}_{0},z)$ is deterministic. Then by Fokker-Plank equation, we have

[TABLE]

Then the time evolution of $\rho_{t}(x)$ is given by

[TABLE]

Now, as shown in [72] we have,

[TABLE]

Then using (79), we have

[TABLE]

Now we use the fact that $\nabla\cdot(ax)=ax\cdot\nabla a+a\nabla\cdot x$ where $a$ is a scalar, and $x$ is a vector:

[TABLE]

Now as $\rho_{t}(x)\left(\nabla\log\left(\frac{\rho_{t}(x)}{\pi(x)}\right)+\mathbf{E}_{0,z|t}\left[g_{\nu,b}({\hat{x}}_{0},z)|{\hat{x}}_{t}=x\right]-\nabla f(x)\right)\log\left(\frac{\rho_{t}(x)}{\pi(x)}\right)$ decays to [math] as $x$ goes to infinity, we have,

[TABLE]

Then we get,

[TABLE]

The second equality above follows from (26), and in the last line we have substituted $x_{t}$ in place of $x$. Now we will upper bound the second term above.

[TABLE]

Now, from Lemma 1.1, we have,

[TABLE]

where $C_{1}=\frac{4(d+5)\sigma^{2}}{b}+\frac{3\nu^{2}M^{2}(d+3)^{3}}{2}$. We also have, with $\tau_{0}\sim N(0,\mathbf{I_{d}})$

[TABLE]

Combining (81), (82), and (83), for $t\leq 1/(2M)$ we get

[TABLE]

We get the last inequality using Lemma 12 of [72]. Now combining, (80), and (84), we get,

[TABLE]

Using (27), we get

[TABLE]

Taking $t\leq h$, we get,

[TABLE]

Multiplying both sides with $e^{\frac{3\lambda t}{2}}$, and integrating from $t=0$ to $h$, we get

[TABLE]

As in [72], in the penultimate step we use the fact $e^{a}\leq 1+2a$ for $0<a=\frac{3\lambda h}{2}$, and $h\leq\frac{2}{3\lambda}$. Hence, we have

[TABLE]

Choosing $b\geq\frac{384M^{2}(d+5)}{\lambda^{2}}$, and $h\leq\frac{\lambda}{12M^{2}}$, we get,

[TABLE]

Then we have,

[TABLE]

Observe that when ${\hat{x}}_{0}=x_{n}$, $\rho_{0}$ is same as $\varpi_{n}$, and then $\rho_{h}(x)$ is same as $\varpi_{n+1}$. Then

[TABLE]

Choosing $n=N=\frac{1}{\lambda h}\log\left(\frac{\epsilon^{2}}{H_{\pi}(\varpi_{0})}\right)$, and using $1-e^{-\lambda h}\geq\frac{\lambda h}{2}$, for $h\leq\frac{1}{\lambda}$, we get

[TABLE]

Choosing $b$, $\nu$, and $h$ as in (29), we get

[TABLE]

Using, (28), we get,

[TABLE]

10 Proofs for Section 4

Proof. [of Lemma 4.1] First note that, we have

[TABLE]

Hence, we have

[TABLE]

Now note that, using independence of $\xi_{i},\xi_{i}^{\prime}$, and $u_{i}$, we have $\forall~{}i$

[TABLE]

We also have, $\forall~{}i\neq j$,

[TABLE]

Using, Lemma 1.1, we hence have,

[TABLE]

Furthermore, we have

[TABLE]

Combining, (87), (86), and (88), we obtain Lemma 4.1.

Proof. [of Theorem 4.2] Using Lemma 4.1, (45) changes to,

[TABLE]

Now the last term involves $\nu$ in the denominator. To counter the effect we have to increase the sample size to $b=\frac{d}{\epsilon^{2}}$.

Proof. [of Lemma 4.2] First note that, we have

[TABLE]

Hence, we have

[TABLE]

Using, Lemma 1.1, we have,

[TABLE]

Furthermore, note that

[TABLE]

Combining, (89), and (90), we get the result stated in Lemma 4.1.

11 Proofs for Section 5

Proof. [of Theorem 5.3] First, we have,

[TABLE]

where $a^{\prime}=a-M\nu\sqrt{s}\leq a-\|\nabla f(\theta)-\nabla f_{\nu}(\theta)\|$ is a lower bound for $|[\nabla f_{\nu}(\theta)]_{j}|$. Next we utilize concentration inequalities to give a bound for the tail of approximation error $\zeta_{j}$. Denote $[g_{\nu,1}]_{j}=\frac{f(\theta+\nu u)-f(\theta)}{\nu}u_{j}\stackrel{{\scriptstyle def}}{{=}}\phi(\nu,u)u_{j}$, where $\phi(\nu,u)$ is sub-exponential with

[TABLE]

where $\|\cdot\|_{\Psi_{1}}=\sup_{p\geq 1}p^{-1}{\bf E}[|\cdot|^{p}]^{1/p}$ and $\|\cdot\|_{\Psi_{2}}=\sup_{p\geq 1}p^{-1/2}{\bf E}[|\cdot|^{p}]^{1/p}$ are the sub-exponential and sub-Gaussian norm respectively (see, for example [70] for more details). In the last inequality we require that $\nu\leq\frac{R}{M\|u\|_{\Psi_{2}}}$. Note that $u\sim N(0,\bm{I}_{d})$ can be replaced by $\sum_{k\in S^{*}}u_{k}e_{k}\sim N(0,\bm{I}_{s})$ due to Assumption 5.1. Moreover, we have the following estimate.

[TABLE]

which implies that $\|\phi(\nu,u)\|_{\Psi_{1}}\leq 2RC_{2}\sqrt{s},\;\|u_{1}\|_{\Psi_{2}}\leq C_{1}$. We now state the following concentration inequality proved in [10].

Lemma 11.1

Let $(X_{i},Y_{i})$, $i=1,\ldots,n$ be $n$ independent copies of random variables $X$ and $Y$. Let $X$ be a sub-Gaussian random variable with $\|X\|_{\psi_{2}}\leq\Upsilon_{1}$, and $Y$ be a sub-exponential random variable with $\|Y\|_{\psi_{1}}\leq\Upsilon_{2}$ for some constants $\Upsilon_{1}$ and $\Upsilon_{2}$. Then for any $t\geq K\cdot\max\{\Upsilon_{1}^{3},\Upsilon_{1}\}\cdot\Upsilon_{2}$, we have

[TABLE]

where $K$ and $K_{1}$ are absolute constants.

From Lemma 11.1, for $n\geq\max\left\{K_{1}\frac{2RC\sqrt{s}}{\tau},\left(\frac{2RC\sqrt{s}}{\tau}\right)^{4}\right\}$, we have:

[TABLE]

where $C=C_{1}C_{2}=\sqrt{\frac{8}{3\log 2(1-\log 2)}},K_{1},K_{2}$ are absolute constants. Therefore, by setting the threshold $\tau=a^{\prime}/2$, the probability of error is bounded by

[TABLE]

Given a pre-specified error rate $\epsilon>0$, it suffices to have $\nu\leq\frac{a}{2M\sqrt{s}}\wedge\frac{R}{MC_{2}\sqrt{s}}$ and

[TABLE]