1 Introduction↩︎

Sampling from a distribution whose density is only known up to a normalization constant is a fundamental problem in statistics. Formally, given some potential \(V: \mathbb{R}^d \to \mathbb{R}\) such that \(\int e^{-V(x)} \mathrm{d}x < \infty\), the sampling problem consists in obtaining a sample from some distribution \(p\) such that \(p\) is \(\epsilon\)-close to the target \(\mu \propto e^{-V}\) with respect to some divergence while maintaining the complexity, i.e., the number of queries to \(V\) and possibly to its derivatives, as low as possible. Depending on the shape of the distribution, the typical complexity of existing sampling algorithms can significantly differ.

1.0.0.1 Log-concave and "PL-like" distributions

As in Euclidean optimization, a common assumption in the sampling literature is to assume that \(\mu\) is log-concave and log-smooth or equivalently, that \(V\) is convex and smooth. Specifically, when the potential \(V\) is assumed to be \(\alpha\)-strongly convex and to have an \(L\)-Lipschitz gradient, the popular Unajusted Langevin Algorithm (ULA) is known to achieve fast convergence [1], [2]. Because strong-log concavity implies uni-modality, thus excluding many distributions of interest, this assumption was further relaxed to \(\mu\) verifying an \({\alpha}^{-1}\)-log-Sobolev inequality, later interpreted as a Polyak–Łojasiewicz type condition on \(\text{KL}(\cdot|\mu)\) for the Wasserstein geometry [3]. Under these conditions, ULA was shown to achieve \(\epsilon\)-error in Kullback-Leibler (\(\mathop{\mathrm{KL}}\nolimits\)) divergence in \(\tilde{O}(L^2 \alpha^{-2} d \epsilon^{-1})\) queries to \(\nabla V\) [4]. While these polynomial guarantees do go beyond the uni-modal setting, we show in the next paragraph that most existing algorithms still fail to sample from truly multi-modal distributions.

1.0.0.2 Multi-modal distributions

Designing sampling algorithms for multi-modal distributions is an active area of research. However, most existing sampling algorithms are limited in at least two of the following ways:

1. The query complexity is exponential in the parameters of the problem. For instance, when the global \(\alpha\)-strong convexity of \(V\) is relaxed with \(\alpha\)-strong convexity outside a ball of radius \(R\), the log-Sobolev constant of \(\mu\) degrades to \(O(e^{8LR^2}\alpha)\) [5]. In practice, this does not simply translate to poor worst case bounds: practitioners are well aware that when dealing with multi-modal distributions, ULA-based algorithms suffer from metastability, where they get stuck in local modes, leading to slow convergence [6].

2. Guarantees are obtained under a log-smoothness assumption; we show in Sec. 3.2 that even for Gaussian Mixtures, which is arguably the most basic multi-modal model, this assumption may not be verified.

3. Explicit a priori knowledge on the target distribution is required for the algorithm to converge. For instance, importance sampling requires a proposal distribution whose "effective" support must cover sufficiently well the one of the target; however, since \(V\) can only be evaluated point-wise, the access to such a support is unclear. In fact, this limitation goes beyond the multi-modal setting in some cases: no matter the log-Sobolev constant of the target, similarly to euclidean gradient descent, ULA still requires an a priori bound on the inverse smoothness constant of \(V\) and sufficiently small step size in comparison in order to converge.

Generally speaking, finding a sampling algorithm that addresses even the first limitation is not possible. Recent results provide lower-bounds on the complexity of sampling from a multi-modal distribution: they are exponential in the dimension, as shown in [7], [8] and [9]. However, when the dimensionality is fixed, the question of whether the three limitations can be addressed remains. In fixed dimension, can we sample from a broad class of multi-modal distributions with a polynomial number of queries in the parameters of the problem, (e.g. \(L, R, \alpha\)) and without prior knowledge on these parameters?

1.0.0.3 Our contributions

Our work answers this question positively. We provide a sampling algorithm that addresses the three limitations outlined above. First, we show that our algorithm has polynomial complexity in all the problem parameters but the dimension, thus enabling to efficiently sample from highly multi-modal distributions in fixed dimension. Second, we show that unlike most existing results, our guarantees hold under relaxed regularity assumptions that cover general Gaussian mixtures. Third, our algorithm yield guarantees without prior knowledge on the parameters of the distribution.

The workhorse of our algorithm is made of two key ingredients: the reverse diffusion scheme, that transfers the sampling problem into a score estimation problem with different levels of noise, and self-normalized importance sampling, to estimate these noisy scores using only query access to \(V\) and additional Gaussian samples. Both reverse diffusion [10]–[12] and self-normalized importance sampling [13]–[17] were already used separately to sample in a Bayesian setting, where only the log-density of the target is available. However, these works either offer no theoretical guarantees or fail to improve upon existing results. In this article, we combine both reverse diffusion and self-normalized importance sampling in a single algorithm, allowing us to recover the first polynomial time guarantees to sample from highly multi-modal distributions in fixed dimension.

This paper is structured as follows. First, we survey related work on multi-modal sampling in 2 and present our main result in 3. Then, we detail our sampling algorithm in 4 and provide the key ingredients for the proof in 5.

2 Related work↩︎

Before detailing our sampling algorithm in Section 4, we review the main alternatives to ULA when sampling from a multi-modal distribution and the guarantees they offer. We show that existing approaches suffer from at least two of three drawbacks outlined in the introduction: exponential query complexity, a restrictive smoothness assumption that excludes general Gaussian mixtures (as shown later in Sec. 3.2), and unavailable prior knowledge of the distribution.

2.0.0.1 Proposal-based algorithms

At a high level, proposal-based algorithms use a proposal distribution that is easy to sample from and that can either be directly used as a proxy for the target, or, alternatively, whose samples will be rejected (respectively re-weighted) to obtain approximate samples from the target; the corresponding algorithm for the latter case is the well-known rejection sampling (respectively importance sampling) scheme. Guarantees for these methods typically assume that the target distribution is log-smooth. Furthermore, they must use carefully designed proposals that require prior knowledge of the target distribution that is often unavailable in practice.

For instance, when the target is \(L\)-log-smooth and with finite second moment \(m_2\), one can design and sample from a proxy distribution that is \(\epsilon\)-close in \(\mathop{\mathrm{TV}}\nolimits\) to the target using \(\tilde{O}((L m_2\epsilon^{-1})^{O(d)})\) queries to the potential function \(V\) and its gradient \(\nabla V\) [9]. While this bound is indeed polynomial in fixed dimension, the design of the proxy requires an \(\epsilon\)-approximation of the global minimum of \(V\) which can only be achieved if \(m_2\), therefore the location of the mass, is explicitly available; this is rarely the case in practice.

Similarly, if we further assume that the target is \(\alpha\)-strongly log-concave outside of a ball of radius \(R\), one can achieve an \(\epsilon\)-precise approximation of the target distribution in \(\mathop{\mathrm{TV}}\nolimits\) with polynomial number of queries to \(\nabla V\) when \(d\) is fixed via an importance sampling scheme [8]. In this case, the proposal is such that it coincides with the target when the potential of the target is above some cutoff, and is flat elsewhere. While there exists a cutoff value ensuring the log-concavity of the proposal, thus allowing its efficient sampling, this value depends on the unknown constants \(L, R\) and \(\alpha\) [8].

2.0.0.2 Tempering-based algorithms

Instead of directly sampling from \(\mu \propto e^{-V}\), these algorithms start by sampling from a flattened version of \(\mu\) given by \(\mu^\beta \propto e^{-\beta V}\) for small \(\beta\), and gradually increase \(\beta\) to \(1\), a strategy sometimes also referred to as annealing. When \(\mu\) is assumed to be a finite mixture of the same shifted \(\alpha\)-strongly log-concave and \(L\)-log-smooth distribution, [7] proved that for a well-chosen (stochastic) sequence of flattened distributions \(\mu^{\beta_i}\), sampling up to precision \(\epsilon\) in \(\mathop{\mathrm{KL}}\nolimits\) can be achieved in \(poly(L, \alpha^{-1}, w_{min}^{-1}, d, \epsilon^{-1}, R)\) queries to \(\nabla V\), where \(R\) is the location of the furthest mode and \(w_{\min}\) the minimum weight in the mixture. However, this setting is quite restrictive: for Gaussian mixtures for instance, it only handles the case where the covariance matrices are identical for all components. Furthermore, the algorithm requires explicit knowledge of \(R\) which is unavailable.

