Sampling from multi-modal distributions with polynomial query complexity in fixed dimension via reverse diffusion

Related work on Gaussian Mixtures -- As noted by the reviewer, the authors of [1] assume that samples of the mixture are available, whereas we assume query access to the potential $V$. While these are two very different settings (Generative modeling vs classical sampling), we will include [1] as a reference in our paragraph on Generative Modeling.

Sampling from a Gaussian Mixture is easy when $V$ is known -- We would like to clarify that we do not assume that $V$ is known analytically, only that we can evaluate $V(x)$ at any given point $x \in \mathbb{R}^d$. In particular, even if $V$ is the potential of a Gaussian Mixture, we do not know the corresponding weights, the means and the covariance matrices.

Result presented before algorithm -- We followed a standard way of presenting results in TCS (theoretical computer science) [1, 5], where the emphasis is put on the consequences (in this case, polynomial complexity for sampling from a multimodal distribution) rather than on the algorithm itself, which is just a means to an end. However, we fully recognize that clarity is important and when introducing a concept early, we will explicitly state in which later section of the paper it will be elaborated.

Comparison with Saremi et al., 2024 -- Even though we use the same estimator of the score, note that we do not use it within the same pipeline: in Saremi et al., 2024 the scores are used in a chain of MCMC at a single noise level. In our work, we use this estimator within the reverse diffusion framework, with multiple noise levels. Hence, our overall algorithms differ. Furthermore, note that unlike Saremi et al., 2024, we provide a non-asymptotic error for this estimator (Proposition 4) and we fully control the overall error of our algorithm which is, in our opinion, the main contribution of our work.

Third contribution is an overclaim -- We are unsure of the link between the SNIS estimator of the score and our "third contribution": in our introduction, we stated our third contribution as the fact that our algorithm does not require prior knowledge on the target to provably work. Is this what you had in mind?

Dissipativity beyond Gaussian Mixtures -- First, we emphasize that our work being the first one to handle Gaussian Mixtures, this seemingly toy case is not an easy one. Furthermore, dissipativity is a quite broad setting: any distribution that is strongly log-concave outside a compact is dissipative and so are all distributions of the form $p(x) \propto e^{- \| x \|^m}$ with $m \geq 2$ (and by stability of dissipativity by mixtures, so are all mixtures of the two previous kinds).

Comparison with vanilla score matching -- The asymptotic sample complexity of vanilla score matching scales with the square of the Poincaré constant of the target density (Theorem 2 in [3]). For multi-modal distributions such as Gaussian mixtures, this constant increases exponentially with the between-mode distance (Section 4.1 in [4]). In contrast, our result shows polynomial scaling with respect to the between-mode distance (see our Section 3), which we believe is a key theoretical strength of our work.

We hope we addressed all of your questions and hope you will consider revising your score. Please let us know if any further clarifications are needed.

[1] Chen et al., Learning general Gaussian mixtures with efficient score matching, Arxiv 2024

[2] Saremi et al. Chain of log-concave markov chains, ICLR 2024

[3] Koehler et al., Statistical Efficiency of Score Matching: The View from Isoperimetry, ICLR 2023

[4] Schitling, Poincare and Log–Sobolev Inequalities for Mixtures, Entropy 2019

[5] Jason Altschuler and Sinho Chewi, Faster high-accuracy log-concave sampling via algorithmic warm starts, Journal of the ACM 2024.

Simple combination of previous works --- Indeed, our algorithm simply combines reverse diffusion + a self normalized estimator of the score. However, that combination having never been explored before, that makes it a new algorithm in its own right. Furthermore, we must insist that our contribution is not the algorithm in itself, which is not very hard to come up with, but rather to study its error and most importantly quantitatively prove that, in terms of query complexity, it is state-of-the-art compared to previous works as it is robust to multi-modality.

Additional references --- we have read through the additional references and agree that they are relevant. We appreciate the reviewer's clear summary and DSI / TSI classification. We will add these references and discuss them, especially as they use estimators of the scores which is a central part of our analysis.

Toy experiment --- The Figure 1. should not be seen as an experiment but rather as a visual illustration of the performance of our algorithm in comparison with previous works. We highlight the fact that related works have been published at this venue and similar venues without experimental results [1,2,3,4], as the primary contribution is theoretical. That being said, we ran an experiment that was suggested by the reviewer. The distribution to sample is a 2D mixture of three Gaussians with equal weights $1/3$, equidistant modes located at a distance $R$ from the origin, and different covariances $(\mathrm{Id}, \mathrm{Id}/2, \mathrm{Id}/4)$. We allowed each 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. Finally, we generated $500$ samples for each three algorithms that we compared in Wasserstein distance to $500$ samples generated from ground truth. As expected, as the between-mode distance grows with $R$, our algorithm yields the lowest error.

Doubts on poor performance of RDMC on 2D illustration -- It is not surprising that RDMC fails on this 2D experiment as the dimension is not the relevant parameter for RDMC, but rather the isoperimetric constants (Log-Sobolev and/or Poincaré) of the target. Indeed, as we explicitly state in Section 3, RDMC uses ULA to sample from $q_{t, z}(x) \propto e^{-V(x)} e^{- \| x e^{-t} - z \|^2/(2 (1-e^{-2t}))}$ which is almost as hard as sampling from the target $e^{-V}$ as $t$ grows. In particular, in our example, since the target is highly multi-modal, we expect this inner step to fail when $t = O(1)$ and as a consequence, the intermediate scores $\nabla \log(p_t)$ will be poorly estimated.