2.0.0.3 Föllmer Flows

Instead of directly sampling from the target distribution \(\mu \propto e^{-V}\), Föllmer Flows start from a simple (e.g. Gaussian or Dirac mass) distribution that is progressively interpolated to the target using a Schrödinger bridge. When the initial distribution is a Dirac mass at \(0\), this bridge solves a closed-from SDE [18] which can thus be discretized to generate samples from the target. In the context of Bayesian inference, where one only has access to the unnormalized density, [19] estimated the drift of the resulting SDE via neural methods; in particular no guarantees on the sampling quality are provided. In [13]–[16], the drift is estimated via a Monte-Carlo method. While in appearance, these works provide strong polynomial guarantees, a closer look shows that these guarantees only hold if the function \(f(x) = e^{-V(x) + \| x\|^2/2}\) is Lipschitz, smooth and bounded from below; we show in Appendix 9.1.2 that this assumption is quite restrictive.

2.0.0.4 Diffusion-based methods

Over the past few years, diffusion-based algorithms, and especially reverse diffusion, that we shall review in details in Sec. 3, have emerged as solid candidates for multi-modal sampling. In essence, they allow to transfer the sampling problem into the problem of estimating the scores of intermediate distributions that are given by the convolution of the initial distribution with increasing levels of Gaussian noise. Under an \(\epsilon\)-oracle of these intermediate scores, it has been shown that diffusion-based methods could yield an \(\epsilon\)-approximate sample of the target in \(\text{poly}(\epsilon^{-1})\) time under milder and milder assumptions [20]–[26]. This framework has been applied with tremendous empirical success in generative modeling, where numerous samples of the target are already available and one seeks to produce new samples from the target, for several years now [27]–[30]. However, it was only quite recently that this framework has been applied in a Bayesian context, where only an unnormalized density of the target, instead of samples, is available [10]–[12], [31], [32]. In a work closely related to ours, [11] showed that if the intermediate scores remain \(L\)-log-smooth for any noise level, which in particular implies that the target itself is \(L\)-log-smooth, their algorithm could reach a complexity of \(O(e^{L^3\log^3((Ld+m_2)/\epsilon)})\) with \(m_2\) the second order moment of the target. In [12], the smoothness assumption is relaxed with a sub-quadratic growth assumption \(V(x) - V(x^*) \leq L \| x - x^*\|^2\), where \(x^*\) is any global minimizer of \(V\), yet at the price of an oracle access to \(V(x^*)\); under these assumptions, the authors manage to obtain a complexity that is at best \(O(L^{d/2} \epsilon^{-d} e^{L \| x^*\|^2 + \| x_N \|^2})\) where \(x_N\) is the final output sample. Under the reasonable (and desirable) assumption that \(\mathbb{E}[\| x_N||^2] \approx m_2\), Jensen’s inequality yields an overall \(O(L^{d/2} \epsilon^{-d} e^{L \| x^*\|^2 + m_2})\) complexity. In particular, both these works suffer at least two of the limitations mentioned in the introduction, making them ill-suited for multi-modal sampling.

3 Presentation of the main result and application to Gaussian Mixtures↩︎

3.1 Main result↩︎

As in the works mentioned above, we rely on the recent advances on reverse diffusion and focus on the task of estimating the intermediate scores. Using an estimator that is described in Sec. 4, we recover polynomial sampling guarantees for densities verifying the assumptions described hereafter.

This assumption shall be referred to as semi-log-convexity, by analogy with the functional analysis literature [33]. Note that is has also been referred to as 1-sided Lipschitzness for \(\nabla V\) [25]. It relaxes the classical log-smoothness assumption which implies the additional lower bound \(\nabla^2 V \succeq - \beta I_d\). In particular, unlike the latter, a mixture of semi-log-convex densities remains semi-log-convex [34]; we provide a quantitative version of this statement in Sec. 3.2 in the case of Gaussian mixtures.

This assumption is referred to as dissipativity as common in the sampling and optimization literature [35]–[37]. Note that this assumption relaxes strong convexity outside of a ball which can be equivalently re-written as \(\langle \nabla V(x) - \nabla V(y), x- y \rangle \geq a \| x - y \|^2 - b\) for all pairs \((x,y) \in \mathbb{R}^d\), also referred to as strong dissipativity [38], [39]. We show in Sec. 3.2 that unlike strong-dissipativity, a mixture of dissipative distributions remain dissipative. distribution is also log-concave outside a ball.

Our algorithm addresses the three limitations outlined in the sections above. When the dimension \(d\) is fixed, we obtain a polynomial query complexity. This guarantee does not assume log-smoothness and applies to general Gaussian mixtures, as will be shown in the next subsection. Moreover, this guarantee does not require running the algorithm with any explicit knowledge of the target distribution’s constants \(a, b, \beta\).

3.1.0.1 Overall comparison

We summarize in Table 1 how our algorithm compares to previous approaches in terms of assumptions, oracles required to run the algorithm and resulting complexity. Along with the work of [9], our algorithm is the only one that is polynomial in the parameters of the distribution when \(d\) is fixed. Furthermore, while our dissipativity assumption is stronger than finite second moment, we relax the log-smoothness assumption by semi-log-convexity which notably covers general Gaussian mixtures. Finally, as mentioned above, because their algorithm requires an \(\epsilon\)-approximation of the global minimum of \(V\), they require an explicit upper-bound on the second order moment of \(\mu\) which may not be available in practice.

3.1.0.2 Numerical illustration

In Figure 2, we consider a standard task in the literature: sampling from a mixture of 16 equally weighted Gaussians with unit variance and centers uniformly distributed in \([-40, 40]^2\) [40]. We compare our algorithm against Unadjusted Langevin Algorithm (ULA) and to the reverse diffusion algorithm of [10] (RDMC). We also implemented the zeroth-order method of [12], but it failed to converge. ULA was initialized from \(\mathcal{N}(0, I_2)\) and run for \(5 \times 10^4\) steps. All three methods used the same discretization step size \(h = 0.01\). Our reverse diffusion algorithm and RDMC were both run with \(500\) reverse diffusion steps and both used \(100\) samples to estimate the intermediate scores. As discussed in Sec. 4, while our samples are simply drawn from a Gaussian distribution, the samples used in RDMC are drawn from a auxiliary, multi-modal distribution generated via an inner ULA step. In this experiment, this inner ULA was initialzed with a standard Gaussian and was run for \(100\) steps; in particular, we emphasize that while ULA and our algorithm were roughly given the same computational budget, the one given to RDMC was a hundred times superior. As a result, ULA and our algorithm took approximately one minute to run on a computer locally and the RDMC method required over an hour. We observe that unlike the two others, our algorithm successfully recovers all the modes.

In Figure 1, we monitor the convergence of the same three algorithms as a function of the problem difficulty, measured by the distance between the modes. Here, the target distribution is a mixture of three Gaussians in two dimensions. It has equal weights \(1/3\), equidistant modes located at a distance \(R\) from the origin, and different covariances \((I, I/2, I/4)\). The final error is measured in Wasserstein distance, because it can be easily approximated using samples: we use \(500\) samples generated by the sampling algorithm and \(500\) samples from the true target distribution. We allowed each of the three algorithms \(10^6\) queries to the potential \(V\) or to its gradient: we performed \(10^6\) iterations for ULA, we used \(100\) steps of reverse diffusion for our algorithm and RDMC. We used \(100\) inner steps of ULA for RDMC that was initialized with a standard Gaussian. We used \(100\) particles for score estimation for RDMC and \(10000\) for our algorithm. As expected, as the between-mode distance grows with \(R\), our algorithm yields the lowest error.

3.2 Application to Gaussian mixtures↩︎

We now apply our results to derive provable sampling guarantees for general Gaussian mixtures. Apart from the very recent exception of [41] that we discuss below, note that despite their wide popularity, no sampling guarantees were yet derived for general Gaussian mixtures. Indeed, unless they verify some specific assumptions such as identical covariance matrices among components [26], general Gaussian mixtures do not satisfy the log-smoothness assumption, thus they do not fit the framework of many previously discussed works. In fact, their gradient may not even be Hölder continuous: consider the simple counter-example of a two-dimensional mixture \(\mu = 0.5 \mathcal{N}(0, \Sigma_1) + 0.5 \mathcal{N}(0, \Sigma_2)\) with covariances \(\Sigma_1 = diag(1, 0.5)\) and \(\Sigma_2 = diag(0.5, 1)\). On the diagonal \(x=y\), the score is \(\nabla \log(\mu)(x, x) = -3/2(x, x)\), while near the diagonal, right above it for instance, the score behaves asymptotically as \(\nabla \log(\mu)(x,x) \sim_{ x \to +\infty} -(2x, x)\); we provide a rigorous analysis in Appendix 9.2. Fortunately though, Gaussian mixtures do verify Assumptions 1-2, as we next show.

The proof is deferred to Appendix 9.3. Combined with Theorem 1, Proposition 2 shows that we can sample any Gaussian mixture with average precision \(\epsilon\) in \(\mathop{\mathrm{KL}}\nolimits\) in \(O(\mathrm{poly}(\kappa,R, d, \lambda_{min}^{-d}, \epsilon^{-d}))\) queries to \(V\) where \(R = \max_{i}\| \mu_i \|\) and \(\kappa = \lambda_{max}/\lambda_{min}\). There exists a relatively recent literature seeking to relax the log-smoothness assumption. For instance,  [37], [42], [43] work under weak smoothness, i.e. \(\alpha\)-Hölder continuous gradient of the potential for \(\alpha\) in \([0,1]\) (recall that \(\alpha\)=1 recovers the smooth case). However, as shown above, this relaxation is not sufficient yet to cover general Gaussian mixtures. The only reference we know of that does is the work of [41], who used a regularized version of Langevin to relax the global smoothness condition with local Lipschitz smoothness and polynomial growth of the Lipschitz constant. Yet crucially, their complexity bound scales as a polynomial of the Poincaré constant of \(\mu\). We prove in Appendix 9.4 that for the mixture \(\mu = 0.5 \mathcal{N}(Ru, \lambda_{\max} I_d) + 0.5 \mathcal{N}(-Ru, \lambda_{\min} I_d)\) with \(u\) any unit vector and \(0< \lambda_{\min} \leq \lambda_{\max}\), this constant is at least \((R^2e^{R^2/(2\lambda_{\max})})/2\). In particular, this method still degrades exponentially with the multi-modality parameter \(R\).

4 Our sampling algorithm↩︎

In this section, we introduce the reverse diffusion framework and explain how it reduces the sampling problem to that of estimating the scores along the forward Ornstein-Uhlenbeck (OU) process.

4.1 Reverse diffusion: from sampling to score estimation↩︎

Reverse diffusion methods emerged as an alternative to Langevin-based samplers in order to overcome metastability and were first introduced to the ML community in [44]. They rely on the so-called forward process \[\label{eq:def95ou} \begin{cases} & \mathrm{d}X_t = - X_t \mathrm{d}t+ \sqrt{2} \mathrm{d}B_t \,\\ & X_0 \sim \mu \, , \end{cases}\tag{1}\] which corresponds to the standard OU process initialized at \(\mu\), that is a specific case of a Langevin diffusion targeting a standard Gaussian, that we will denote \(\pi\). Note that since the target of this process, the standard Gaussian, is \(1\)-strongly log-concave, the resulting process converges exponentially fast to the equilibrium. In order to sample from \(\mu\), reverse diffusion algorithms rely on the semi-discretized backward process: given a horizon \(T\) that we discretize as \(0=t_0 \leq t_1 \leq \cdots t_{N-1} \leq t_N = T\), the latter writes \[\label{eq:def95backward} \mathrm{d}Y_t = Y_t \mathrm{d}t + 2\nabla \log(p_{t_{N-k}})(Y_t) \mathrm{d}t + \sqrt{2}\mathrm{d}B_t \, , ~~ t \in ]T- t_{N-k}, T- t_{N-(k+1)}] \, ,\tag{2}\] with \(Y_0 \sim p_T\) and where \(p_{t}\) is the distribution of the forward process Eq. 1 at time \(t\). Note that this reverse process cannot be readily implemented for two reasons: first, it requires the knowledge of the intermediate scores \(\nabla \log(p_{t_k})\) which are not available in closed form. Second, it requires sampling from the distribution \(p_T\). Nevertheless, if one can access a proxy \(s_{t_k}\) of the scores \(\nabla \log(p_{t_k})\), and considering \(T\) large enough so that \(p_T\approx\pi\), we can implement instead \[\label{eq:def95backward95approx} \mathrm{d}Y_t = Y_t \mathrm{d}t + 2s_{t_k}(Y_k) \mathrm{d}t + \sqrt{2}\mathrm{d}B_t \, , ~~ t \in ]T- t_{N-k}, T- t_{N-(k+1)}] \, ,\tag{3}\] with \(Y_0 \sim \pi\) and where all iterations can be solved in closed form. Because the forward process Eq. 1 converges exponentially fast, we can expect the initialization error \(Y_0 \sim \pi\) instead of \(Y_0 \sim p_T\) to be small after a short time \(T\). Furthermore, if the proxies \(s_{t_k}\) are sufficiently accurate, one can expect that the process output by the approximate scheme Eq. 3 has a distribution that is close to the target \(\mu\). Over the past three years, several works provided quantitative bounds of the error induced by the discretization, the use of an approximate score and the initialization error with respect to different divergences and under various assumptions [20], [21], [24], [29], [45]. Yet, we shall rely exclusively on the following theorem as it is the most suited to our framework.

The previous theorem shows that under mild assumptions that notably allow for multi-modality, the problem of sampling from \(\mu\) can be transferred into a score approximation problem along the forward process. In the next subsection, we present an estimator for these intermediate scores that is tractable given the knowledge of the unnormalized density \(\mu\propto e^{-V}\).

4.2 Derivation of an estimator↩︎

The key observation to derive an estimator of the intermediate scores is that the forward process Eq. 1 is nearly available in closed form. However, this closed form may be written in different manners.

4.2.0.1 Different expressions of the scores

Consider Eq. 1 integrates to \(X_t = \sqrt{\lambda_t} X_0 + \sqrt{1-\lambda_t} Z \, ,\) where \(\lambda_t = e^{-2t}\), \(Z\) is a standard \(d\)-dimensional Gaussian, and \(X_0\) is a random variable simulating the target distribution. The corresponding density convolves the target density with a Gaussian, as \[\begin{align} p_t^V(x) &= \frac{1}{\sqrt{1-\lambda_t}} \pi \left( \frac{x}{\sqrt{1 - \lambda_t}} \right) \ast \frac{1}{\sqrt{\lambda_t}} \mu \left( \frac{x}{\sqrt{\lambda_t}} \right) \, , \end{align}\] where \(\pi\) is a standard Gaussian and \(\mu(\cdot) \propto \exp(-V(\cdot))\) is the target distribution. From this expression, we can obtain different formulas for the scores. We retain

1. \(\nabla \log p_t^V(x) = \frac{1}{\sqrt{\lambda_t}} \mathbb{E}_{y | x} \left[ \nabla \log \mu \left( \frac{x - y}{\sqrt{\lambda_t}} \right) \right] , \quad y | x \sim \propto \mu \left( \frac{x - y}{\sqrt{\lambda_t}} \right) \times \mathcal{N} ( y; 0, (1 - \lambda_t) I )\)   ,

2. \(\nabla \log p_t^V(x) = \frac{1}{\sqrt{\lambda_t}} \frac{ \mathbb{E}_{y} \left[ \nabla \log \mu \left( \frac{x - \sqrt{1 - \lambda_t} y}{\sqrt{\lambda_t}} \right) \mu \left( \frac{x - \sqrt{1 - \lambda_t} y}{\sqrt{\lambda_t}} \right) \right] }{ \mathbb{E}_{y} [ \mu \left( \frac{x - \sqrt{1 - \lambda_t} y}{\sqrt{\lambda_t}} \right) ] } , \quad y \sim \mathcal{N}(y; 0, I)\)   ,

3. \(\nabla \log p_t^V(x) = \frac{1}{1 - \lambda_t} \left( \mathbb{E}_{y | x} [y] - x \right) , \quad y | x \sim \mu \left( \frac{y}{\sqrt{\lambda_t}} \right) \times \mathcal{N}(y; x, (1 - \lambda_t) I)\)   ,