Initialization of Langevin for RDMC --- For RDMC, the inner Langevin step is initialized with a standard gaussian $\mathcal{N}(0, I_2)$. While some improvements are possible [5, 6], they may require unknown parameters of the target distribution (e.g. its smoothness) and the task of sampling from a multi-modal distribution with an inner ULA step remains fundamentally difficult.

Why our illustration is in 2D --- Our theory shows that time-reversed diffusion, when combined with a suitable SNIS score estimator, achieves polynomial query complexity only in low dimensions. This motivates our choice of a 2D illustration. We emphasize that the target distribution is not a "simple 2D problem": in the absence of prior knowledge of the mode locations, sampling is highly nontrivial. In fact, even in 1D, sampling from this unnormalized density using ULA, for instance, typically incurs exponential query complexity in its parameters (e.g. smoothness), whereas our method retains polynomial complexity.

Ability to recover mode weights --- as is typical in theoretical work, our guarantees are stated in terms of distributional closeness in KL divergence. Explicitly assessing the recovery of mode weights would require a different statistical divergence that is specifically sensitive to mode weights. This is unfortunately beyond the scope of our analysis as well as most of the theoretical literature, to the best of our knowledge. Nevertheless, this is an interesting question for future research.

Positioning our contribution --- given their success, time-reversed diffusions have recently gained traction for sampling from unnormalized densities. The reviewer lists many algorithms that are "simple combinations of previous works" yet are still noteworthy. Our method fits this category but 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. We hope our responses have addressed the reviewer’s concerns and that they will consider revising their score.

[1] Chen et al., Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions, ICLR 2023.

[2] Lee et al., Convergence for score-based generative modeling with polynomial complexity, NeurIPS 2022.

[3] Conforti et al., KL Convergence Guarantees for Score Diffusion Models under Minimal Data Assumptions, SIAM 2025.

[4] Chen et al., Learning general Gaussian mixtures with efficient score matching, Arxiv 2024

[5] Grenioux et al., Stochastic Localization via Iterative Posterior Sampling, ICML 2024

[6] Huang et al., Faster Sampling without Isoperimetry via Diffusion-based Monte Carlo, COLT 2024

We thank the reviewer for his careful review and positive comments on our paper.

Clarifications in the text --- thank you for picking up on these. We will clarify the caption of Figure 1 with the following explanation. The color scheme indicates the probability density value of the distribution we want to sample from (dark is low probability density, bright is high probability). The blue dots are the samples produced by our algorithm. We will also regroup notations in a paragraph of table as suggested by the reviewer.

Experimental validation --- our paper aims to provide strong theoretical guarantees for sampling using time-reversed diffusions. Related works have been published at this venue and similar venues without experimental results [1,2,3,4], as the primary contribution is theoretical. Figure 1 is intended as a visual aid to support the theory. That being said, we ran an experiment that was suggested by reviewer 29WY. The distribution to sample is a 2D mixture of three Gaussians with equal weights $1/3$, equidistant modes located at a distance $R$ from the origin, and different covariances $(\mathrm{Id}, \mathrm{Id}/2, \mathrm{Id}/4)$. We compare convergence in Wasserstein distance between ULA, RDMC, and our algorithm. As expected, as the between-mode distance grows with $R$, our algorithm yields much lower error than ULA and RDMC.

Removing the constraint on the dimension --- As was stated in our introduction, the exponential dependence in the dimension cannot be avoided under our hypotheses because of existing lower-bounds [5,6]. Intuitively, this is because finding the modes of the target distribution would require exploring the entire space, and that search problem grows exponentially with the dimensionality. Our results surprisingly show that, unlike many sampling algorithms, the exponential complexity is only in the dimension, not in the other parameters (e.g. smoothness) of the distribution.

[1] Chen et al., Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions, ICLR 2023.

[2] Lee et al., Convergence for score-based generative modeling with polynomial complexity, NeurIPS 2022.

[3] Conforti et al., KL Convergence Guarantees for Score Diffusion Models under Minimal Data Assumptions, SIAM 2025.

[4] Chen et al., Learning general Gaussian mixtures with efficient score matching, Arxiv 2024

[5] Lee et al., Beyond Log-concavity: Provable Guarantees for Sampling Multi-modal Distributions using Simulated Tempering Langevin Monte Carlo, NeurIPS 2018

[6] Chak, On theoretical guarantees and a blessing of dimensionality for nonconvex sampling, Arxiv 2024

Contributions and future work --- We appreciate the reviewer’s positive assessment. Time-reversed diffusions have significantly impacted sample generation from multi-modal distributions, and we believe this work offers a timely contribution by providing a strong theoretical foundation for their use to sampling from unnormalized densities. In response to the reviewer’s comment, interesting future directions include developing theoretical guarantees for multi-modal distributions with non-Gaussian tails.