4. \(\nabla \log p_t^V(x) = \frac{1}{\sqrt{1 - \lambda_t}} \frac{ \mathbb{E}_{y} \left[ y \mu \left( \frac{x + \sqrt{1 - \lambda_t} y}{\sqrt{\lambda_t}} \right) \right] }{ \mathbb{E}_{y} \left[ \mu \left( \frac{x + \sqrt{1 - \lambda_t} y}{\sqrt{\lambda_t}} \right) \right] } , \quad y \sim \mathcal{N} ( y; 0, I )\, .\)

Derivations can be found in Appendix 8. The second [32], third [10], [12], [31], and fourth [17] identities have been used to build Monte Carlo estimators of the score. Hence, one natural way to group these identities is by how difficult it is to draw the samples \(y\). The second and fourth identities sample from a Gaussian distribution, whereas the first and third identities sample from an auxiliary multi-modal distribution from generating \(y|x\). As \(t\) gets closer to zero, this auxiliary distribution resembles the original target, so sampling from it progressively becomes as difficult as sampling from the original target density. This observation explains why, unsurprisingly, [10], [12] do not improve upon existing results. Another way to group these identities is whether or not they require evaluating the score of the target distribution: identities one and two do and are referred to as TSI estimators, referring to the Target Score Identity (TSI) used to obtain them [46]; identities three and four do not and are referred to as DSI estimators, referring to a Denoising Score Identity (DSI) used to obtain them [47]. Combinations of TSI and DSI estimators have also been considered [47], [48].

4.2.0.2 A self-normalized estimator

We now focus on the estimator of the scores that we use. It is a rewrite of the fourth identity and is a ratio of expectations under Gaussians \[\nabla \log(p_t^V)(z) = \frac{-1}{1-e^{-2t}}\frac{\mathbb{E}[Y_te^{-V(e^{t}(z-Y_t))}]}{\mathbb{E}[e^{-V(e^{t}(z-Y_t))}]} \, ,\] with \(Y_t \sim \mathcal{N}(0, (1-e^{-2t})I_d)\) which is easy to simulate. While conventional statistical wisdom may suggest using independent samples to estimate both the numerator and the denominator, we voluntarily choose to correlate them and implement instead

\[\label{eq:def95estimator} \hat{s}_{t,n}(z) \mathrel{\vcenter{:}}= \frac{-1}{1-e^{-2t}} \frac{\sum_{i=1}^n y_ie^{-V(e^{t}(z-y_i))}}{\sum_{i=1}^n e^{-V(e^{t}(z-y_i))}} \, ,\tag{4}\] where the \(y_i\) are independent Gaussians such that \(y_i \sim \mathcal{N}(0, (1-e^{1-2t})I_d)\) ; we refer to this estimator as self-normalized as common in the sampling literature [49]. The key property of self-normalized estimators is that they remain nearly bounded: in our case, it holds uniformly in \(z\) that \[\mathbb{E}[\| \hat{s}_{t,n}(z)] \| \leq \frac{\mathbb{E}[\max_i \| y_i \|]}{1-e^{-2t}} \sim \sqrt{\frac{d\log(n)}{1-e^{-2t}}} \, .\] This boundedness will allow us to derive a non-asymptotical control on the quadratic error that we present in the next section.

We now explain how this differs from previous work. While this estimator was already considered in [17], the authors did not use it within the reverse diffusion pipeline and more importantly, we are the first work to derive quantitative guarantees for this estimator, which is one of our core contributions. Similarly, an estimator close to this one was considered in the context of Föllmer Flows [14]–[16] in order to approximate the shift in the corresponding Schrödinger bridge. We show in Appendix 9.1.1 that while this shift is itself close to the intermediate scores that we seek to approximate, the resulting self-normalized estimator is more degenerate. Namely, as discussed above, their guarantees are significantly weaker.

5 Sketch of proof of the main result↩︎

The proof is decomposed in three steps: (i) we derive a non-asymptotic bound on the quadratic error of the estimator presented in Eq. 4 (ii) we show that under Assumptions 1-2, the integrated error of this estimator (that appears in the bound of 3) can be fully controlled by the zeroth and second order moments of the ratio \(\Phi_t = p_t^{2V}/p_t^{V}\), where \(p_t^{2V}\) (respectively \(p_t^{V}\)) is the density of the forward process defined in Eq. 1 initialized at \(\mu^2\) (respectively \(\mu\)) (iii) we provide a quantitative bound on these moments as well as other relevant quantities and we conclude with Theorem 3.

The complete proof is left to Appendix 10 yet we briefly sketch the main arguments. We split the expectation on the event \(A\) where the empirical denominator \(\hat{D}\) of Eq. 4 (respectively the empirical numerator \(\hat{N}\)) is not too small (respectively not too large) with respect to its expectation \(D\) (respectively \(\| N \|\)) and on the complementary \(\bar{A}\). Over \(A\), we use a second-order Taylor expansion to make the variances of both the numerator and the denominator appear, and compute them explicitly. Conversely, the quadratic error of the estimator remains almost bounded on \(\bar{A}\). We use Chebyshev’s inequality to upper-bound \(\mathbb{P}(\bar{A})\) and make the variances of the numerator and of the denominator appear again, which concludes the proof.

Then, we need to control the Laplacian of the forward processes as well as their gradient appearing in \(\theta_t\) in the former proposition. As mentioned in the previous section, the intermediate scores can be re-written as \[\nabla \log(p_t^V)(z) = \frac{z - e^{-t} \mathbb{E}_{q_{t,z}}[Y]}{1-e^{-2t}} \, ,\] with \(q_{t,z}(x) \propto e^{-V(x)} e^{-\frac{\| e^{-t} x - z \|^2}{2(1-e^{-2t})}}\). Now we note that if \(\mu \propto e^{-V}\) is dissipative, so is \(q_{t,z}\); in particular we can quantitatively bound its second order moment and a fortiori, upper-bound \(\|\nabla \log(p_t^V)(z)\|^2\).

The complete proof is left in Appendix 11.2. This result implies that the \(\theta_t\) term defined in Proposition 4 is at most of order of order \(\theta_t(z)\sim (1+ \| z \|^2)\) w.r.t. \(z\). In particular, the average integrated error \(\mathbb{E}[\| \hat{s}_{t,n}(z) - \nabla \log(p_t^V) \|_{L^2(p_t^V)}^2]\) can be upper-bounded with respect to the zeroth and the second order moments of the ratio \(\Phi_t = p_t^{2V}/p_t^{V}\). In the next proposition, we show that these moments are bounded under the semi-log-convexity assumption. By a slight abuse of notation, we shall denote \(m_i(\Phi_t) = \int \| x \|^{i} \Phi_t(z) \mathrm{d}z\) the \(i\)-th moment of \(\Phi_t\).

The proof of this lemma is deferred to Appendix 12. It relies on a key result in [33] where it is shown that for \(\beta\)-semi-log convex distributions, the following bound holds: \[\nabla^2 \log(\pi) - \nabla^2 \log(p_t^V) \preceq \frac{(\beta-1) e^{-2t}}{(1-e^{-2t})(\beta-1) + 1}I_d \, .\] In particular, the right-hand-side remains bounded w.r.t. \(\beta\) whenever \(t>0\), which allows to avoid an exponential dependence in \(\beta\) in our final bounds. It suffices now to control \(\mathcal{I}(\mu, \pi), m_2\) and \(\mu(0)\) in order to apply Theorem 3 and conclude.

The proof is left to Appendix 11.1. We can now state our main result.

The proof is deferred to Appendix 13 and is mainly an application of the results collected above.

6 Conclusion↩︎

In this article, we successfully applied the reduction from sampling to intermediate score estimation, initiated over the past three years by [20]–[22], [24], to the problem of low dimensional multi-modal sampling. Using the self-normalized estimator of the scores, our results provide polynomial query complexity guarantees in fixed dimension, apply to general Gaussian mixtures, and do not require prior knowledge of the target distribution’s constants. Interesting future directions include extending theoretical guarantees to more general multi-modal distributions with heavy tails for instance.

We note that our sampling algorithm is based on time-reversed diffusions, which have recently gained traction for sampling from unnormalized densities. Our method stands out by offering rigorous, non-asymptotic theoretical guarantees on query complexity, especially in the presence of score estimation error that we precisely quantify. Such theoretical guarantees are scarce and we believe our results are therefore a meaningful and timely contribution to the field.

7 Acknowledgements↩︎

This work was supported by the Agence Nationale de la Recherche (ANR) through the JCJC WOS project and the PEPR PDE-AI project (ANR-23-PEIA-0004). We thank Pierre Monmarché for his help in proving Lemma 3.

8 Expressions the intermediate scores↩︎

We next provide derivations of the Monte Carlo estimators of the intermediate scores discussed in section 4.2.

The probability law along the reverse diffusion path has density \[\begin{align} p_t^V(x) &= \frac{1}{\sqrt{1-\lambda_t}} \pi \left( \frac{x}{\sqrt{1 - \lambda_t}} \right) \ast \frac{1}{\sqrt{\lambda_t}} \mu \left( \frac{x}{\sqrt{\lambda_t}} \right) \end{align}\] where \(\pi\) is a standard Gaussian and \(\mu(\cdot) \propto \exp(-V(\cdot))\) is the target distribution

8.0.0.1 Case 1: write the convolution as an integral against the proposal distribution

We have \[\begin{align} p_t^V(x) &= \frac{1}{Z_{\lambda_t}} \int_{y \in \mathbb{R}^d} \mu \left( \frac{x - y}{\sqrt{\lambda_t}} \right) \pi \left( \frac{y}{\sqrt{1 - \lambda_t}} \right) dy = \frac{1}{Z_{\lambda_t}} \int f(x, y) dy \end{align}\] where we denoted the integrand by \(f(x, y)\). We can compute the score

\[\begin{align} \nabla \log p_t^V(x) &= \frac{\nabla p_t^V(x)}{p_t^V(x)} = \frac{\int \nabla f(x, y) dy}{\int f(x, y) dy} = \frac{1}{\sqrt{\lambda_t}} \int \nabla \log \mu \left( \frac{x - y}{\sqrt{\lambda_t}} \right) \frac{f(x, y)}{\int f(x, y) dy} dy \enspace . \end{align}\] From this, we can either define a Monte Carlo estimator as \[\begin{align} \nabla \log p_t^V(x) &= \frac{1}{\sqrt{\lambda_t}} \mathbb{E}_{y | x} \left[ \nabla \log \mu \left( \frac{x - y}{\sqrt{\lambda_t}} \right) \right] , \quad y | x \sim \propto f(x, y) \end{align}\]

where \(f(x, y)\) is a smoothened version of the target distribution distribution. Replacing \(\pi\) with a standard Gaussian. Or else, by unpacking \(f(x, y)\) and using the proposal as the sampling distribution as the integration variable \(y\) appears in it, \[\begin{align} \nabla \log p_t^V(x) &= \frac{1}{\sqrt{\lambda_t}} \frac{ \mathbb{E}_{y} \left[ \nabla \log \mu \left( \frac{x - y}{\sqrt{\lambda_t}} \right) \mu \left( \frac{x - y}{\sqrt{\lambda_t}} \right) \right] }{ \mathbb{E}_{y} [ \mu \left( \frac{x - y}{\sqrt{\lambda_t}} \right) ] } , \quad y \sim \pi \left( \frac{y}{\sqrt{1 - \lambda_t}} \right) \end{align}\]

8.0.0.2 Case 2: write the convolution as an integral against the target distribution

We now have \[\begin{align} p_t^V(x) &= \frac{1}{Z_{\lambda_t}} \int_{y \in \mathbb{R}^d} \mu \left( \frac{y}{\sqrt{\lambda_t}} \right) \pi \left( \frac{x - y}{\sqrt{1 - \lambda_t}} \right) dy = \frac{1}{Z_{\lambda_t}} \int f(x, y) dy \end{align}\] where we now denote the integrand by \(f(x, y)\). We can compute the score \[\begin{align} \nabla \log p_t^V(x) &= \frac{\nabla p_t^V(x)}{p_t^V(x)} = \frac{\int \nabla f(x, y) dy}{\int f(x, y) dy} \\ &= \frac{1}{\sqrt{1 - \lambda_t}} \int \nabla \log \pi \left( \frac{x - y}{\sqrt{1 - \lambda_t}} \right) \frac{f(x, y)}{\int f(x, y) dy} dy \enspace . \end{align}\] From this, we can either define a Monte Carlo estimator as \[\begin{align} \nabla \log p_t^V(x) &= \frac{1}{\sqrt{1 - \lambda_t}} \mathbb{E}_{y | x} \left[ \nabla \log \pi \left( \frac{x - y}{\sqrt{1 - \lambda_t}} \right) \right] , \quad y | x \sim \propto f(x, y) \end{align}\] where \(f(x, y)\) is a smoothened version of the target distribution distribution. Or else, by upacking \(f(x, y)\) and using the target as the sampling distribution as the integration variable \(y\) appears in it, \[\begin{align} \nabla \log p_t^V(x) &= \frac{1}{\sqrt{1 - \lambda_t}} \frac{ \mathbb{E}_{y} \left[ \nabla \log \pi \left( \frac{x - y}{\sqrt{1 - \lambda_t}} \right) \pi \left( \frac{x - y}{\sqrt{1 - \lambda_t}} \right) \right] }{ \mathbb{E}_{y} [ \pi \left( \frac{x - y}{\sqrt{1 - \lambda_t}} \right) ] } , \quad y \sim \mu \left( \frac{y}{\sqrt{\lambda_t}} \right) \end{align}\] which is not very useful given that sampling the target is hard in the first place. However, we can use the proposal distribution as the sampling distribution as the integration variable \(y\) appears in it as well. To see this, recall that \[\begin{align} \pi \left( \frac{x - y}{\sqrt{1 - \lambda_t}} \right) &= \mathcal{N} \left( \frac{x - y}{\sqrt{1 - \lambda_t}} ; 0, I_d \right) = \sqrt{1 - \lambda_t} \mathcal{N} \left( x - y ; 0, (1 - \lambda_t) I_d \right) \\ &= \sqrt{1 - \lambda_t} \mathcal{N} \left( y ; x, (1 - \lambda_t) I_d \right) \end{align}\] This leads to \[\begin{align} \nabla \log p_t^V(x) &= \frac{1}{1 - \lambda_t} \frac{ \mathbb{E}_{y} \left[ (y - x) \mu \left( \frac{y}{\sqrt{\lambda_t}} \right) \right] }{ \mathbb{E}_{y} \left[ \mu \left( \frac{y}{\sqrt{\lambda_t}} \right) \right] } , \quad y \sim \mathcal{N} ( y; x, (1 - \lambda_t) I_d ) \\ &= \frac{1}{\sqrt{1 - \lambda_t}} \frac{ \mathbb{E}_{z} \left[ z \mu \left( \frac{x + \sqrt{1 - \lambda_t} z}{\sqrt{\lambda_t}} \right) \right] }{ \mathbb{E}_{z} \left[ \mu \left( \frac{x + \sqrt{1 - \lambda_t} z}{\sqrt{\lambda_t}} \right) \right] } , \quad z \sim \mathcal{N} ( z; 0, I_d ). \end{align}\] Using a change of variables, we obtain \[\begin{align} \nabla \log(p_t^V)(x) = \frac{-1}{1-\lambda_t}\frac{ \mathbb{E} \big[ Y_t \exp(-V( \frac{x - Y_t}{\sqrt{\lambda_t}})) \big] }{ \mathbb{E}\big[ \exp(-V( \frac{x - Y_t}{\sqrt{\lambda_t}})) \big] } \end{align}\] with \(Y_t \sim \mathcal{N}(0, (1-\lambda_t)I_d)\).

9 Additional discussions↩︎

9.1 Föllmer Flows↩︎

9.1.1 Comparison with our work↩︎

Instead of implementing a reverse diffusion process, Föllmer Flows [13]–[15] seek to implement the following Schrödinger bridge: \[\mathrm{d}X_t = b(X_t, t) \mathrm{d}t + \mathrm{d}W_t, ~~ X_0 = 0 \, ,\] with \(b(x, t) = \nabla \log \mathbb{E}[f(x + W_{1-t})]\) and where \(f\) is given by \[f \mathrel{\vcenter{:}}= \frac{\mathrm{d}\mu}{\mathrm{d}\mathcal{N}(0, I_d)} \, ,\] with \(\mu\) the target density; for such a process, we have that \(X_1 \sim \mu\). We refer the reader to the notations of Proposition 11 to observe that the shift \(b(x,t)\) is given by \[\begin{align} b(x,t) & = \nabla \log Q_{-\log(t)/2}(x/\sqrt{t}) \, ,\\ & = \frac{1}{\sqrt{t}}(\nabla \log p_{-\log(t)/2}^V(x/\sqrt{t}) - \nabla \log \pi(x)) \, , \end{align}\] with \(\pi\) is the density of the standard Gaussian. Hence, up to the Gaussian term and a time-rescaling, the shift is exactly given by the intermediate scores of the reverse diffusion. However, because of the Gaussian correction, an extra exponential term appears in the associated self-normalized estimator of the shift: in [13]–[15], the authors implement \[\hat{b}(x, t) = \frac{\sum_{i=1}^n (x + \sqrt{1-t}z_i - \nabla V(x + \sqrt{1-t}z_i))e^{- V(x + \sqrt{1-t}z_i) + \| x + \sqrt{1-t}z_i\|^2/2}}{\sum_{i=1}^ne^{- V(x + \sqrt{1-t}z_i) + \| x + \sqrt{1-t}z_i\|^2/2}} \, ,\] with \(z_i\) \(n\)-i.i.d. samples from the standard Gaussian distribution. We suspect that because of the extra exponential terms \(e^{\| x + \sqrt{1-t}z_i\|^2/2}\), as discussed below, their estimator provide appealing guarantees only under very stringent assumptions.

9.1.2 Theoretical guarantees↩︎

The theoretical complexity bounds derived in the works of [13]–[15] quantitatively rely on the assumptions that \(f =\frac{\mathrm{d}\mu}{\mathrm{d}\mathcal{N}(0, I_d)}\) is Lipschitz, has Lipschitz gradient and is bounded from below. In particular, assume for instance that \(\mu\) is a standard Gaussian centered at some point \(c\in \mathbb{R}^d\). In this case, the ratio \(f\) reads \[f(x) = e^{2 x^\top c + \| c \|^2} \, ,\] which verifies none of the assumptions above if \(c \neq 0\). More broadly, even if \(f\) verifies these assumptions, the resulting quantities degrade exponentially with the mismatch between \(x \mapsto-\log(\mu)(x)\) and \(x \mapsto \| x \|^2/2\). As noted in [19], this limitation is not only a theoretical artifact; in practice, these methods are very unstable and fail to convergence on simple examples.

9.2 Details on the non-smooth example↩︎

Consider the mixture \(\mu = 0.5 \mathcal{N}(0, \Sigma_1) + 0.5 \mathcal{N}(0, \Sigma_2)\) with \(\Sigma_1 = diag(1, 0.5)\) and \(\Sigma_1 = diag(0.5, 1)\). For \((x,y) \in \mathbb{R}^2\), it holds that \[- \nabla \log(\mu)(x,y) = \frac{(x, 2y) e^{-x^2/2 - y^2} + (2x, y) e^{-x^2 - y^2/2}}{e^{-x^2/2 - y^2} + e^{-x^2 - y^2/2}} \, .\] Hence, when \(x=y\), we have \(- \nabla \log(\mu)(x,y) = 3/2(x,x)\). Now for \(y = x + \eta\), the score reads \[\begin{align} - \nabla \log(\mu)(x,x + \eta) & = \frac{(x, 2x+ 2\eta) e^{-x^2/2 - x^2 - 2\eta x - \eta^2} + (2x, x + \eta) e^{-x^2 - x^2/2 - \eta x - \eta^2/2}}{e^{-x^2/2 - x^2 - 2\eta x - \eta^2} + e^{-x^2 - x^2/2 - \eta x - \eta^2/2}} \\ & = \frac{(x, 2x+ 2\eta) e^{ - \eta x - \eta^2/2} + (2x, x + \eta) }{e^{ - \eta x - \eta^2/2}+1} \\ & \sim_{x\to + \infty} (2x, x) \, . \end{align}\] In particular, \(- \nabla \log(\mu)\) is not Hölder.

9.3 Proof of Proposition 2↩︎

Recall we denote \(\lambda_{I_dn} := I_dn_i \lambda_{I_dn}(\Sigma_i), \quad \lambda_{\max} := \max_i \lambda_{\max}(\Sigma_i)\). We write: \[\mu(x) = \sum_{i=1}^p \tilde{w}_i \exp\left( -\frac{1}{2} (x - \mu_i)^\top \Sigma_i^{-1} (x - \mu_i) \right):=\sum_{i=1}^p \tilde{w}_i \phi_i(x)\] with \(\tilde{w}_i = w_i (2\pi)^{-d/2} \det(\Sigma_i)^{-1/2}\).

9.3.0.1 Bound on the Hessian.

The Hessian of \(\mu\) writes: \[\nabla^2 \log \mu(x) = \frac{1}{\mu(x)} \nabla^2 \mu(x) - \frac{1}{\mu(x)^2} \nabla \mu(x) \nabla \mu(x)^\top\] Since \[\begin{align} \nabla \phi_i(x)&= -\Sigma_i^{-1}(x - \mu_i) \, \phi_i(x)\\ \nabla^2 \phi_i(x)& = \left[ \Sigma_i^{-1} (x - \mu_i)(x - \mu_i)^\top \Sigma_i^{-1} - \Sigma_i^{-1} \right] \phi_i(x) \end{align}\] We have \[\begin{align} \nabla \mu(x) &=\sum_i \tilde{w}_i \nabla \phi_i(x) = -\sum_i \tilde{w}_i \phi_i(x) \Sigma_i^{-1}(x - \mu_i)\\ \nabla^2 \mu(x) &= \sum_{i=1}^p \tilde{w}_i \nabla^2 \phi_i(x) = \sum_{i=1}^p \tilde{w}_i \phi_i(x) \left[ \Sigma_i^{-1} (x - \mu_i)(x - \mu_i)^\top \Sigma_i^{-1} - \Sigma_i^{-1} \right] \end{align}\] Denoting \(\gamma_i(x) := \frac{\tilde{w}_i \phi_i(x)}{\mu(x)} \quad\) and \(s_i(x)=- \Sigma_i^{-1} (x - \mu_i)\) we get \[\begin{align} \nabla^2 \log \mu(x) &= \sum_i \gamma_i(x) \left[ s_i(x) s_i(x)^\top - \Sigma_i^{-1} \right] - \left( \sum_i \gamma_i(x) s_i(x) \right) \left( \sum_j \gamma_j(x) s_j(x) \right)^\top\\ & = \operatorname{Cov}_{\gamma(x)}[s_i(x)] - \sum_i \gamma_i(x) \Sigma_i^{-1} \end{align}\] where the first term is a covariance matrix of the vectors \(s_i(x)\) under the weights \(\gamma_i(x)\), hence which is a positive semi-definite matrix. Therefore: \[\begin{align} -\nabla^2 \log \mu(x) &\preceq \sum_i \gamma_i(x) \Sigma_i^{-1}\preceq \sum_i \gamma_i(x) \cdot \frac{I_d.}{\lambda_{I_dn}} = \frac{I_d.}{\lambda_{I_dn}} . \end{align}\]

9.3.0.2 Bound on the drift.

For the (negative) score of the mixture we have: \[-\nabla \log \mu(x) = \sum_i \gamma_i(x) \Sigma_i^{-1}(x - \mu_i)\] Hence \[\begin{align} \langle- \nabla \log \mu(x), x \rangle &= \sum_{i=1}^p \gamma_i(x) \langle \Sigma_i^{-1}(x - \mu_i), x \rangle\\ &= \sum_{i=1}^p \gamma_i(x) \left( x^\top \Sigma_i^{-1} x - \mu_i^\top \Sigma_i^{-1} x \right) \ge \sum_{i=1}^p \gamma_i(x) \left( \frac{\|x\|^2}{\lambda_{\max}} - \mu_i^\top \Sigma_i^{-1} x \right), \end{align}\] where we have used \(\Sigma_i^{-1} \succeq \frac{1}{\lambda_{\max}} I_d\) for the first term. Now for the second term, using Cauchy-Schwartz and Young’s inequality: \[\begin{align} \mu_i^\top \Sigma_i^{-1} x \leq \|\mu_i\| \cdot \|\Sigma_i^{-1} x\| \leq \frac{\|\mu_i\| \|x\|}{\lambda_{I_dn}} \leq \frac{1}{2\lambda_{I_dn}} (\|x\|^2 + \|\mu_i\|^2). \end{align}\] Hence, using that \(\sum_i \gamma_i(x) = 1\) and that \(\sum_i \gamma_i(x) \|\mu_i\|^2 \leq \max_i \|\mu_i\|^2\) we have: \[\begin{align} \langle - \nabla \log \mu(x), x \rangle &\geq \sum_{i=1}^p \gamma_i(x) \left( \frac{\|x\|^2}{\lambda_{\max}} - \frac{1}{2\lambda_{I_dn}}(\|x\|^2 + \|\mu_i\|^2) \right) \\ &= \left( \frac{1}{\lambda_{\max}} - \frac{1}{2\lambda_{I_dn}} \right) \|x\|^2 - \frac{1}{2\lambda_{I_dn}} \sum_{i=1}^p \gamma_i(x) \|\mu_i\|^2.\\ &\geq \left( \frac{1}{\lambda_{\max}} - \frac{1}{2\lambda_{I_dn}} \right) \|x\|^2 - \frac{1}{2\lambda_{I_dn}} \max_i \|\mu_i\|^2 \\& \geq \frac{\|x\|^2}{2\lambda_{\max}} - \lambda_{\max} \cdot \max_i \left( \frac{\|\mu_i\|}{\lambda_{I_dn}} \right)^2 . \end{align}\]

9.4 Discussion on  [41]↩︎

To justify our claims at the end of 3.2, we first need some preliminary background on Poincaré constants. Let \(q \in \mathcal{P}_{\rm ac}(\mathbb{R}^d)\). We say that \(q\) satisfies the Poincaré inequality with constant \(C_P\ge0\) if for all \(f \in W_0^{1,2}(q)\) (functions with zero-mean with respect to \(q\) and whose gradient is squared integrable with respect to \(q\)): \[\label{eqn:poincare} \int f^2(x)dq(x) \leqslant C_P\|\nabla f\|^2_{L^2(q)},\tag{5}\] and let \(C_P(q)\) be the best constant in Eq. 5 , or \(+\infty\) if it does not exist. We show that in the case where \(\mu = 1/2 \mathcal{N}(c, \lambda_{max}) + 1/2 \mathcal{N}(-c, \lambda_{min})\) the Poincaré constant is at least \(2e^{\frac{\| c \|^2}{2\lambda_{max}}}/\|c\|^2\) therefore the resulting sampling complexity in [41] is at least \(O(\text{poly}(\|c\|^2e^{\frac{\|c\|^2}{2\lambda_{max}}}, \epsilon^{-1}))\).

9.5 Discussion on [49]↩︎

In [49], for \(\nu\) dominated by \(\pi\), denoting \(g\) the unnormalised density ratio \[\frac{\mathrm{d}\nu}{\mathrm{d}\pi} (u) \mathrel{\vcenter{:}}= \frac{g(u)}{\int g(u) \mathrm{d}\pi(u)} \, ,\] the authors derived error bounds for self-normalised estimators: for any test function \(\phi: \mathbb{R}^d \to \mathbb{R}\), the expectation of \(\phi\) under \(\nu\), denoted \(\nu(\phi)\), is estimated via \[\nu^n(\phi) = \frac{\sum_{i=1}^n \phi(x_i) g(x_i)}{\sum_{i=1}^n g(x_i)} \, ,\] with \((x_i)\) \(n\)-iid samples drawn from \(\pi\). In Theorem 2.3, they provide the following quadratic error bound: \[\label{eq:ub95agapiou} \mathbb{E}[(\nu(\phi) - \nu^n(\phi))^2] \lesssim \frac{3\pi(|\phi g|^{2d})^{1/d}m_{2e}(g)^{1/e}}{n\pi(g)^4} \, ,\tag{6}\] with \(m_t(h) = \pi(|h(\cdot) - \pi(h)|^t)\) and \(d,e \in (1, +\infty[\) such that \(1/d+1/e = 1\). Now recall that, denoting \(Y_t \sim \mathcal{N}(0, (1-e^{-2t})I_d)\) the score along the forward is given by \[\begin{align} \nabla \log(p_t^V)(z) &= \frac{-1}{1-e^{-2t}} \frac{\mathbb{E}[Y_t e^{-V(e^{t}(z-Y_t))}]}{\mathbb{E}[e^{-V(e^{t}(z-Y_t))}]} \, ,\\ & = \frac{-1}{1-e^{-2t}}\frac{\int y e^{-V(e^{t}(z-y))} e^{-\frac{\|y\|^2}{2(1-e^{-2t})}}\mathrm{d}y}{\int e^{-V(e^{t}(z-y))} e^{-\frac{\|y\|^2}{2(1-e^{-2t})}} \mathrm{d}y} \, ,\\ & = \frac{-1}{1-e^{-2t}} \mathbb{E}_{Z\sim\nu}[Z]\, ,\\ \end{align}\] where \(\nu(x) \propto e^{-V(e^{t}(z-x))} e^{-\frac{\|x\|^2}{2(1-e^{-2t})}}\). In particular, denoting \(\pi(x) \propto e^{-\frac{\|x\|^2}{2(1-e^{-2t})}}\), the quadratic error of the self-normalized estimator reads \[\begin{align} \|\hat{s}_{t,n}(z) - \nabla \log p_t^V(z)\|^2 & = \frac{1}{(1-e^{-2t})^2} \biggl \| \frac{\sum_{i=1}^n y_i g(y_i)}{\sum_{i=1}^n g(y_i)} - \mathbb{E}_{Z\sim \nu}[Z] \biggr \|^2 \, , \\ & = \frac{1}{(1-e^{-2t})^2} \sum_{k=1}^d \biggl ( \frac{\sum_{i=1}^n \phi_k(y_i) g(y_i)}{\sum_{i=1}^n g(y_i)} - \mathbb{E}_{Z\sim \nu}[\phi_k(Z)] \biggr)^2 \, . \end{align}\] with \(g = \mathrm{d}\nu / \mathrm{d}\pi\), \((y_i)\) \(n\)-iid samples drawn from \(\pi\) and \(\phi_k: \mathbb{R}^d \to \mathbb{R}\) the projector on the \(k\)-th coordinate. Hence, the bound Eq. 6 of [49] provides the following upper-bound \[\|\hat{s}_{t,n}(z) - \nabla \log p_t^V(z)\|^2 \leq \sum_{k=1}^d \frac{3\pi(|\phi_k g|^{2d})^{1/d}m_{2e}(g)^{1/e}}{n\pi(g)^4} \, .\] We show in the next paragraph that, when integrated against \(p_t^V\), this upper-bound is vacuous. By Jensen’s inequality, it holds that \[\begin{cases} & m_{2e}(g)^{1/(2e)} \geq m_{2}(g)^{1/2} \, ,\\ & \pi(|\phi_k g|^{2d})^{1/(2d)} \geq \pi(|\phi_k g|^{2})^{1/2} \, . \end{cases}\] Hence, the previous right-hand side of is lower-bounded as \[\sum_{k=1}^d \frac{3\pi(|\phi_k g|^{2d})^{1/d}m_{2e}(g)^{1/e}}{n\pi(g)^4} \geq \sum_{k=1}^d \frac{3\pi(|\phi_k g|^{2})m_{2}(g)}{n\pi(g)^4} = \frac{3m_{2}(g)}{n\pi(g)^4} \sum_{k=1}^d \pi(|\phi_k g|^{2})\, .\] We now compute explicitly each of the quantities in the right-hand-side above using Proposition 9. First we have \[\begin{align} \pi(g) & =(2\pi(1-e^{-2t}))^{-d/2} \int e^{-V(e^{t}(z-x))} e^{-\frac{\|x\|^2}{2(1-e^{-2t})}} \mathrm{d}x\, , \\ & = e^{-td}Z_Vp_t^{V}(z) \, . \end{align}\] Then, the variance term reads \[\begin{align} m_{2}(g) & = (2\pi(1-e^{-2t}))^{-d/2}\int (e^{-V(e^{t}(z-x))} - e^{-td}Z_Vp_t^V(z))^2 e^{-\frac{\| x\|^2}{2(1-e^{-2t})}} \mathrm{d}x \, ,\\ & = e^{-td}Z_{2V}p_t^{2V}(z) - 2e^{-2td}Z_{2V}Z_Vp_t^V(z)p_t^{2V}(z) + (e^{-td}Z_V)^2p_t^V(z)^2 \, . \end{align}\] Finally, the second order moment term reads \[\begin{align} \sum_{k=1}^d \pi(|\phi_k g|^{2}) & = (2\pi(1-e^{-2t}))^{-d/2}\int \| x \|^2 e^{-2V(e^{t}(z-x))} e^{-\frac{\| x\|^2}{2(1-e^{-2t})}} \mathrm{d}x \, ,\\ & = e^{-td}p_t^{2V}(z)\Theta(z) \, , \end{align}\] where we denoted \[\Theta(z) \mathrel{\vcenter{:}}= Z_{2V}(1-e^{-2t})^2(\Delta \log(p_t^{2V})(z) + \| \nabla \log(p_t^{2V})(z) \|^2 + d/(1-e^{-2t}))\, .\] Hence, once integrated against \(p_t^V\), our lower bound reads \[\begin{align} \int \frac{3m_{2}(g)}{n\pi(g)^4} \sum_{k=1}^d \pi(|\phi_k g|^{2})p_t^V(z) \mathrm{d}z = & ~e^{2td}\frac{Z_{2V}^2}{Z_V^4}\int \Theta(z)\frac{p_t^{2V}(z)^2}{p_t^V(z)^3} \mathrm{d}z - 2e^{td}\frac{Z_{2V}^2}{Z_V^3}\int \Theta(z)\frac{p_t^{2V}(z)^2}{p_t^V(z)^2} \mathrm{d}z \\ & + e^{td}\frac{Z_{2V}}{Z_V^2}\int \Theta(z)\frac{p_t^{2V}(z)}{p_t^V(z)} \mathrm{d}z \, . \end{align}\]

Combining Proposition 6 and Lemma 1, we have that the third term is always finite. For the specific case \(V(x) = \frac{\| x \|^2}{2}\), we show that the second term also remains finite while the first diverges thus making the overall lower bound diverge. First, we recall that for this choice of potential, we have \[\begin{cases} & p_t^V(z) \propto e^{-\frac{\|z\|^2}{2}} \, , \\ & p_t^{2V}(z) \propto e^{-\frac{\|z\|^2}{2(1-e^{-2t}/2)}} \, , \end{cases}\] from what we obtain \[\begin{cases} & \Theta(z) = Z_{2V}\big(\| z\|^2 + \frac{d(1-e^{-2t})[2(1-e^{-2t})+e^{-4t}]}{2-e^{-2t}}\big) \, , \\ & \frac{p_t^{2V}(z)^2}{p_t^{V}(z)^2} \propto e^{-\frac{\|e^{-t}z\|^2}{2(1-e^{-2t}/2)}} \, , \\ & \frac{p_t^{2V}(z)^2}{p_t^{V}(z)^3} \propto e^{\frac{\|z\|^2(1-e^{-2t})}{2(1-e^{-2t}/2)}} \, . \end{cases}\] In particular, we obtain that the second term is bounded and that the first term diverges to infinity making the overall bound of [49] vacuous.

10 Proof of Proposition 4↩︎

Before starting the proof, we recall the following identities.

For the rest of the proof we shall drop the dependence in \(z\) and write \(\hat{s}_{t,n}(z) = \frac{\hat{N}}{\hat{D}}\) with the empirical numerator \(\hat{N} = -\frac{\sum_{i=1}^n y_i e^{-V(e^t(z-y_i))}}{1-e^{-2t}}\) and denominator \(\hat{D} = \sum_{i=1}^n e^{-V(e^t(z-y_i))}\) where we recall \(y_i \sim \mathcal{N}(0, (1-e^{-2t})I_d)\). In the following proposition, we explicitly compute the variances of \(\hat{N}\) and \(\hat{D}\) using the formulas above. In what follows, we shall denote \(N = \mathbb{E}[N]\) and \(D = \mathbb{E}[D]\)

We can now prove Proposition 4.

11 Proof of Proposition  6 and Lemma 2↩︎

Before starting the proofs, we recall the following usefull lemma that bounds the second order moment of dissipative distributions.

11.1 Proof of Lemma 2↩︎

11.2 Proof of Proposition  6↩︎

Recall that the intermediate scores read \[\nabla \log(p_t^V)(z) = \frac{z - e^{-t} \mathbb{E}_{q_{t,z}}[Y]}{1-e^{-2t}} \, ,\] with \(q_{t,z}(x) \propto e^{-V(x)} e^{-\frac{\| e^{-t} x - z \|^2}{2(1-e^{-2t})}}\). In particular, if \(V\) is dissipative with constant \(a, b\) then it holds that \[\begin{align} \langle -\nabla \log(q_{t,z})(x), x \rangle &= \langle \nabla V(x) + \frac{e^{-t}(e^{-t}x -z)}{1-e^{-2t}}, x \rangle \, , \\ & \geq a \| x \|^2 - b + \frac{e^{-t}}{1-e^{-2t}}(e^{-t} \|x\|^2 - \langle z, x \rangle). \end{align}\] Now recall that \(e^{-t} \|x\|^2 - \langle z, x \rangle \geq -e^{t} \| z \|^2/4\) which yields \[\langle -\nabla \log(q_{t,z})(x), x \rangle \geq a \| x\|^2 - (b+\|z \|^2/(4(1-e^{-2t})))\, .\] Hence, using Lemma 3, it holds that \[\mathbb{E}_{q_{t,z}}[\| Y \|^2] \leq \frac{2b+ \|z \|^2/(2(1-e^{-2t})) + d}{a} \, .\] Hence we recover that \[\begin{align} \| \nabla \log(p_t^V)(z) \|^2 & \leq \frac{2\|z \|^2}{(1-e^{-2t})^2} + \frac{2e^{-2t} \mathbb{E}_{q_{t,z}}[\| Y \|^2]}{(1-e^{-2t})^2} \, ,\\ & \leq \frac{2\|z \|^2}{(1-e^{-2t})^2} + \frac{e^{-2t}}{(1-e^{-2t})^2}\biggl(\frac{2(2b+d)}{a} + \frac{\| z \|^2}{(1-e^{-2t})a}\biggr) \, , \\ & = \frac{\|z \|^2}{(1-e^{-2t})^2}\biggl( 2 + \frac{e^{-2t}}{a(1-e^{-2t})}\biggr) + \frac{2e^{-2t}(2b +d)}{a(1-e^{-2t})^2} \, . \end{align}\] Similarly, recall that \[\nabla^2 \log(p_t^V)(z) - \nabla^2 \log(\pi)(z) = \frac{e^{-2t}}{1-e^{-2t}} \left(\frac{\text{Cov}_{q_{t, z}}(X)}{1-e^{-2t}} - I_d \right) \, ,\] from which we can deduce \(\Delta \log(p_t^{2V}) \leq \frac{e^{-2t}}{(1-e^{-2t})^2} \mathbb{E}_{q_{t,z}}[\| Y \|^2]\) which yields again \[\Delta \log(p_t^{2V}) \leq \frac{e^{-2t} }{(1-e^{-2t})^2} \biggl( \frac{\|z \|^2}{2a(1-e^{-2t})} + \frac{2b + d}{a}\biggr) \, .\]

12 Proof of Lemma 1↩︎

Before starting the proof, we recall the result of [33].

Before proving Lemma 1, we introduce this preliminary result on the evolution of \(\Phi_t\).

We now provide the proof of Lemma 1.

In order to recover a bound on the second moment of \(\Phi_t\), we need several intermediate results. We first prove that the maximum of the ratio decreases through time.

We then derive an upper-bound on the maximum of a log-smooth distribution.

Using the previous result, we can derive an upper-bound on the integrated squared gradient at \(0\).

We can now bound the first order moment of \(\Phi_t\)

We can now derive our upper-bound on \(m_2(\Phi_t)\).

13 Proof of Theorem 7↩︎

References↩︎