On the Mode-Seeking Properties of Langevin Dynamics

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 11:06:49

Abstract: The Langevin Dynamics framework, which aims to generate samples from the score function of a probability distribution, is widely used for analyzing and interpreting score-based generative modeling. While the convergence behavior of Langevin Dynamics under unimodal distributions has been extensively studied in the literature, in practice the data distribution could consist of multiple distinct modes. In this work, we investigate Langevin Dynamics in producing samples from multimodal distributions and theoretically study its mode-seeking properties. We prove that under a variety of sub-Gaussian mixtures, Langevin Dynamics is unlikely to find all mixture components within a sub-exponential number of steps in the data dimension. To reduce the mode-seeking tendencies of Langevin Dynamics, we propose \emph{Chained Langevin Dynamics}, which divides the data vector into patches of constant size and generates every patch sequentially conditioned on the previous patches. We perform a theoretical analysis of Chained Langevin Dynamics by reducing it to sampling from a constant-dimensional distribution. We present the results of several numerical experiments on synthetic and real image datasets, supporting our theoretical results on the iteration complexities of sample generation from mixture distributions using the chained and vanilla Langevin Dynamics.

Main paper: 
 introduction: 
 A central task in unsupervised learning involves learning the underlying probability distribution of training data and efficiently generating new samples from the distribution. Score-based generative modeling (SGM) (Song et al., 2020c) has achieved state-of-the-art performance in various learning tasks including image generation (Song and Ermon, 2019, 2020; Ho et al., 2020; Song et al., 2020a; Ramesh et al., 2022; Rombach et al., 2022) , audio synthesis (Chen et al., 2020; Kong et al., 2020) , and video generation (Ho et al., 2022; Blattmann et al., 2023) . In addition to the successful empirical results, the convergence analysis of SGM has attracted significant attention in the recent literature (Lee et al., 2022 (Lee et al., , 2023;; Chen et al., 2023; Li et al., 2023 Li et al., , 2024)) .Stochastic gradient Langevin dynamics (SGLD) (Welling and Teh, 2011) , as a fundamental methodology to implement and interpret SGM, can produce samples from the (Stein) score function of a probability density, i.e., the gradient of the log probability density function with respect to data. It has been widely recognized that a pitfall of SGLD is its slow mixing rate (Wooddard et al., 2009; Raginsky et al., 2017; Lee et al., 2018) . Specifically, Song and Ermon (2019) shows that under a multi-modal data distribution, the samples from Langevin dynamics may have an incorrect relative density across the modes. Based on this finding, Song and Ermon (2019) proposes anneal Langevin dynamics, which injects different levels of Gaussian noise into the data distribution and samples with SGLD on the perturbed distribution. While outputting the correct relative density across modes can be challenging for SGLD, a natural question is whether SGLD would be able to find all the modes of a multi-modal distribution.Submitted to 38th Conference on Neural Information Processing Systems (NeurIPS 2024). Do not distribute.In this work, we study this question by analyzing the mode-seeking properties of SGLD. The notion of mode-seekingness (Bishop, 2006; Ke et al., 2021; Li and Farnia, 2023) refers to the property that a generative model captures only a subset of the modes of a multi-modal distribution. We note that a similar problem, known as metastability, has been studied in the context of Langevin diffusion, a continuous-time version of SGLD described by stochastic differential equation (SDE) (Bovier et al., 2002 (Bovier et al., , 2004;; Gayrard et al., 2005) . Specifically, Bovier et al. (2002) gave a sharp bound on the mean hitting time of Langevin diffusion and proved that it may require exponential (in the space dimensionality d) time for transition between modes. Regarding discrete SGLD, Lee et al. (2018) constructed a probability distribution whose density is close to a mixture of two well-separated isotropic Gaussians, and proved that SGLD could not find one of the two modes within an exponential number of steps. However, further exploration of mode-seeking tendencies of SGLD and its variants such as annealed Langevin dynamics for general distributions is still lacking in the literature.In this work, we theoretically formulate and demonstrate the potential mode-seeking tendency of SGLD. We begin by analyzing the convergence under a variety of Gaussian mixture probability distributions, under which SGLD could fail to visit all the mixture components within sub-exponential steps (in the data dimension). Subsequently, we generalize this result to mixture distributions with sub-Gaussian modes. This generalization extends our earlier result on Gaussian mixtures to a significantly larger family of mixture models, as the sub-Gaussian family includes any distribution over an ℓ 2 -norm-bounded support set. Furthermore, we extend our theoretical results to anneal Langevin dynamics with bounded noise scales.To reduce SGLD's large iteration complexity shown under a high-dimensional input vector, we propose Chained Langevin Dynamics (Chained-LD). Since SGLD could suffer from the curse of dimensionality, we decompose the sample x ∈ R d into d/Q patches x (1) , • • • , x (d/Q) , each of constant size Q, and sequentially generate every patch x (q) for all q ∈ [d/Q] statistically conditioned on previous patches, i.e., P (x (q) | x (0) , • • • x (q-1) ). The combination of all patches generated from the conditional distribution faithfully follows the probability density P (x), while learning each patch requires less cost due to the reduced dimension. We also provide a theoretical analysis of Chained-LD by reducing the convergence of a d-dimensional sample to the convergence of each patch.Finally, we present the results of several numerical experiments to validate our theoretical findings.For synthetic experiments, we consider moderately high-dimensional Gaussian mixture models, where the vanilla and annealed Langevin dynamics could not find all the components within a million steps, while Chained-LD could capture all the components with correct frequencies in O(10 4 ) steps.For experiments on real image datasets, we consider a mixture of two modes by using the original images from MNIST/Fashion-MNIST training dataset (black background and white digits/objects) as the first mode and constructing the second mode by i.i.d. flipping the images (white background and black digits/objects) with probability 0.5. Following from Song and Ermon (2019) , we trained a Noise Conditional Score Network (NCSN) to estimate the score function. Our numerical results indicate that vanilla Langevin dynamics can fail to capture the two modes, as also observed by Song and Ermon (2019) . On the other hand, Chained-LD was capable of finding both modes regardless of initialization. We summarize the contributions of this work as follows:• Theoretically studying the mode-seeking properties of vanilla and annealed Langevin dynamics,• Proposing Chained Langevin Dynamics (Chained-LD), which decomposes the sample into patches and sequentially generates each patch conditioned on previous patches,• Providing a theoretical analysis of the convergence behavior of Chained-LD,• Numerically comparing the mode-seeking properties of vanilla, annealed, and chained Langevin dynamics.Notations: We use [n] to denote the set {1, 2, • • • , n}. Also, in the paper, ∥•∥ refers to the ℓ 2 norm.We use 0 n and 1 n to denote a 0-vector and 1-vector of length n. We use I n to denote the identity matrix of size n × n. In the text, TV stands for the total variation distance. Roberts and Tweedie, 1996; Bakry and Émery, 1983; Bakry et al., 2008) . Langevin dynamics, also known as Langevin Monte Carlo, is a discretization of Langevin diffusion typically modeled as a Markov Chain Monte Carlo (Welling and Teh, 2011) . For unimodal distributions, e.g., the probability density function that is log-concave or satisfies log-Sobolev inequality, the convergence of Langevin dynamics is provably fast (Dalalyan, 2017; Durmus and Moulines, 2017; Vempala and Wibisono, 2019) . However, for multimodal distributions, the non-asymptotic convergence analysis is much more challenging (Cheng et al., 2018) . Raginsky et al. (2017) gave an upper bound on the convergence time of Langevin dynamics for arbitrary non-log-concave distributions with certain regularity assumptions, which, however, could be exponentially large without imposing more restrictive assumptions. Lee et al. (2018) studied the special case of a mixture of Gaussians of equal variance and provided heuristic analysis of sampling from general non-log-concave distributions.The investigation of the mode-seekingness of generative models starts with different generative adversarial network (GAN) (Goodfellow et al., 2014) model formulations and divergence measures, from both the practical (Goodfellow, 2016; Poole et al., 2016) and theoretical (Shannon et al., 2020; Li and Farnia, 2023) perspectives. In the context of Langevin dynamics, mode-seekingness is closely related to a lower bound on the transition time between two modes, e.g., two local maximums. Bovier et al. (2002 Bovier et al. ( , 2004)) ; Gayrard et al. (2005) studied the mean hitting time of the continuous Langevin diffusion. Lee et al. (2018) proved the existence of a mixture of two Gaussian distributions whose covariance matrices differ by a constant factor, Langevin dynamics cannot find both modes in polynomial time.Score-based Generative Modeling: Since Song et al. (2020b) proposed sliced score matching which can train deep models to learn the score functions of implicit probability distributions on highdimensional data, score-based generative modeling (SGM) has been going through a spurt of growth.Annealed Langevin dynamics (Song and Ermon, 2019) estimates the noise score of the probability density perturbed by Gaussian noise and utilizes stochastic gradient Langevin dynamics to generate samples from a sequence of decreasing noise scales. Song and Ermon (2020) conducted a heuristic analysis of the effect of noise levels on the performance of annealed Langevin dynamics. Denoising diffusion probabilistic model (DDPM) (Ho et al., 2020) incorporates a step-by-step introduction of random noise into data, followed by learning to reverse this diffusion process in order to generate desired data samples from the noise. Song et al. (2020c) unified anneal Langevin dynamics and DDPM via a stochastic differential equation. A recent line of work focuses on the non-asymptotic convergence guarantees for SGM with an imperfect score estimation under various assumptions on the data distribution (Block et al., 2020; De Bortoli et al., 2021; Lee et al., 2022; Chen et al., 2023; Benton et al., 2023; Li et al., 2023 Li et al., , 2024)) . 
limitations: 
 Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] Justification: The limitations of this work are discussed in Section 7.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.• The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.• The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.• The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.• If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
langevin dynamics: 
 Generative modeling aims to produce samples such that their distribution is close to the underlying true distribution P . For a continuously differentiable probability density P (x) on R d , its score function is defined as the gradient of the log probability density function (PDF) ∇ x log P (x). Langevin diffusion is a stochastic process defined by the stochastic differential equation (SDE)dx t = -∇ x log P (x t ) dt + √ 2 dw t ,where w t is the Wiener process on R d . To generate samples from Langevin diffusion, Welling and Teh (2011) proposed stochastic gradient Langevin dynamics (SGLD), a discretization of the SDE for T iterations. Each iteration of SGLD is defined asEQUATIONwhere δ t is the step size and ϵ t ∼ N (0 d , I d ) is Gaussian noise. It has been widely recognized that Langevin diffusion could take exponential time to mix without additional assumptions on the probability density (Bovier et al., 2002 (Bovier et al., , 2004;; Gayrard et al., 2005; Raginsky et al., 2017; Lee et al., 2018) . To combat the slow mixing, Song and Ermon (2019) proposed annealed Langevin dynamics by perturbing the probability density with Gaussian noise of variance σ 2 , i.e.,EQUATIONand running SGLD on the perturbed data distribution P σt (x) with gradually decreasing noise levelsEQUATIONwhere δ t is the step size and ϵ t ∼ N (0 d , I d ) is Gaussian noise. When the noise level σ is vanishingly small, the perturbed distribution is close to the true distribution, i.e., P σ (x) ≈ P (x). Since we do not have direct access to the (perturbed) score function, Song and Ermon (2019) proposed the Noise Conditional Score Network (NCSN) s θ (x, σ) to jointly estimate the scores of all perturbed data distributions, i.e.,∀σ ∈ {σ t } t∈[T ] , s θ (x, σ) ≈ ∇ x log P σ (x).To train the NCSN, Song and Ermon (2019) adopted denoising score matching, which minimizes the following lossL θ; {σ t } t∈[T ] := 1 2T t∈[T ] σ 2 t E x∼P E x∼N (x,σ 2 t I d ) s θ (x, σ t ) - x -x σ 2 t 2 .Assuming the NCSN has enough capacity, s θ * (x, σ) minimizes the loss L θ;{σ t } t∈[T ] if and only if s θ * (x, σ t ) = ∇ x log P σt (x) almost surely for all t ∈ [T ].Our work focuses on multi-modal distributions. We use P = i∈[k] w i P (i) to represent a mixture of k modes, where each mode P (i) is a probability density with frequency w i such that w i > 0for all i ∈ [k] and i∈[k] w i = 1.In our theoretical analysis, we consider Gaussian mixtures and sub-Gaussian mixtures, i.e., every component P (i) is a Gaussian or sub-Gaussian distribution. A probability distribution p(z) of dimension d is defined as a sub-Gaussian distribution with parameter ν 2 if, given the mean vector µ := E z∼p [z], the moment generating function (MGF) of p satisfies the following inequality for every vector α ∈ R d :EQUATION)We remark that sub-Gaussian distributions include a wide variety of distributions such as Gaussian distributions and any distribution within a bounded ℓ 2 -norm distance from the mean µ. From equation 2 we note that the perturbed distribution is the convolution of the original distribution and a Gaussian random variable, i.e., for random variables z ∼ p and t ∼ N (0 d , I d ), their sum z+t ∼ p σ follows the perturbed distribution with noise level σ. Therefore, a perturbed (sub)Gaussian distribution remains (sub)Gaussian. We formalize this property in Proposition 1 and defer the proof to Appendix A for completeness.Proposition 1. Suppose the perturbed distribution of a d-dimensional probability distribution p with noise level σ is p σ , then the mean of the perturbed distribution is the same as the original distribution,i.e., E z∼pσ [z] = E z∼p [z]. If p = N (µ, Σ) is a Gaussian distribution, p σ = N (µ, Σ + σ 2 I d ) is also a Gaussian distribution. If p is a sub-Gaussian distribution with parameter ν 2 , p σ is a sub-Gaussian distribution with parameter (ν 2 + σ 2 ).4 Theoretical Analysis of the Mode-Seeking Properties of Langevin DynamicsIn this section, we theoretically investigate the mode-seeking properties of vanilla and annealed Langevin dynamics. We begin with analyzing Langevin dynamics in Gaussian mixtures.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
langevin dynamics in gaussian mixtures: 
 Assumption 1. Consider a data distribution P := k i=0 w i P (i) as a mixture of Gaussian distributions, where 1 ≤ k = o(d) and w i > 0 is a positive constant such that k i=0 w i = 1. Suppose thatP (i) = N (µ i , ν 2 i I d ) is a Gaussian distribution over R d for all i ∈ {0} ∪ [k] such that for all i ∈ [k], ν i < ν 0 and ∥µ i -µ 0 ∥ 2 ≤ ν 2 0 -ν 2 i 2 log ν 2 i ν 2 0 - ν 2 i 2ν 2 0 + ν 2 0 2ν 2 i d. Denote ν max := max i∈[k] ν i .Regarding the first requirement ν i < ν 0 , we first note that the probability density p(z) of a Gaussian distribution N (µ, ν 2 I d ) decays exponentially in terms of ∥z-µ∥ 2 ν 2. When a state z is sufficiently far from all modes (i.e., ∥z∥ ≫ ∥µ i ∥), the Gaussian distribution with the largest variance (i.e., P (0) in Assumption 1) dominates all other modes because ∥z-µ0∥ 2ν 2 0 ≈ ∥z∥ 2 ν 2 0 ≫ ∥z∥ 2 ν 2 i ≈ ∥z-µi∥ 2 ν 2 i. We call such mode P (0) the universal mode. Therefore, if z is initialized far from all modes, it can only converge to the universal mode because the gradient information of other modes is masked. Once z enters the universal mode P (0) , if the step size δ t of Langevin dynamics is small (i.e., δ t ≤ ν 2 0 ), it would take exponential steps to escape the local mode P (0) ; while if the step size is large (i.e., δ t > ν 2 0 ), the state z would again be far from all modes and thus the universal mode P (0) dominates all other modes. Hence, z can only visit the universal mode unless the stochastic noise ϵ t miraculously leads it to the region of another mode. In addition, it can be verified that logν 2 i ν 2 0 - ν 2 i 2ν 2 0 + ν 2 0 2ν 2 i is apositive constant for ν i < ν 0 , thus the second requirement of Assumption 1 essentially represents∥µ i -µ 0 ∥ 2 ≤ O(d).We formalize the intuition in Theorem 1 and defer the proof to Appendix A.1.Theorem 1. Consider a data distribution P satisfying Assumption 1. We follow Langevin dynamics for T = exp(O(d)) steps. Suppose the sample is initialized in P (0) , then with probability at least1 -T • exp(-Ω(d)), we have ∥x t -µ i ∥ 2 > ν 2 0 +ν 2 max 2 d for all t ∈ {0} ∪ [T ] and i ∈ [k].We note that ∥x t -µ i ∥ 2 > ν 2 0 +ν 2 max 2d is a strong notion of mode-seekingness, since the probability density of modeP (i) = N (µ i , ν 2 i I d ) concentrates around the ℓ 2 -norm ball z : ∥z -µ i ∥ 2 ≤ ν 2 i d .This notion can also easily be translated into a lower bound in terms of other distance measures such as total variation distance and Wasserstein 2-distance. Moreover, in Theorem 2 we extend the result to annealed Langevin dynamics with bounded noise level, and the proof is deferred to Appendix A.2.Theorem 2. Consider a data distribution P satisfying Assumption 1. We follow annealed Langevin dynamics for T = exp(O(d)) steps with noise levelsc σ ≥ σ 0 ≥ • • • ≥ σ T ≥ 0 for constant c σ > 0.In addition, assume for all i ∈[k], ∥µ i -µ 0 ∥ 2 ≤ ν 2 0 -ν 2 i 2 log ν 2 i +c 2 σ ν 2 0 +c 2 σ - ν 2 i +c 2 σ 2ν 2 0 +c 2 σ + ν 2 0 +c 2 σ 2ν 2 i +c 2 σ d.Suppose that the sample is initialized in P (0) σ0 , then with probability at least 1 -T • exp(-Ω(d)), wehave ∥x t -µ i ∥ 2 > ν 2 0 +ν 2 max +2σ 2 t 2 d for all t ∈ {0} ∪ [T ] and i ∈ [k].We further generalize our results to sub-Gaussian mixtures. We impose the following assumptions on the mixture. It is worth noting that these assumptions automatically hold for Gaussian mixtures.Assumption 2. Consider a data distribution P := k i=0 w i P (i) as a mixture of sub-Gaussian distributions, where1 ≤ k = o(d) and w i > 0 is a positive constant such that k i=0 w i = 1.Suppose thatP (0) = N (µ 0 , ν 2 0 I d ) is Gaussian and for all i ∈ [k], P (i) satisfies i. P (i) is a sub-Gaussian distribution of mean µ i with parameter ν 2 i ,ii. P (i) is differentiable and∇P (i) (µ i ) = 0 d ,iii. the score function ofP (i) is L i -Lipschitz such that L i ≤ c L ν 2 i for some constant c L > 0, iv. ν 2 0 > max 1, 4(c 2 L +cν c L ) cν (1-cν ) ν 2 max 1-cν for constant c ν ∈ (0, 1), where ν max := max i∈[k] ν i , v. ∥µ i -µ 0 ∥ 2 ≤ (1-cν )ν 2 0 -ν 2 i 2(1-cν ) log cν ν 2 i (c 2 L +cν c L )ν 2 0 - ν 2 i 2(1-cν )ν 2 0 + (1-cν )ν 2 0 2ν 2 i d.We EQUATIONt , • • • , x (q-1) t + √ δ t ϵ t , where ϵ t ∼ N (0 Q , I Q ) 5:end for 6:x(q) 0 ← x (q) T Q/d 7: end for 8: return x T Q/dTheorem 3. Consider a data distribution P satisfying Assumption 2. We follow Langevin dynamics for T = exp(O(d)) steps. Suppose the sample is initialized in P (0) , then with probability at least1 -T • exp(-O(d)), we have ∥x t -µ i ∥ 2 > ν 2 0 2 + ν 2 max 2(1-cν ) d for all t ∈ {0} ∪ [T ] and i ∈ [k].Finally, we slightly modify Assumption 2 and extend our results to annealed Langevin dynamics under sub-Gaussian mixtures in Theorem 4. The details of Assumption 3 and the proof of Theorem 4 are deferred to Appendix A.4.Theorem 4. Consider a data distribution P satisfying Assumption 3. We follow annealed Langevindynamics for T = exp(O(d)) steps with noise levels c σ ≥ σ 0 ≥ • • • ≥ σ T ≥ 0. Suppose the sample is initialized in P (0)σ0 , then with probability at least 1 -T • exp(-O(d)), we have∥x t -µ i ∥ 2 > ν 2 0 +σ 2 t 2 + ν 2 max +σ 2 t 2(1-cν ) d for all t ∈ {0} ∪ [T ] and i ∈ [k]. 
chained langevin dynamics: 
 To reduce the mode-seeking tendencies of vanilla and annealed Langevin dynamics, we propose Chained Langevin Dynamics (Chained-LD) in Algorithm 1. While vanilla and annealed Langevin dynamics apply gradient updates to all coordinates of the sample in every step, we decompose the sample into patches of constant size and generate each patch sequentially to alleviate the exponential dependency on the dimensionality. More precisely, we divide a sample x into d/Q patchesx (1) , • • • x (d/Q) of some constant size Q, and apply annealed Langevin dynamics to sample each patch x (q) (for q ∈ [d/Q]) from the conditional distribution P (x (q) | x (1) , • • • x (q-1) ).An ideal conditional score function estimator s θ could jointly estimate the scores of all perturbed conditional patch distribution, i.e., ∀σ∈ {σ t } t∈[T Q/d] , q ∈ [d/Q], s θ x (q) | σ, x (1) , • • • , x (q-1) ≈ ∇ x (q) log P σ (x (q) | x (1) , • • • x (q-1) ).Following from Song and Ermon (2019) , we use the denoising score matching to train the estimator.For a given σ, the denoising score matching objective isℓ(θ; σ) := 1 2 E x∼P E x∼N (x,σ 2 I d ) q∈[d/Q] s θ x (q) | σ, x (1) , • • • , x (q-1) - x(q) -x (q) σ 2 2 .Then, combining the objectives gives the following lossL θ; {σ t } t∈[T Q/d] := d T Q t∈[T Q/d] σ 2 t ℓ(θ; σ t ).As shown in Vincent ( 2011), an estimator s θ with enough capacity minimizes the loss L if and only if s θ outputs the scores of all perturbed conditional patch distribution almost surely. Ideally, if a sampler perfectly generates every patch, combining all patches gives a sample from the original distributionsince P (x) = q∈[d/Q] P (x (q) | x (1) , • • • x (q-1)). In Theorem 5 we give a linear reduction from producing samples of dimension d using Chained-LD to learning the distribution of a Q-dimensional variable for constant Q. The proof of Theorem 5 is deferred to Appendix A.5.Theorem 5. Consider a sampler algorithm taking the first q -1 patches x (1) , • • • , x (q-1) as input and outputing a sample of the next patch x (q) with probability Px (q) | x (1) , • • • , x (q-1) for all q ∈ [d/Q].Suppose that for every q ∈ [d/Q] and any given previous patches x (1) , • • • , x (q-1) , the sampler algorithm can achieveTV P x (q) | x (1) , • • • , x (q-1) , P x (q) | x (1) , • • • , x (q-1) ≤ ε • Q d in τ (ε, d) iterations for some ε > 0. Then, equipped with the sampler algorithm, the Chained-LD algorithm in d Q • τ (ε, d) iterations can achieve TV P (x), P (x) ≤ ε.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes] Justification: We submitted our code in supplementary materials for the readers to reproduce our numerical results.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details.• While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).• The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https://nips.cc/public/guides/CodeSubmissionPolicy) for more details.• The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.• The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.• At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).• Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
numerical results: 
 In this section, we empirically evaluated the mode-seeking tendencies of vanilla, annealed, and chained Langevin dynamics. We performed numerical experiments on synthetic Gaussian mixture models and real image datasets including MNIST (LeCun, 1998) and Fashion-MNIST (Xiao et al., 2017) . Details on the experiment setup are deferred to Appendix B. andP (2) = N (-1 d , I d ).The frequencies of the three modes are 0.2, 0.4 and 0.4, i.e.,P = 0.2P (0) + 0.4P (1) + 0.4P (2) = 0.2N (0 d , 3I d ) + 0.4N (1 d , I d ) + 0.4N (-1 d , I d ).As shown in Figure 1 , vanilla and annealed Langevin dynamics cannot find mode 1 or 2 within 10 6 iterations if the sample is initialized in mode 0, while chained Langevin dynamics can find the other two modes in 1000 steps and correctly recover their frequencies as gradually increasing the number of iterations. In Appendix B.1 we present additional experiments on samples initialized in mode 1 or 2, which also verify the mode-seeking tendencies of vanilla and annealed Langevin dynamics.Image datasets: We construct the distribution as a mixture of two modes by using the original images from MNIST/Fashion-MNIST training dataset (black background and white digits/objects) as the first mode and constructing the second mode by i.i.d. randomly flipping an image (white background and black digits/objects) with probability 0.5. Regarding the neural network architecture of the score function estimator, for vanilla and annealed Langevin dynamics we use U-Net (Ronneberger et al., 2015) following from Song and Ermon (2019) . For chained Langevin dynamics, we proposed to use Recurrent Neural Network (RNN) architectures. We note that for a sequence of inputs, the output of RNN from the previous step is fed as input to the current step. Therefore, in the scenario of chained Langevin dynamics, the hidden state of RNN contains information about the previous patches and allows the network to estimate the conditional score function ∇ x (q) log P (x (q) | x (1) , • • • x (q-1) ).More implementation details are deferred to Appendix B.2. Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: All details of the numerical experiments are listed in Appendix B.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.• The full details can be provided either with the code, in appendix, or as supplemental material. 
conclusion: 
 In this work, we theoretically and numerically studied the mode-seeking properties of vanilla and annealed Langevin dynamics sampling methods under a multi-modal distribution. We characterized Gaussian and sub-Gaussian mixture models under which Langevin dynamics are unlikely to find all the components within a sub-exponential number of iterations. To reduce the mode-seeking tendency of vanilla Langevin dynamics, we proposed Chained Langevin Dynamics (Chained-LD) and analyzed its convergence behavior. Studying the connections between Chained-LD and denoising diffusion models will be an interesting topic for future exploration.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: For the results of synthetic data, we report 1000 samples for every experiment.For the results of image datasets, we report 49 or 100 samples for every experiment.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).• If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: Information on the computational resources of the numerical experiments is provided in Appendix B.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.• The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: We read the guidelines and wrote the paper according to the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: As discussed in Appendix C, this work has potential societal consequences, none of which needs to be specially highlighted.Guidelines: 
 
 Appendix: Our RNN-based implementation of Chained-LD is currently limited to image data generation tasks.An interesting future direction is to extend the application of Chained-LD to other domains such as audio and text data. Another future direction could be to study the convergence of Chained-LD under an imperfect score estimation which we did not address in our analysis.We begin by introducing some well-established lemmas used in our proof. We first provide the proof of Proposition 1 for completeness:Proof of Proposition 1. By the definition in equation 2, we havep σ (z) = p(t)N (z | t, σ 2 I d ) dt = p(t)N (z -t | 0 d , σ 2 I d ) dt.For random variables t ∼ p and y ∼ N (0 d , I d ), their sum z = t + y ∼ p σ follows the perturbed distribution with noise level σ. Therefore,E z∼pσ [z] = E (t+y)∼pσ [t + y] = E t∼p [t] + E y∼N (0 d ,I d ) [y] = E t∼p [t]. If t ∼ p = N (µ, Σ) follows a Gaussian distribution, we have z = t + y ∼ p σ = N (µ, Σ + σ 2 I d ).If p is a sub-Gaussian distribution with parameter ν 2 , we have z = t + y ∼ p σ is a sub-Gaussian distribution with parameter (ν 2 + σ 2 ). Hence we obtain Proposition 1.We use the following lemma on the tail bound for multivariate Gaussian random variables.Lemma 1 (Lemma 1, Laurent and Massart (2000) ). Suppose that a random variable z ∼ N (0 d , I d ).Then for any λ > 0,P ∥z∥ 2 ≥ d + 2 √ dλ + 2λ ≤ exp(-λ), P ∥z∥ 2 ≤ d -2 √ dλ ≤ exp(-λ).We also use a tail bound for one-dimensional Gaussian random variables and provide the proof here for completeness.Lemma 2. Suppose a random variable Z ∼ N (0, 1). Then for any t > 0,P(Z ≥ t) = P(Z ≤ -t) ≤ exp(-t 2 /2) √ 2πt .Proof of Lemma 2. Since z t ≥ 1 for all z ∈ [t, ∞), we haveP(Z ≥ t) = 1 √ 2π ∞ t exp - z 2 2 dz ≤ 1 √ 2π ∞ t z t exp - z 2 2 dz = exp(-t 2 /2) √ 2πt .Since the Gaussian distribution is symmetric, we have P(Z ≥ t) = P(Z ≤ -t). Hence we obtain the desired bound.Without loss of generality, we assume that µ 0 = 0 d for simplicity. Let r and n respectively denote the rank and nullity of the vector space{µ i } i∈[k] , then we have r + n = d and 0 ≤ r ≤ k = o(d).Denote R ∈ R d×r an orthonormal basis of the vector space {µ i } i∈[k] , and denote N ∈ R d×n an orthonormal basis of the null space of {µ i } i∈ [k] . Now consider decomposing the sample x t by r t := R T x t , andn t := N T x t ,where r t ∈ R r , n t ∈ R n . Then we havex t = Rr t + Nn t .Similarly, we decompose the noise ϵ t into ϵ (r)t := R T ϵ t , and ϵ (n) t := N T ϵ t ,where ϵ (r)t ∈ R r , ϵ (n) t ∈ R n .Then we haveϵ t = Rϵ (r) t + Nϵ (n) t .Since a linear combination of a Gaussian random variable still follows Gaussian distribution, byϵ t ∼ N (0 d , I d ), R T R = I r , and N T N = I n we obtain ϵ (r)t ∼ N (0 r , I r ), and ϵ(n) t ∼ N (0 n , I n ).By the definition of Langevin dynamics in equation 1, the two components of x t follow from the update rule:EQUATIONr t = r t-1 + δ t 2 R T ∇ x log P (x t-1 ) + δ t ϵ (r) t .It is worth noting that sinceN T µ i = 0 n . To show ∥x t -µ i ∥ 2 > ν 2 0 +ν 2 max 2 d, it suffices to prove ∥n t ∥ 2 > ν 2 0 + ν 2 max 2 d.We start by proving that the initialization of the state x 0 has a large norm on the null space with high probability in the following proposition.Proposition 2. Suppose that a sample x 0 is initialized in the distribution P (0) , i.e., x 0 ∼ P (0) , then for any constant ν max < ν 0 , with probability at least 1 -exp(-Ω(d)), we have∥n 0 ∥ 2 ≥ 3ν 2 0 +ν 2 max 4 d. Proof of Proposition 2. Since x 0 ∼ P (0) = N (0 d , ν 2 0 I d ) and N T N = I n , we know n 0 = N T x 0 ∼ N (0 n , ν 2 0 I n ).Therefore, by Lemma 1 we can boundP ∥n 0 ∥ 2 ≤ 3ν 2 0 + ν 2 max 4 d = P   ∥n 0 ∥ 2 ν 2 0 ≤ d -2 d • ν 2 0 -ν 2 max 8ν 2 0 2 d   ≤ P   ∥n 0 ∥ 2 ν 2 0 ≤ n -2 n ν 2 0 -ν 2 max 8ν 2 0 2 d 2   ≤ exp - ν 2 0 -ν 2 max 8ν 2 0 2 d 2 ,where the second last step follows from the assumption d -n = r = o(d). Hence we complete the proof of Proposition 2.Then, with the assumption that the initialization satisfies ∥n 0 ∥ 2 ≥ 3ν 2 0 +ν 2 max 4 d, the following proposition shows that ∥n t ∥ remains large with high probability.Proposition 3. Consider a data distribution P satisfies the constraints specified in Theorem 1.We follow the Langevin dynamics for T = exp(O(d)) steps. Suppose that the initial samplesatisfies ∥n 0 ∥ 2 ≥ 3ν 2 0 +ν 2 max 4d, then with probability at least 1 -T • exp(-Ω(d)), we have that∥n t ∥ 2 > ν 2 0 +ν 2 max 2 d for all t ∈ {0} ∪ [T ].Proof of Proposition 3. To establish a lower bound on ∥n t ∥, we consider different cases of the step size δ t . Intuitively, when δ t is large enough, n t will be too noisy due to the introduction of random noise √ δ t ϵ(n) t in equation 5. While for small δ t , the update of n t is bounded and thus we can iteratively analyze n t . We first handle the case of large δ t in the following lemma.Lemma 3. If δ t > ν 2 0 , with probability at least 1 -exp(-Ω(d)), for n t satisfying equation 5, wehave ∥n t ∥ 2 ≥ 3ν 2 0 +ν 2 max 4d regardless of the previous state x t-1 .Proof of Lemma 3. Denote v := n t-1 + δt 2 N T ∇ x log P (x t-1 ) for simplicity. Note that v is fixed for any given x t-1 . We decompose ϵ(n) t into a vector aligning with v and another vector orthogonalto v. Consider an orthonormal matrix M ∈ R n×(n-1) such that M T v = 0 n-1 and M T M = I n-1 . By denoting u := ϵ (n) t -MM T ϵ (n) twe have M T u = 0 n-1 , thus we obtain∥n t ∥ 2 = v + δ t ϵ (n) t 2 = v + δ t u + δ t MM T ϵ (n) t 2 = v + δ t u 2 + δ t MM T ϵ (n) t 2 ≥ δ t MM T ϵ (n) t 2 ≥ ν 2 0 M T ϵ (n) t 2 . Since ϵ (n) t ∼ N (0 n , I n ) and M T M = I n-1 , we obtain M T ϵ (n) t ∼ N (0 n-1 , I n-1). Therefore, by Lemma 1 we can boundP ∥n t ∥ 2 ≤ 3ν 2 0 + ν 2 max 4 d ≤ P M T ϵ (n) t 2 ≤ 3ν 2 0 + ν 2 max 4ν 2 0 d = P   M T ϵ (n) t 2 ≤ d -2 d • ν 2 0 -ν 2 max 8ν 2 0 2 d   ≤ P   M T ϵ (n) t 2 ≤ (n -1) -2 (n -1) ν 2 0 -ν 2 max 8ν 2 0 2 d 2   ≤ exp - ν 2 0 -ν 2 max 8ν 2 0 2 d 2 ,where the second last step follows from the assumption d -n = r = o(d). Hence we complete the proof of Lemma 3.We then consider the case when δ t ≤ ν 2 0 . Let r := R T x and n := N T x, then x = Rr + Nn. We first show that when ∥n∥ 2 ≥ν 2 0 +ν 2 max 2 d, P (i) (x) is exponentially smaller than P (0) (x) for all i ∈ [k]in the following lemma.Lemma 4. Given that ∥n∥ 2 ≥ ν 2 0 +ν 2 max 2 d and ∥µ i ∥ 2 ≤ ν 2 0 -ν 2 i 2 log ν 2 i ν 2 0 - ν 2 i 2ν 2 0 + ν 2 0 2ν 2 i d for all i ∈ [k], we have P (i) (x) P (0) (x) ≤ exp(-Ω(d)) for all i ∈ [k]. Proof of Lemma 4. For all i ∈ [k], define ρ i (x) := P (i) (x) P (0) (x) , then ρ i (x) = P (i) (x) P (0) (x) = (2πν 2 i ) -d/2 exp -1 2ν 2 i ∥x -µ i ∥ 2 (2πν 2 0 ) -d/2 exp -1 2ν 2 0 ∥x∥ 2 = ν 2 0 ν 2 i d/2 exp 1 2ν 2 0 ∥x∥ 2 - 1 2ν 2 i ∥x -µ i ∥ 2 = ν 2 0 ν 2 i d/2 exp 1 2ν 2 0 - 1 2ν 2 i ∥Nn∥ 2 + ∥Rr∥ 2 2ν 2 0 - ∥Rr -µ i ∥ 2 2ν 2 i = ν 2 0 ν 2 i d/2 exp 1 2ν 2 0 - 1 2ν 2 i ∥n∥ 2 + ∥r∥ 2 2ν 2 0 - r -R T µ i 2 2ν 2 i ,where the last step follows from the definition that R ∈ R d×r an orthonormal basis of the vector space{µ i } i∈[k] and N T N = I n . Since ν 2 0 > ν 2 i , the quadratic term ∥r∥ 2 2ν 2 0 - ∥r-R T µi∥ 2 2ν 2 i is maximized at r = ν 2 0 R T µi ν 2 0 -ν 2 i . Therefore, ∥r∥ 2 2ν 2 0 - r -R T µ i 2 2ν 2 i ≤ ν 4 0 R T µ i 2 2ν 2 0 (ν 2 0 -ν 2 i ) 2 - 1 2ν 2 i ν 2 0 ν 2 0 -ν 2 i -1 2 R T µ i 2 = ∥µ i ∥ 2 2(ν 2 0 -ν 2 i ).Hence, for ∥n∥ 2 ≥ ν 2 0 +ν 2 max 2 d and ∥µ i ∥ 2 ≤ ν 2 0 -ν 2 i 2 log ν 2 i ν 2 0 - ν 2 i 2ν 2 0 + ν 2 0 2ν 2 i d, we have ρ i (x) = ν 2 0 ν 2 i d/2 exp 1 2ν 2 0 - 1 2ν 2 i ∥n∥ 2 + ∥r∥ 2 2ν 2 0 - r -R T µ i 2 2ν 2 i ≤ ν 2 0 ν 2 i d/2 exp 1 2ν 2 0 - 1 2ν 2 i ν 2 0 + ν 2 i 2 d + ∥µ i ∥ 2 2(ν 2 0 -ν 2 i ) = exp -log ν 2 i ν 2 0 - ν 2 i 2ν 2 0 + ν 2 0 2ν 2 i d 2 + ∥µ i ∥ 2 2(ν 2 0 -ν 2 i ) ≤ exp -log ν 2 i ν 2 0 - ν 2 i 2ν 2 0 + ν 2 0 2ν 2 i d 4 .Notice that for functionf (z) = log z -z 2 + 1 2z , we have f (1) = 0 and d dz f (z) = 1 z -1 2 -1 2z 2 = -1 2 1 z -1 2 < 0 when z ∈ (0, 1). Thus, log ν 2 i ν 2 0 - ν 2 i 2ν 2 0 + ν 2 0 2ν 2 i is a positive constant for ν i < ν 0 , i.e., ρ i (x) = exp(-Ω(d)). Therefore we finish the proof of Lemma 4.Lemma 4 implies that when ∥n∥ is large, the Gaussian mode P (0) dominates other modes P (i) . To bound ∥n t ∥, we first consider a simpler case that ∥n t-1 ∥ is large. Intuitively, the following lemma proves that when the previous state n t-1 is far from a mode, a single step of Langevin dynamics with bounded step size is not enough to find the mode.Lemma 5. Suppose δ t ≤ ν 2 0 and ∥n t-1 ∥ 2 > 36ν 2 0 d, then for n t following from equation 5, we have ∥n t ∥ 2 ≥ ν 2 0 d with probability at least 1 -exp(-Ω(d)).Proof of Lemma 5. From the recursion of n t in equation 5 we haveEQUATIONBy Lemma 4, we haveEQUATIONOn the other hand, from ϵ(n) t ∼ N (0 n , I n ) we know ⟨nt-1,ϵ (n) t ⟩ ∥nt-1∥∼ N (0, 1) for any fixed n t-1 ̸ = 0 n , hence by Lemma 2 we haveEQUATIONCombining equation 6, equation 7 and equation 8 gives that∥n t ∥ 2 ≥ 1 3 2 ∥n t-1 ∥ 2 -2ν 0 |⟨n t-1 , ϵ (n) t ⟩| ≥ 1 9 ∥n t-1 ∥ 2 - ν 0 √ d 2 ∥n t-1 ∥ ≥ 1 9 • 36ν 2 0 d - ν 0 √ d 2 • 6ν 0 √ d = ν 2 0 d with probability at least 1 -8 √ 2πd exp -d 32 = 1 -exp(-Ω(d)). This proves Lemma 5.We then proceed to bound ∥n t ∥ iteratively for ∥n t-1 ∥ 2 ≤ 36ν 2 0 d. Recall that equation 5 givesn t = n t-1 + δ t 2 N T ∇ x log P (x t-1 ) + δ t ϵ (n) t .We notice that the difficulty of solving n t exhibits in the dependence of log P (x t-1 ) on r t-1 . SinceP = k i=0 w i P (i) = k i=0 w i N (µ i , ν 2 i I d ),we can rewrite the score function asEQUATIONNow, instead of directly working with n t , we consider a surrogate recursion nt such that n0 = n 0 and for all t ≥ 1,EQUATIONThe advantage of the surrogate recursion is that nt is independent of r, thus we can obtain the closed-form solution to nt . Before we proceed to bound nt , we first show that nt is sufficiently close to the original recursion n t in the following lemma.Lemma 6. For any t ≥ 1, given that δ j ≤ ν 2 0 andν 2 0 +ν 2 max 2 d ≤ ∥n j-1 ∥ 2 ≤ 36ν 2 0 d for all j ∈ [t]and∥µ i ∥ 2 ≤ ν 2 0 -ν 2 i 2 log ν 2 i ν 2 0 - ν 2 i 2ν 2 0 + ν 2 0 2ν 2 i d for all i ∈ [k], we have ∥n t -n t ∥ ≤ t exp(Ω(d)) √ d.Proof of Lemma 6. Upon comparing equation 5 and equation 10, by equation 9 we have that for allEQUATIONEQUATIONP (i) (x j-1 ) P (0) (x j-1 ) 1 ν 2 i - 1 ν 2 0 6ν 0 √ d.By Lemma 4, we haveP (i) (xj-1) P (0) (xj-1) ≤ exp(-Ω(d)) for all i ∈ [k],hence we obtain a recursive bound∥n j -n j ∥ ≤ ∥n j-1 -n j-1 ∥ + 1 exp(Ω(d)) √ d.Finally, by n0 = n 0 , we have∥n t -n t ∥ = j∈[t] (∥n j -n j ∥ -∥n j-1 -n j-1 ∥) ≤ t exp(Ω(d)) √ d.Hence we obtain Lemma 6.We then proceed to analyze nt , The following lemma gives us the closed-form solution of nt . We slightly abuse the notations here, e.g.,c2 i=c1 1 -δi 2ν 2 0 = 1 and c2 j=c1 δ j = 0 for c 1 > c 2 . Lemma 7. For all t ≥ 0, nt ∼ N t i=1 1 -δi 2ν 2 0 n 0 , t j=1 t i=j+1 1 -δi 2ν 2 0 2 δ j I n ,where the mean and covariance satisfyt i=1 1 -δi 2ν 2 0 2 + 1 ν 2 0 t j=1 t i=j+1 1 -δi 2ν 2 0 2 δ j ≥ 1.Proof of Lemma 7. We prove the two properties by induction. When t = 0, they are trivial. Suppose they hold for t -1, then for the distribution of nt , we havent = nt-1 - δ t 2ν 2 0 nt-1 + δ t ϵ (n) t ∼ N   1 - δ t 2ν 2 0 t-1 i=1 1 - δ i 2ν 2 0 n 0 , 1 - δ t 2ν 2 0 2 t-1 j=1 t-1 i=j+1 1 - δ i 2ν 2 0 2 δ j I n + δ t I n   = N   t i=1 1 - δ i 2ν 2 0 n 0 , t j=1 t i=j+1 1 - δ i 2ν 2 0 2 δ j I n   .For the second property,t i=1 1 - δ i 2ν 2 0 2 + 1 ν 2 0 t j=1 t i=j+1 1 - δ i 2ν 2 0 2 δ j = 1 - δ t 2ν 2 0 2   t-1 i=1 1 - δ i 2ν 2 0 2 + 1 ν 2 0 t-1 j=1 t-1 i=j+1 1 - δ i 2ν 2 0 2 δ j   + 1 ν 2 0 δ t ≥ 1 - δ t 2ν 2 0 2 + 1 ν 2 0 δ t = 1 + δ 2 t 4ν 4 0 ≥ 1.Hence we finish the proof of Lemma 7.Armed with Lemma 7, we are now ready to establish the lower bound on ∥n t ∥. For simplicity, denote α := Proof of Lemma 8. By nt = αn 0 + √ βν 0 ϵ we have∥n t ∥ 2 = α 2 ∥n 0 ∥ 2 + βν 2 0 ∥ϵ∥ 2 + 2α βν 0 ⟨n 0 , ϵ⟩By Lemma 1 we can boundP ∥ϵ∥ 2 ≤ 3ν 2 0 + ν 2 max 4ν 2 0 d = P   ∥ϵ∥ 2 ≤ d -2 d • ν 2 0 -ν 2 max 8ν 2 0 2 d   ≤ P   ∥ϵ∥ 2 ≤ (n -1) -2 (n -1) ν 2 0 -ν 2 max 8ν 2 0 2 d 2   ≤ exp - ν 2 0 -ν 2 max 8ν 2 0 2 d 2 ,where the second last step follows from the assumption d-n = r = o(d). Since ϵ ∼ N (0 n , I n ),we know ⟨n0,ϵ⟩ ∥n0∥ ∼ N (0, 1). Therefore by Lemma 2,P ⟨n 0 , ϵ⟩ ∥n 0 ∥ ≤ - ν 2 0 -ν 2 max 4ν 0 3ν 2 0 + ν 2 max √ d ≤ 4ν 0 3ν 2 0 + ν 2 max √ 2π(ν 2 0 -ν 2 max ) √ d exp - (ν 2 0 -ν 2 max ) 2 d 32ν 2 0 (3ν 2 0 + ν 2 max ) Conditioned on ∥n 0 ∥ 2 ≥ 3ν 2 0 +ν 2 max 4 d, ∥ϵ∥ 2 > 3ν 2 0 +ν 2 max 4ν 2 0 d and 1 ∥n0∥ ⟨n 0 , ϵ⟩ > - ν 2 0 -ν 2 max 4ν0 √ 3ν 2 0 +ν 2 max √ d,since Lemma 7 gives α 2 + β ≥ 1 we have∥n t ∥ 2 = α 2 ∥n 0 ∥ 2 + βν 2 0 ∥ϵ∥ 2 + 2α βν 0 ⟨n 0 , ϵ⟩ ≥ α 2 ∥n 0 ∥ 2 + βν 2 0 ∥ϵ∥ 2 -2α βν 0 ∥n 0 ∥ ν 2 0 -ν 2 max 4ν 0 3ν 2 0 + ν 2 max √ d ≥ α 2 ∥n 0 ∥ 2 + βν 2 0 ∥ϵ∥ 2 -2α βν 0 ∥n 0 ∥ ∥ϵ∥ • ν 2 0 -ν 2 max 6ν 2 0 + 2ν 2 max ≥ 1 - ν 2 0 -ν 2 max 6ν 2 0 + 2ν 2 max α 2 ∥n 0 ∥ 2 + βν 2 0 ∥ϵ∥ 2 ≥ 5ν 2 0 + 3ν 2 max 6ν 2 0 + 2ν 2 max α 2 + β • 3ν 2 0 + ν 2 max 4 d ≥ 5ν 2 0 + 3ν 2 max 8 d.Hence by union bound, we complete the proof of Lemma 8.Upon having all the above lemmas, we are now ready to establish Proposition 3 by induction. Suppose the theorem holds for all T values of 1, • • • , T -1. We consider the following 3 cases:• If there exists some t ∈ [T ] such that δ t > ν 2 0 , by Lemma 3 we know that with probability ∥n T -n T ∥ < 5ν 2 0 + 3ν 2 max 8 - ν 2 0 + ν 2 max 2 √ d.By Lemma 8 we have that with probability at least 1 -exp(-Ω(d)),∥n T ∥ 2 ≥ 5ν 2 0 + 3ν 2 max 8 d.Combining the two inequalities implies the desired bound∥n T ∥ ≥ ∥n T ∥ -∥n T -n T ∥ > ν 2 0 + ν 2 max 2 d.Hence by induction we obtain ∥n t ∥ 2 > ν 2 0 +ν 2 max 2 d for all t ∈ [T ] with probability at least(1 -(T -1) exp(-Ω(d))) • (1 -exp(-Ω(d))) ≥ 1 -T exp(-Ω(d)).Therefore we complete the proof of Proposition 3.Finally, combining Propositions 2 and 3 finishes the proof of Theorem 1.To establish Theorem 2, we first note from Proposition 1 that perturbing a Gaussian distributionN (µ, ν 2 I d ) with noise level σ results in a Gaussian distribution N (µ, (ν 2 + σ 2 )I d ). Therefore, for a Gaussian mixture P = k i=0 w i P (i) = k i=0 w i N (µ i , ν 2 i I d ), the perturbed distribution of noise level σ is P σ = k i=0 w i N (µ i , (ν 2 i + σ 2 )I d ).Similar to the proof of Theorem 1, we decomposex t = Rr t + Nn t , andϵ t = Rϵ (r) t + Nϵ (n)t , where R ∈ R d×r an orthonormal basis of the vector space {µ i } i∈[k] and N ∈ R d×n an orthonormal basis of the null space of {µ i } i∈ [k] . Now, we prove Theorem 2 by applying the techniques developed in Appendix A.1 via substituting ν 2 with ν 2 + σ 2 t at time step t.First, by Proposition 2, suppose that the sample is initialized in the distribution P (0) σ0 , then with probability at least 1 -exp(-Ω(d)), we haveEQUATIONThen, with the assumption that the initialization satisfies ∥n 0 ∥ 2 ≥ 3ν 2 0 +ν 2 max +4σ 2 0 4 d, the following proposition similar to Proposition 3 shows that ∥n t ∥ remains large with high probability.Proposition 4. Consider a data distribution P satisfies the constraints specified in Theorem 2.We follow annealed Langevin dynamics for T = exp(O(d)) steps with noise level c σ ≥ σ 0 ≥σ 1 ≥ σ 2 ≥ • • • ≥ σ T ≥ 0 for some constant c σ > 0.Suppose that the initial sample satisfies∥n 0 ∥ 2 ≥ 3ν 2 0 +ν 2 max +4σ 2 0 4 d, then with probability at least 1 -T • exp(-Ω(d)), we have that ∥n t ∥ 2 > ν 2 0 +ν 2 max +2σ 2 t 2 d for all t ∈ {0} ∪ [T ].Proof of Proposition 4. We prove Proposition 4 by induction. Suppose the theorem holds for all T values of 1, • • • , T -1. We consider the following 3 cases:• If there exists some t ∈ [T ] such that δ t > ν 2 0 + σ 2 t , by Lemma 3 we know that with probability at least 1 -exp(-Ω(d)), we have∥n t ∥ 2 ≥ 3(ν 2 0 +σ 2 t )+(ν 2 max +σ 2 t ) 4 d = 3ν 2 0 +ν 2 max +4σ 2 t 4 d,thus the problem reduces to the two sub-arrays n 0 , • • • , n t-1 and n t , • • • , n T , which can be solved by induction. = nt-1 - δ t 2ν 2 0 + 2σ 2 t nt-1 + δ t ϵ (n) t .Since ν 0 > ν i and c σ ≥ σ t for all t ∈ {0} ∪ [T ], we haveν 2 i +c 2 σ ν 2 0 +c 2 σ ≥ ν 2 i +σ 2 t ν 2 0 +σ 2 t . Notice that for function f (z) = log z -z 2 + 1 2z , we have d dz f (z) = 1 z -1 2 -1 2z 2 = -1 2 1 z -1 2 ≤ 0.Thus, by the assumption∥µ i -µ 0 ∥ 2 ≤ ν 2 0 -ν 2 i 2 log ν 2 i + c 2 σ ν 2 0 + c 2 σ - ν 2 i + c 2 σ 2ν 2 0 + c 2 σ + ν 2 0 + c 2 σ 2ν 2 i + c 2 σ d,we have that for all t ∈ [T ],∥µ i -µ 0 ∥ 2 ≤ ν 2 0 -ν 2 i 2 log ν 2 i + σ 2 t ν 2 0 + σ 2 t - ν 2 i + σ 2 t 2ν 2 0 + σ 2 t + ν 2 0 + σ 2 t 2ν 2 i + σ 2 t d. Conditioned on ∥n t-1 ∥ 2 > ν 2 0 +ν 2 max +2σ 2 t-1 2 d for all t ∈ [T ],by Lemma 6 we have that forT = exp(O(d)), ∥n T -n T ∥ < 5ν 2 0 + 3ν 2 max + 8σ 2 T 8 - ν 2 0 + ν 2 max + 2σ 2 T 2 √ d.By Lemma 8 we have that with probability at least 1 -exp(-Ω(d)),∥n T ∥ 2 ≥ 5ν 2 0 + 3ν 2 max + 8σ 2 T 8 d.Combining the two inequalities implies the desired bound∥n T ∥ ≥ ∥n T ∥ -∥n T -n T ∥ > ν 2 0 + ν 2 max + 2σ 2 T 2 d.Hence by induction we obtain∥n t ∥ 2 > ν 2 0 +ν 2 max +2σ 2 t 2d for all t ∈ {0} ∪ [T ] with probability at least(1 -(T -1) exp(-Ω(d))) • (1 -exp(-Ω(d))) ≥ 1 -T exp(-Ω(d)).Therefore we complete the proof of Proposition 4.Finally, combining equation 11 and Proposition 4 finishes the proof of Theorem 2.The proof framework is similar to the proof of Theorem 1. To begin with, we validate Assumption 2.v. in the following lemma:Lemma 9. For constants ν 0 , ν i , c ν , c L satisfying Assumptions 2.iii. and 2.iv., we have(1-cν )ν 2 0 -ν 2 i 2(1-cν ) > 0 and log cν ν 2 i (c 2 L +cν c L )ν 2 0 - ν 2 i 2(1-cν )ν 2 0 + (1-cν )ν 2 0 2ν 2 i > 0 are both positive constants. Proof of Lemma 9. From Assumption 2.iv. that ν 2 0 > ν 2 max 1-cν ≥ ν 2 i 1-cν , we easily obtain (1-cν )ν 2 0 -ν 2 i 2(1-cν )> 0 is a positive constant. For the second property, let f (z) := logcν ν 2 i (c 2 L +cν c L )z - ν 2 i 2(1-cν )z + (1-cν )z 2ν 2 i .For any z >ν 2 i 1-cν , the derivative of f (z) satisfies d dz f (z) = - 1 z + ν 2 i 2(1 -c ν )z 2 + 1 -c ν 2ν 2 i = ν 2 i 2(1 -c ν ) 1 -c ν ν 2 i - 1 z 2 > 0.Therefore, when4(c 2 L +cν c L ) cν (1-cν ) ≤ 1, we have f (ν 2 0 ) > f ν 2 i 1 -c ν = log c ν (1 -c ν ) c 2 L + c ν c L ≥ log 4 > 0. When 4(c 2 L +cν c L ) cν (1-cν ) > 1, we have f (ν 2 0 ) > f 4(c 2 L + c ν c L ) c ν (1 -c ν ) ν 2 i 1 -c ν = 2 log c ν (1 -c ν ) 2(c 2 L + c ν c L ) - c ν (1 -c ν ) 8(c 2 L + c ν c L ) + 2(c 2 L + c ν c L ) c ν (1 -c ν ) ≥ 2 -2 log 2 - 2(c 2 L + c ν c L ) c ν (1 -c ν ) - c ν (1 -c ν ) 8(c 2 L + c ν c L ) + 2(c 2 L + c ν c L ) c ν (1 -c ν ) > 2 -2 log 2 - 1 2 > 0.Thus we obtain Lemma 9.Without loss of generality, we assume µ 0 = 0 d . Similar to the proof of Theorem 1, we decomposex t = Rr t + Nn t , andϵ t = Rϵ (r) t + Nϵ (n)t , where R ∈ R d×r an orthonormal basis of the vector space {µ i } i∈ [k] and N ∈ R d×n an orthonormal basis of the null space of{µ i } i∈[k] . To show ∥x t -µ i ∥ 2 > ν 2 0 2 + ν 2 max 2(1-cν ) d, it suffices to prove ∥n t ∥ 2 > ν 2 0 2 + ν 2 max 2(1-cν ) d.By Proposition 2, if x 0 is initialized in the distribution P (0) , i.e.,x 0 ∼ P (0) , since ν 2 0 > 1 1-cν ν 2 max , with probability at least 1 -exp(-Ω(d)) we haveEQUATIONThen, conditioned on Proof of Proposition 5. Firstly, by Lemma 3, if δ t > ν 2 0 , since ν 2 0 > ν 2 max 1-cν , we similarly have that∥n 0 ∥ 2 ≥ 3ν 2 0 4 + ν 2 max 4(1-cν ) d,∥n t ∥ 2 ≥ 3ν 2 0 4 + ν 2 max4(1-cν ) d with probability at least 1 -exp(-Ω(d)) regardless of the previous state x t-1 . We then consider the case when δ t ≤ ν 2 0 . Intuitively, we aim to prove that the scorefunction is close to -x ν 2 0 when ∥n∥ 2 ≥ ν 2 0 2 + ν 2 max 2(1-cν ) d.Towards this goal, we first show thatP (0) (x) is exponentially larger than P (i) (x) for all i ∈ [k]in the following lemma:Lemma 10. Suppose P satisfies Assumption 2. Then for any ∥n∥ 2 ≥ν 2 0 2 + ν 2 max2(1-cν ) d, we haveP (i) (x) P (0) (x) ≤ exp(-Ω(d))and∥∇xP (i) (x)∥ P (x) ≤ exp(-Ω(d)) for all i ∈ [k].Proof of Lemma 10. We first give an upper bound on the sub-Gaussian probability density. For any vector v ∈ R d , by considering some vector m ∈ R d , from Markov's inequality and the definition in equation 4 we can boundP z∼P (i) m T (z -µ i ) ≥ m T (v -µ i ) ≤ E z∼P (i) exp m T (z -µ i ) exp (m T (v -µ i )) ≤ exp ν 2 i ∥m∥ 2 2 -m T (v -µ i ) .Upon optimizing the last term at m = v-µi ν 2 i , we obtainEQUATION) Denote B := z : (v -µ i ) T (v -z) ≤ 0 . To bound P z∼P (i) (z ∈ B), we first note that log P (i) (v) -log P (i) (z) = 1 0 ⟨v -z, ∇ log P (i) (v + λ(z -v))⟩ dλ = ⟨v -z, ∇ log P (i) (v)⟩ + 1 0 ⟨v -z, ∇ log P (i) (v + λ(z -v)) -∇ log P (i) (v)⟩ dλ ≤ ∥v -z∥ ∇ log P (i) (v) + 1 0 ∥v -z∥ ∇ log P (i) (v + λ(z -v)) -∇ log P (i) (v) dλ ≤ ∥v -z∥ • L i ∥v -µ i ∥ + 1 0 ∥v -z∥ • L i ∥λ(z -v)∥ dλ (14) ≤ L i c ν 2c L ∥v -µ i ∥ 2 + c L + c ν 2c ν L i ∥v -z∥ 2 ,where equation 14 follows from Assumption 2.ii. that ∇ log P (i) (µ i ) = 0 d and Assumption 2.iii.that the score function ∇ log P (i) is L i -Lipschitz. Therefore we obtainEQUATION)By observing that g :EQUATION)Hence, by combining equation 13, equation 15, and equation 16, we obtainexp - ∥v -µ i ∥ 2 2ν 2 i ≥ P z∼P (i) (v -µ i ) T (v -z) ≤ 0 ≥ P (i) (v) exp - L i c ν 2c L ∥v -µ i ∥ 2 • 1 2 2πc ν (c L + c ν )L i d 2 . By Assumption 2.iii. that L i ≤ c L ν 2 iwe obtain the following bound on the probability density:EQUATION)Then we can bound the ratio of P (i) and P (0) . For all i ∈ [k], define ρ i (x) := P (i) (x) P (0) (x) , then we haveρ i (x) = P (i) (x) P (0) (x) ≤ 2(2πc ν ν 2 i /(c 2 L + c ν c L )) -d/2 exp -(1 -c ν ) ∥x -µ i ∥ 2 /2ν 2 i (2πν 2 0 ) -d/2 exp -∥x∥ 2 /2ν 2 0 = 2 (c 2 L + c ν c L )ν 2 0 c ν ν 2 i d 2 exp ∥x∥ 2 2ν 2 0 - (1 -c ν ) ∥x -µ i ∥ 2 2ν 2 i = 2 (c 2 L + c ν c L )ν 2 0 c ν ν 2 i d 2 exp 1 2ν 2 0 - 1 -c ν 2ν 2 i ∥Nn∥ 2 + ∥Rr∥ 2 2ν 2 0 - (1 -c ν ) ∥Rr -µ i ∥ 2 2ν 2 i = 2 (c 2 L + c ν c L )ν 2 0 c ν ν 2 i d 2 exp 1 2ν 2 0 - 1 -c ν 2ν 2 i ∥n∥ 2 + ∥r∥ 2 2ν 2 0 - (1 -c ν ) r -R T µ i 2 2ν 2 i ,where the last step follows from the definition that R ∈ R d×r an orthogonal basis of the vector space{µ i } i∈[k] and N T N = I n . Since ν 2 i < (1 -c ν )ν 2 0 , the quadratic term ∥r∥ 2 2ν 2 0 - (1-cν )∥r-R T µi∥ 2 2ν 2 i is maximized at r = (1-cν )ν 2 0 R T µi (1-cν )ν 2 0 -ν 2 i. Therefore, we obtain∥r∥ 2 2ν 2 0 - (1 -c ν ) r -R T µ i 2 2ν 2 i ≤ (1 -c ν ) ∥µ i ∥ 2 2((1 -c ν )ν 2 0 -ν 2 i ) . Hence, for ∥µ i -µ 0 ∥ 2 ≤ (1-cν )ν 2 0 -ν 2 i 2(1-cν ) log cν ν 2 i (c 2 L +cν c L )ν 2 0 - ν 2 i 2(1-cν )ν 2 0 + (1-cν )ν 2 0 2ν 2 i d and ∥n∥ 2 ≥ ν 2 0 2 + ν 2 max 2(1-cν ) d, we have ρ i (x) ≤ 2 (c 2 L + c ν c L )ν 2 0 c ν ν 2 i d 2 exp 1 2ν 2 0 - 1 -c ν 2ν 2 i ∥n∥ 2 + (1 -c ν ) ∥µ i ∥ 2 2((1 -c ν )ν 2 0 -ν 2 i ) ≤ 2 (c 2 L + c ν c L )ν 2 0 c ν ν 2 i d 2 exp 1 2ν 2 0 - 1 -c ν 2ν 2 i ν 2 0 2 + ν 2 i 2(1 -c ν ) d + (1 -c ν ) ∥µ i ∥ 2 2((1 -c ν )ν 2 0 -ν 2 i ) = 2 exp -log c ν ν 2 i (c 2 L + c ν c L )ν 2 0 - ν 2 i 2(1 -c ν )ν 2 0 + (1 -c ν )ν 2 0 2ν 2 i d 2 + (1 -c ν ) ∥µ i ∥ 2 2((1 -c ν )ν 2 0 -ν 2 i ) ≤ 2 exp -log c ν ν 2 i (c 2 L + c ν c L )ν 2 0 - ν 2 i 2(1 -c ν )ν 2 0 + (1 -c ν )ν 2 0 2ν 2 i d 4 .From Lemma 9, we obtain ρ i (x) ≤ exp(-Ω(d)).To show∥∇xP (i) (x)∥ P (x)≤ exp(-Ω(d)), from Assumptions 2.ii. and 2.iii. we have∇ x P (i) (x) P (i) (x) = ∇ x P (i) (x) P (i) (x) - ∇ x P (i) (µ i ) P (i) (µ i ) = ∇ x log P (i) (x) -∇ x log P (i) (µ i ) ≤ L i ∥x -µ i ∥ ≤ c L ν 2 i ∥x -µ i ∥ .Therefore, we can bound∥∇xP (i) (x)∥ P (x) ≤ c L ν 2 i ρ i (x) ∥x -µ i ∥. When ∥x -µ i ∥ = exp(o(d)) is small, by ρ i (x) ≤ exp(-Ω(d)) we directly have ∥∇xP (i) (x)∥ P (x) ≤ exp(-Ω(d)). When ∥x -µ i ∥ = exp(Ω(d)) is exceedingly large, from equation 17 we have∇ x P (i) (x) P (x) ≤ 2c L ν 2 i (c 2 L + c ν c L )ν 2 0 c ν ν 2 i d 2 exp ∥x∥ 2 2ν 2 0 - (1 -c ν ) ∥x -µ i ∥ 2 2ν 2 i ∥x -µ i ∥ . Since ν 2 0 > ν 2 i 1-cν , when ∥x -µ i ∥ = exp(Ω(d)) ≫ ∥µ i ∥ we have exp ∥x∥ 2 2ν 2 0 - (1 -c ν ) ∥x -µ i ∥ 2 2ν 2 i = exp(-Ω(∥x -µ i ∥ 2 )).Therefore∥∇xP (i) (x)∥ P (x)≤ exp(-Ω(d)). Thus we complete the proof of Lemma 10.Similar to Lemma 5, the following lemma proves that when the previous state n t-1 is far from a mode, a single step of Langevin dynamics with bounded step size is not enough to find the mode.Lemma 11. Suppose δ t ≤ ν 2 0 and ∥n t-1 ∥ 2 > 36ν 2 0 d, then we have ∥n t ∥ 2 ≥ ν 2 0 d with probability at least 1 -exp(-Ω(d)).Proof of Lemma 11. For simplicity, denote v := n t-1 + δt 2 N T ∇ x log P (x t-1 ). Since P = k i=0 w i P (i) and P (0) = N (µ 0 , ν 2 0 I d ), the score function can be written asEQUATIONFor ∥n t-1 ∥ 2 > 36ν 2 0 d by Lemma 10 we have∥∇xP (i) (xt-1)∥ P (xt-1) ≤ exp(-Ω(d)). Since δ t ≤ ν 2 0 , we can bound the norm of v by ∥v∥ = n t-1 + δ t 2 N T ∇ x log P (x t-1 ) = n t-1 - δ t 2ν 2 0 n t-1 + i∈[k] w i δ t 2ν 2 0 P (i) (x t-1 ) P (x t-1 ) n t-1 + i∈[k] w i δ t 2 N T ∇ x P (i) (x t-1 ) P (x t-1 ) ≥   1 - δ t 2ν 2 0 + i∈[k] w i δ t 2ν 2 0 P (i) (x t-1 ) P (x t-1 )   n t-1 - i∈[k] w i δ t 2 ∇ x P (i) (x t-1 ) P (x t-1 ) ≥ 1 2 ∥n t-1 ∥ - i∈[k] w i δ t 2 exp(-Ω(d)) > 2ν 0 √ d.On the other hand, from ϵEQUATIONt ⟩ ∥v∥ ∼ N (0, 1) for any fixed v ̸ = 0 n , hence by Lemma 2 we haveEQUATIONCombining the above inequalities gives∥n t ∥ 2 = v + δ t ϵ (n) t 2 ≥ ∥v∥ 2 -2ν 0 |⟨v, ϵ (n) t ⟩| ≥ ∥v∥ 2 - ν 0 √ d 2 ∥v∥ > ν 2 0 d with probability at least 1 -8 √ 2πd exp -d 32 = 1 -exp(-Ω(d)). This proves Lemma 11.When ∥n t-1 ∥ 2 ≤ 36ν 2 0 d, similar to Theorem 1, we consider a surrogate recursion nt such that n0 = n 0 and for all t ≥ 1,EQUATIONThe following Lemma shows that nt is sufficiently close to the original recursion n t .Lemma 12. For any t ≥ 1, given that for all j ∈ [t], δ j ≤ ν 2 0 andν 2 0 2 + ν 2 max 2(1-cν ) d ≤ ∥n j-1 ∥ 2 ≤ 36ν 2 0 d, if µ i satisfies Assumption 2.v. for all i ∈ [k], we have ∥n t -n t ∥ ≤ t exp(Ω(d)) √ d.Proof of Lemma 12. By equation 18 we have that for all j ∈ [t],∥n j -n j ∥ = nj-1 -n j-1 - δ j 2ν 2 0 nj-1 - δ j 2 N T ∇ x log P (x j-1 ) = nj-1 -n j-1 - i∈[k] w i P (i) (x j-1 ) ν 2 0 P (x j-1 ) n j-1 - i∈[k] w i N T ∇ x P (i) (x j-1 ) P (x j-1 ) ≤ ∥n j-1 -n j-1 ∥ + i∈[k] w i P (i) (x j-1 ) ν 2 0 P (x j-1 ) ∥n j-1 ∥ + i∈[k]w i ∇ x P (i) (x j-1 ) P (x j-1 ) .By Lemma 10, we haveP (i) (xj-1) P (0) (xj-1) ≤ exp(-Ω(d))and∥∇xP (i) (xj-1)∥ P (xj-1) ≤ exp(-Ω(d)) for all i ∈ [k], hence from ∥n j-1 ∥ ≤ 6ν 0 √ d we obtain a recursive bound ∥n j -n j ∥ ≤ ∥n j-1 -n j-1 ∥ + 1 exp(Ω(d)) √ d.Finally, by n0 = n 0 , we have∥n t -n t ∥ = j∈[t] (∥n j -n j ∥ -∥n j-1 -n j-1 ∥) ≤ t exp(Ω(d)) √ d.Hence we obtain Lemma 12.Armed with the above lemmas, we are now ready to establish Proposition 5 by induction. Please note that we also apply some lemmas from the proof of Theorem 1 by substituting ν 2 max withν 2 max 1-cν .Suppose the theorem holds for all T values of 1, • • • , T -1. We consider the following 3 cases:• If there exists some t ∈ [T ] such that δ t > ν 2 0 , by Lemma 3 we know that with probability at least 1 -exp(-Ω(d)), we have∥n t ∥ 2 ≥ 3ν 2 0 4 + ν 2 max 4(1-cν ) d, thus the problem reduces to the two sub-arrays n 0 , • • • , n t-1 and n t , • • • , n T , which can be solved by induction.• Supposeδ t ≤ ν 2 0 for all t ∈ [T ]. If there exists some t ∈ [T ] such that ∥n t-1 ∥ 2 > 36ν 2 0 d, byLemma 11 we know that with probability at least 1 -exp(-Ω(d)), we have∥n t ∥ 2 ≥ ν 2 0 d > 3ν 2 0 4 + ν 2 max 4(1-cν ) d, thus the problem similarly reduces to the two sub-arrays n 0 , • • • , n t-1and n t , • • • , n T , which can be solved by induction.• Supposeδ t ≤ ν 2 0 and ∥n t-1 ∥ 2 ≤ 36ν 2 0 d for all t ∈ [T ]. Conditioned on ∥n t-1 ∥ 2 > ν 2 0 2 + ν 2 max 2(1-cν ) d for all t ∈ [T ],by Lemma 12 we have that forT = exp(O(d)), ∥n T -n T ∥ < 5ν 2 0 8 + 3ν 2 max 8(1 -c ν ) - ν 2 0 2 + ν 2 max 2(1 -c ν ) √ d.By Lemma 8 we have that with probability at least 1 -exp(-Ω(d)),∥n T ∥ 2 ≥ 5ν 2 0 8 + 3ν 2 max 8(1 -c ν ) d. ∥n T ∥ ≥ ∥n T ∥ -∥n T -n T ∥ > ν 2 0 2 + ν 2 max 2(1 -c ν ) d.Hence by induction we obtain∥n t ∥ 2 > ν 2 0 2 + ν 2 max 2(1-cν ) d for all t ∈ [T ] with probability at least (1 -(T -1) exp(-Ω(d))) • (1 -exp(-Ω(d))) ≥ 1 -T exp(-Ω(d)).Therefore we complete the proof of Proposition 5. Suppose that P (0) = N (µ 0 , ν 2 0 I d ) is Gaussian and for all i ∈ [k], P (i) satisfies i. P (i) is a sub-Gaussian distribution of mean µ i with parameter ν 2 i ,ii. P (i) is differentiable and∇P (i) σt (µ i ) = 0 d for all t ∈ {0} ∪ [T ],iii. for all t ∈ {0} ∪ [T ], the score function ofP (i) σt is L i,t -Lipschitz such that L i,t ≤ c L ν 2 i +σ 2 t for some constant c L > 0, iv. ν 2 0 > max 1, 4(c 2 L +cν c L ) cν (1-cν ) ν 2 max +c 2 σ 1-cν -c 2 σ for constant c ν ∈ (0, 1), where ν max := max i∈[k] ν i , v. ∥µ i -µ 0 ∥ 2 ≤ (1-cν )ν 2 0 -ν 2 i -cν c 2 σ 2(1-cν ) log cν (ν 2 i +c 2 σ ) (c 2 L +cν c L )(ν 2 0 +c 2 σ ) - (ν 2 i +c 2 σ ) 2(1-cν )(ν 2 0 +c 2 σ ) + (1-cν )(ν 2 0 +c 2 σ ) 2(ν 2 i +c 2 σ )d.The feasibility of Assumption 3.v. can be validated by substituting ν 2 in Lemma 9 with ν 2 + c 2 σ .To establish Theorem 4, we first note from Proposition 1 that for a sub-Gaussian mixture P = k i=0 w i P (i) , the perturbed distribution of noise level σ isP σ = k i=0 w i P (i)σ , where P (0) = N (µ 0 , (ν 2 i +σ 2 )I d ) and P (i) is a sub-Gaussian distribution with mean µ i and sub-Gaussian parameter (ν 2 i + σ 2 ). Similar to the proof of Theorem 1, we decomposex t = Rr t + Nn t , andϵ t = Rϵ (r) t + Nϵ (n)t , where R ∈ R d×r an orthonormal basis of the vector space {µ i } i∈ [k] and N ∈ R d×n an orthonormal basis of the null space of {µ i } i∈ [k] . Now, we prove Theorem 4 by applying the techniques developed in Appendix A.1 and A.3 via substituting ν 2 and ν 2 1-cν with ν 2 +σ 2 t 1-cν at time step t. Note that for allt ∈ {0} ∪ [T ], Assumption 3.iv. implies ν 2 0 + σ 2 t > max 1, 4(c 2 L +cν c L ) cν (1-cν ) ν 2 max +σ 2 t 1-cν because c σ ≥ σ t .First, by Proposition 2, suppose that the sample is initialized in the distribution P (0) σ0 , then with probability at least 1 -exp(-Ω(d)), we haveEQUATIONThen, with the assumption that the initialization satisfies ∥n 0 ∥ 2 ≥ 3(ν 2 0 +σ 2 0 ) 4 + ν 2 max +σ 2 0 4(1-cν )nt = nt-1 - δ t 2ν 2 0 + 2σ 2 t nt-1 + δ t ϵ (n) t .Since ν 0 > ν i and c σ ≥ σ t for all t ∈ {0} ∪ [T ], we haveν 2 i +c 2 σ ν 2 0 +c 2 σ > ν 2 i +σ 2 t ν 2 0 +σ 2 t . Notice that for function f (z) = log z -z 2 + 1 2z , we have d dz f (z) = 1 z -1 2 -1 2z 2 = -1 2 1 z -1 2 ≤ 0.Thus, by Assumption 3.v. we have that for all t ∈ [T ],∥µ i -µ 0 ∥ 2 ≤ (1 -c ν )ν 2 0 -ν 2 i -c ν c 2 σ 2(1 -c ν ) log c ν (ν 2 i + c 2 σ ) (c 2 L + c ν c L )(ν 2 0 + c 2 σ ) - (ν 2 i + c 2 σ ) 2(1 -c ν )(ν 2 0 + c 2 σ ) + (1 -c ν )(ν 2 0 + c 2 σ ) 2(ν 2 i + c 2 σ ) d ≤ (1 -c ν )ν 2 0 -ν 2 i -c ν σ 2 t 2(1 -c ν ) log c ν (ν 2 i + σ 2 t ) (c 2 L + c ν c L )(ν 2 0 + σ 2 t ) - (ν 2 i + σ 2 t ) 2(1 -c ν )(ν 2 0 + σ 2 t ) + (1 -c ν )(ν 2 0 + σ 2 t ) 2(ν 2 i + σ 2 t ) d Conditioned on ∥n t-1 ∥ 2 > ν 2 0 +σ 2 t-1 2 + ν 2 max +σ 2 t-1 2(1-cν )d for all t ∈ [T ], by Lemma 12 we have that for T = exp(O(d)),∥n T -n T ∥ < 5(ν 2 0 + σ 2 T ) 8 + 3(ν 2 max + σ 2 T ) 8(1 -c ν ) - ν 2 0 + σ 2 T 2 + ν 2 max + σ 2 T 2(1 -c ν ) √ d.By Lemma 8 we have that with probability at least 1 -exp(-Ω(d)),∥n T ∥ 2 ≥ 5(ν 2 0 + σ 2 T ) 8 + 3(ν 2 max + σ 2 T ) 8(1 -c ν ) d.Combining the two inequalities implies the desired bound∥n T ∥ ≥ ∥n T ∥ -∥n T -n T ∥ > ν 2 0 + σ 2 T 2 + ν 2 max + σ 2 T 2(1 -c ν ) d.Hence by induction we obtain∥n t ∥ 2 > ν 2 0 +σ 2 T 2 + ν 2 max +σ 2 T 2(1-cν )d for all t ∈ [T ] with probability at least(1 -(T -1) exp(-Ω(d))) • (1 -exp(-Ω(d))) ≥ 1 -T exp(-Ω(d)).Therefore we complete the proof of Proposition 6.Finally, combining equation 20 and Proposition 6 finishes the proof of Theorem 4. For the score networks of vanilla and annealed Langevin dynamics, following from Song and Ermon (2019) , we use the 4-cascaded RefineNet (Lin et al., 2017) , a modern variant of U-Net (Ronneberger et al., 2015) with residual design. For the score networks of chained Langevin dynamics, we use the official PyTorch implementation of an LSTM network (Sak et al., 2014) followed by a linear layer.For MNIST and Fashion-MNIST datasets, we set the input size of the LSTM as Q = 14, the number of features in the hidden state as 1024, and the number of recurrent layers as 2. The inputs of LSTM include inputting tensor, hidden state, and cell state, and the outputs of LSTM include the next hidden state and cell state, which can be fed to the next input. To estimate the noisy score function, we first input the noise level σ (repeated for Q times to match the input size of LSTM) and all-0 hidden and This paper presents work whose goal is to advance the field of machine learning. No potential societal consequence of this work needs to be highlighted here. Justification: We list the paper's contributions about the mode-seeking tendencies of vanilla, annealed, and chained Langevin Dynamics in the abstract and at the end of Section 1. The scope of this work is also discussed in Section 1.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper.• The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.• The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.• It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 


Decision: Although the reviewers have identified some positive aspects of this work, I have ultimately decided to recommend rejection. My main concern is that this work is primarily a theory paper, but it does not have sufficient technical novelty to be a good NeurIPS theory paper. 

- I think from a generative modeling perspective, the chained Langevin dynamics is not a particularly novel idea; it's just adding an autoregressive layer to the generation process. 
- From a theoretical standpoint, the only analysis given of the chained Langevin dynamics (theorem 5) is basically just the triangle inequality. 
- In the main theoretical contribution, authors give some examples of multimodal settings where "Langevin Dynamics is unlikely to find all mixture components within a sub-exponential number of steps in the data dimension". At a high level, this is already a very well-known and well understood phenomena and the authors even cite some works on metastability which have clearly already demonstrated this high level point. I do not think the example given by the authors adds something truly new to the literature.  
- The failure of simulated annealing in simple multimodal settings would probably be a more interesting phenomena and the authors made some efforts to extend in this direction. But we already know from the theoretical guarantees for diffusion models that if the setup is done properly then a natural annealed sampler does provably work [e.g. from CCLYZ '23]... (The authors should, I think, at least seriously discuss this in the next version of the work.)

I encourage the authors to revise the work and resubmit to a future venue. Perhaps a more serious investigation of the chained sampler could be done.

Reviews: review: paper_summary: paper_summary: This paper studies Langevin-based algorithms for sampling from multimodal distributions, motivated by generative modeling. The main content of the paper are lower bounds on the convergence of both Langevin and annealed Langevin for mixtures of Gaussian and sub-Gaussian distributions, as well as a proposed modification of the annealed Langevin dynamics to operate on coordinate patches one-at-a-time.
strength_weakness: strength_weakness: - The lower bounds hold only for the distance between the sample and the mean, rather than any standard notion of distance or divergence between probability measures. Moreover, I do not expect that these bounds imply such a quantity is large.

- Related to the above point, it is difficult to appreciate the significance of the lower bound since the lower bound does not depend on the separation between the means. In particular, it seems the lower bounds only show that the iterate remains roughly on the order of the larger variance which is, for example, not surprising in the case where the variances are all of the same order.

- The hidden constants in the $\Omega$ notation are important but difficult to find (as they are suppressed in the main text and some of the appendix). In particular, there should be dependence on the mixture weights but this can't be seen from their result.

- It is unclear if the upper bound in Theorem 5 can be instantiated for their algorithm (see question below).
questions: questions: - In light of [1], there is arguably "no mystery" in terms of convergence of the reversed Langevin dynamics in diffusion models: as long as the score function can be accurately estimated, the reversed dynamics will converge to the target. At a high level, why do you then study annealed Langevin in the setting where the scores are known exactly?

- I suggest you stop using the $\Omega$ notation and make the hidden constants more clear as they are very important for understanding your result.

- With regard to Theorem 5, a remark that explains, even if conjecturally, how the run-time of the Langevin algorithm might scale in the dimension and how the bound in Theorem 5 does, or does not, imply that your method could be successful would be greatly helpful.

- On page 2 you write "Regarding discrete SGLD, Lee et al. (2018) constructed a probability distribution whose density is close to a mixture of two well-separated isotropic Gaussians, and proved that SGLD could not find one of the two modes within an exponential number of steps." However, I was not able to find this lower bound in Lee et al. (2018). Which result specifically are you referring to?

[1] Chen, S., Chewi, S., Li, J., Li, Y., Salim, A., and Zhang, A. R. (2023). Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. In International Conference on Learning Representations
limitations: limitations: The limitations of the work have been adequately addressed.
paper_summary: paper_summary: A new algorithm is proposed, called Chained Langevin Dynamics, to improve on the mode-seeking properties of Langevin Dynamics, after annleade Langevin Dynamics had been proposed but did not give significant improvements.
Results about the mode-seeking properties of the three algorithms are obtained.
The results of numerical experiments on synthetic and real image datasets are also shown.
strength_weakness: strength_weakness: No evident connection has been established between experiments and mathematical results.
The description/comment of experiments could have been more accurate (see Questions).
questions: questions: In structured data, does the order of patches matter?

How large is the selected size Q of the patches in the examples? Why was it selected that way? Smaller patches help convergence due to the reduced dimension. What happens for too small patches? Does the algorithm still work?

"Regarding the neural network architecture of the score function estimator, for vanilla and annealed Langevin dynamics we use U-Net (Ronneberger et al., 2015) following from Song and Ermon (2019). For chained Langevin dynamics, we proposed to use Recurrent Neural Network (RNN) architectures." The change in the set up seems a little arbitrary. What did it happen by using U-Net for chained Langevin dynamics or, viceversa, using the RNN for vanilla and annealed Langevin dynamics?

More on the general concepts of the paper, the authors might find interesting that the principle of chained Langevin dynamics seems based on "nucleation" of different modes and spreading of the information of the randomly selected mode through the conditional probability to the entire image, much like the mechanism described in this old paper that allowed reconstruction of images from very little sampling 
Statistical-Physics-Based Reconstruction in Compressed Sensing
F. Krzakala, M. Mézard, F. Sausset, Y. F. Sun, and L. Zdeborová
Phys. Rev. X 2, 021005 – Published 11 May 2012
limitations: limitations: There is a section about limitations in the text.
paper_summary: paper_summary: The authors study Langevin dynamics (as well as its annealed counterpart) for gaussian mixtures and sub-gaussian mixtures. In Sec. 4, they prove that Langevin remains stuck in the "dominant mode" for an at least exponential time, a claim that is often made in the ML literature but which is never formally proved. In Sec. 5, they provide a sequential method to get rid of this dependence.
strength_weakness: strength_weakness: I do not understand why it is sensible to say that initially, $p_0$ should follow $P_0$, one of the component of the mixture, isn't it a rather strong assumption?
questions: questions: Do similar results hold when $p_0$ is initialized at $\mathcal{N}(0, Id)$ and this component is not in the target mixture?
limitations: limitations: It seems like the assumption $p_0 \sim P_0$ is not enough justified. Also, I would have like an insight of the proof of the Theorems in Sec. 4.
paper_summary: paper_summary: The authors consider the Langevin process for sampling from a target distribution $\pi$. This process is known to be slow-converging for multimodal targets: in practice, it has been observed that the process gets "stuck" in some modes of the target, and do not "reach" other modes of the target. The authors provide theoretical results for this behavior. In Theorem 1, they prove that by evolving a particle with the Langevin process during exponential time (in the dimension), the particle will still be far away (in probability) from some modes. They also prove, in Theorem 2, that this negative result holds even when using the popular heuristic of "annealing" the Langevin process using intermediate distributions, obtained by adding different levels of Gaussian noise to the target samples. 

Instead, the authors propose running an alternative sampling process which they call "Chained Langevin dynamics". This consists in running "annealed" Langevin processes for each component of the target distribution, that is, for each $\pi(x_i | x_{-i})$. The authors estimate the score of each of these conditional targets using a score-matching loss, and empirically demonstrate the ability of their process to reach the different target modes in a limited time. Theoretically, they prove their process approximates the target (in TV divergence) in linear time (in the dimension).
strength_weakness: strength_weakness: See questions.
questions: questions: **Q1**: Theorem 5 is stated at the distribution-level (the result is for $p_t$), while Theorems 1 and 2 are stated at the particle-level (the result is for $x_t \sim p_t$ ). The authors hint that Theorems 1 and 2 could also be formulated at the distribution-level: "this notion can also easily be translated into a lower bound in terms of other distance measures such as total variation distance and Wasserstein 2-distance". If this is an easy result to add, could the authors add it to their paper? Lower-bounds are notoriously hard(er than upper bounds) on the convergence. Namely, stating Theorem 2 as a lower bound on Annealed Langevin Dynamics would be an important contribution. 

**Q2**: Theorem 2 is surprising to me. It seems to convey that Langevin Dynamics are *worse* when annealing, in the sense that the particle $x_t$ is *further* away by the additive constant $2 \sigma_t^2$ from the non-dominant modes $\mu_i$. This seems to be in contradiction with theoretical and empirical results that show that annealing helps reach all the modes of the target distribution, as conveyed in Figure 1 of [1] or Figure 3 of [2] or even Figure [3] of the authors' submission where the Annealed Langevin process produces more accurate samples than the Vanilla Langevin process. Can the authors explain why that is?

**Q3**: When introducing the Chained Langevin Dynamics, the authors seem to motivate this process by saying that by sampling one component at a time from $\pi(x_i | x_{-i})$, this reduces the complexity in the dimension. Could the authors elaborate on this? My understand is that if each of the one-dimensional, conditional targets $\pi(x_i | x_{-i})$ were heavily multimodal, then Annealed Langevin dynamics would still struggle to cover all the modes, although it would struggle less than in higher dimensions in light of Theorem 2 taking $d >> 1$.

[1] Dynamical Regimes of Diffusion Models. Biroli et al. Arxiv, 2024. 

[2] Generative Modeling by Estimating Gradients of the Data Distribution. Song et al. NeurIPS 2019.
limitations: limitations: Yes.

Comments: title: Thanks for your feedback
Thanks for your feedback
Thanks for your feedback
Thanks for your thorough response
Answer to authors
Response to rebuttal
Thanks for your feedback
Answer to the authors
comment: We sincerely thank Reviewer nEdB for the thorough summarization and constructive feedback on our response. The answer is affirmative to the raised question "Do you have these lower bounds for **both** Theorems 1 and 2?". For both vanilla and annealed Langevin dynamics (with bounded noise levels), Theorems 1 and 2 can be translated into lower bounds in total variation distance, as described in global rebuttal #1. We will clarify this point in the text.
We thank Reviewer Ave7 for his/her feedback on our response. First, we would like to clarify that our response was based on a multivariate $d$-dimensional $\mathbf{x}_0 \sim \mathcal{N}(\mathbf{0}_d,\mathbf{I}_d)$, where the complexity analysis is in terms of dimension $d$. Therefore, our analysis does not focus on the univariate case $\mathcal{N}(0,1)$ mentioned in the reviewer’s question.

Next, we would like to clarify the assumption of the initial sample being far from the vector space of $\\left\\{ \mathbf{\mu}\_i \\right\\}_{i \in [k]}$. This assumption only requires that the initialized sample $\mathbf{x}_0\in\mathbb{R}^d$ will have a reasonably-large projection onto the $(d-k)$-dimensional subspace orthogonal to vectors $\boldsymbol{\mu}_1,\\ldots , \boldsymbol{\mu}_k$. By only assuming that $d-k$ is moderately large, this assumption is guaranteed to hold according to Proposition 2. This is following the intuition that  a $d$-dimensional Gaussian vector $\mathbf{x}_0 \sim \mathcal{N}(\mathbf{0}_d,\mathbf{I}_d)$ will (with high probability) possess a $\mathcal{O}(\sqrt{d-k})$-magnitude projection onto a fixed $(d-k)$-dimensional space. 


Regarding Assumption 1, we note that as discussed in rebuttal #1, $P^{(0)}$ is supposed to be a large-variance (yet low-probability) mode *surrounding the small-variance and high-probability modes $P^{(1)}, \cdots, P^{(k)}$*. Assumption 1 formalizes this intuition as the center of the high-variance mode $P^{(0)}$ should not be exceedingly far from the other low-variance modes, otherwise $P^{(0)}$ will not be capable of surrounding the extremely far modes and dominating the Langevin dynamics. Please note that as described in the second paragraph of rebuttal #1, the theoretical result can be further extended to a hardness result assuming the large-enough distance between the low-variance modes $P^{(1)}, \cdots, P^{(k)}$, which requires sufficiently large distance between the means of $P^{(1)}, \cdots, P^{(k)}$ (consistent with the reviewer’s intuition).
We thank Reviewer 3VKE for his/her feedback on our response. We are pleased to hear the reviewer finds the rebuttal satisfactory. Regarding the raised point, we think the reviewer means whether Theorem 1 holds for any $T \le \exp(\mathcal{O}(d))$ ($t$ replaced by $d$). If so, we confirm the reviewer's interpretation of the statement, as the bounds in Theorems 1-4 hold for any $T \le \exp(\mathcal{O}(d))$. We will clarify this point in the revised text.
Thanks for your thorough response, especially the clarification for the TV lower bound as well as the intuition for your setting.

One final request:

- I assume that in Theorem 1, the $T = \exp(O(t))$ condition means the result holds for _any_ $T$ such that $T \leqslant \exp(O(t))$.

I have updated my score to take into account the respones.
Thank you authors for these clarifications: they are quite helpful and would be useful to add to the paper. 

As it is currently written, the paper suggests that:
- **Theorem 1: limitation**. Langevin on the full target can fail to find some modes.
- **Theorem 2: limitation**. Annealed Langevin on the full target can fail to find some modes.
- **Theorem 5: solution**. Annealed Langevin on the target conditionals may solve the problem. This is suggested in writing, given that this method is proposed "to reduce the mode-seeking tendencies of vanilla and annealed Langevin dynamics".

**Comments on Theorems 1 and 2**. As you have said in your comment above, Theorems 1 and 2 are more subtle than this. They point out, on the particle-level, that some target modes may be missed with high probability if the initial distribution is not "noisy enough". This is coherent with the observation of Song and Ermon (2020) who related "noisy enough" to the distance between data points. This result at the particle-level suggests that, at the distribution-level, convergence may be hindered. But that is not automatic. To formally claim that "mode-missing" by particles results in "slower convergence of the probability law", we will need the lower bounds in TV that you stated. 

Do you have these lower bounds for **both** Theorems 1 and 2? So for both Vanilla and Annealed Langevin?

**Comment on Theorem 5**. As it is stated, it is not at all clear that Theorem 5 "reduce[s] the mode-seeking tendencies of vanilla and annealed Langevin dynamics". Theorem 5 is formally stated as: if you sample the conditionals with a given precision, then you sample the joint with a given precision. But the difficulty of sampling the conditionals is not clear. Again, this theorem is stated bluntly in the text and is not commented. Your argument above is useful: comparing $R^2$ of order $d$ vs. $Q$ and explaining why Langevin on the target conditionals has a nicer upper-bound than Annealed Langevin on the full target. It should be part of the explanation.
I thank the authors for their feedback.
Actually, the fact that similar results would hold if $P_0 \sim \mathcal{N}(0, 1)$ under the additional assumption that the initial sample will be far from the vector space of $\left\{ \mathbf{\mu}_i \right\}_{i \in [k]}$ worries me a bit about the significance of the work as such assumption cannot be expected to hold in practice. After having a closer look to Assumption 1, the lower bound feels somewhat odd as we could expect it to depend on the distance between the modes and not assume a priori that they are close.
We thank Reviewer nEdB for his/her time and feedback on our response. Regarding the raised points:

**1- Theoretical results on annealed Langevin dynamics** 

**Re:** Thank you for letting us know about the unclarity in interpreting Theorem 2 on annealed Langevin dynamics (ALD). We would like to clarify that Theorem 2 does not aim to highlight a serious limitation of the ALD approach. Please note that our analysis focuses on the diversity of generated samples under a multi-modal distribution. Especially, Theorem 2 in our work suggests that if we consider an upper-bound $c_\sigma$ on the noise levels $\sigma = O(1)$ in ALD that remains a constant in dimension $d$, then the generated samples over a sub-exponential iterations (in dimension $d$) could miss low-variance modes separated from the initialization mode $P^{(0)}$. Of course, we assume a constant bound (dimension-independent) on the noise level $c_\sigma=O(1)$ in ALD, which does not imply the ALD method would suffer from mode dropping in the general case. 

We note that the proper selection of noise level in ALD has been acknowledged in the literature. For example, Song and Ermon (2020) [2] explain in their paper: "it is necessary for $\sigma_0$ to be numerically comparable to the maximum pairwise distances of data to facilitate transitioning of Langevin dynamics and hence improving sample diversity" (on page 4 of [2]). 

In the revision, we will be clear about the constant noise level assumption in our writing, and use the term “annealed Langevin dynamics with *bounded noise level*” to ensure the result will not be misinterpreted as a general limitation of ALD.

**2- Iteration complexity of Chained Langevin dynamics** 

**Re:** On a high level, Langevin dynamics performs a *noisy local search* to generate samples around the peaks of the likelihood function. When sampling from $P(\mathbf{x})$ for a high-dimensional $\mathbf{x}\in\mathbb{R}^d$, the algorithm has to randomly explore a large volume growing exponentially in $d$ to find the high-probability yet low-variance modes $P^{(1)},\ldots , P^{(k)}$. On the other hand, when sampling from the conditional distribution $P(\mathbf{x}^{(q)} \mid \mathbf{x}^{(1)}, \cdots, \mathbf{x}^{(q-1)})$, the algorithms needs to only search over a $Q$-dimensional space ($Q$ being the patch size where $Q\ll d$) to find the peaks of the resulting multi-modal conditional density. Therefore, one can expect a faster convergence to the support set of target modes $P^{(1)},\ldots , P^{(k)}$.

Regarding the reviewer’s question on applying Theorem 1 in [1], we want to clarify that this theorem states
$$
\tau(\varepsilon) \le \mathcal{O} \left( \exp(32LR^2) \kappa^2 \frac{d}{\varepsilon^2} \ln \left( \frac{d}{\varepsilon^2} \right) \right). 
$$

The above bound is exponential in $R^2$ ($R$ is the radius of strong convexity of $\log(P)$ where $P$ is the density function), which in the Gaussian mixture case of our theorems will scale linearly with dimension $d$. For example, if we directly apply the theorem to the joint target of the Gaussian mixture case in our synthetic experiments (Section 6), we have $R^2 \ge d$, which means Langevin dynamics is expected to require $\mathcal{O}(\exp(32Ld))$ iterations (Theorem 1 in [1]). On the other hand, chained Langevin dynamics breaks the sample into patches of $Q$-dimension, for which case $R^2 = Q$ (for constant $Q$) and thus making the term $\mathcal{O}(\exp(32LQ))$ to be independent of $d$. 



[1] Ma, Y. A., Chen, Y., Jin, C., Flammarion, N., & Jordan, M. I. (2019). Sampling can be faster than optimization. Proceedings of the National Academy of Sciences, 116(42), 20881-20885.

[2] Song, Y. and Ermon, S. (2020). Improved techniques for training score-based generative models. Advances in neural information processing systems, 33:12438–12448.
I thank the authors for their response. However, two points are still quite unclear to me. 

**Explaining the limitation of Annealed Langevin dynamics**. I thank the authors for their clarification on the bounds for Vanilla and Annealed Langevin Dynamics. I understand these bounds have hidden constants and therefore one cannot easily claim that Annealed Langevin Dynamics performs worse. That being said, the authors seem to claim a serious limitation of Annealed Langevin Dynamics in Theorem 2, which I interpret to be that particles following this process will always be "far enough" from some mode and this hinders convergence. As I stated above: this is surprising, given the widespread success of Annealed Langevin Dynamics for efficiently sampling from multimodal target distributions. Reviewer Ave7 also picked up on this, saying "unlike what was believed, annealed Langevin also fails." Because this is surprising, Theorem 2 deserves some commentary, some high-level argument explaining why this popular method may fail. Pointing to the proof is not enough. As it stands, Theorem 2 is formally stated at the end of section 4.1 with no explanation. 

**Explaining the benefits of Chained Langevin dynamics**. While I understand the authors' result, I still find it unclear why using the conditionals of the target $p(x_i | x_{-i})$ makes the sampling problem easier than sampling from the joint target $p(x)$. The authors mention Theorem 1 of [1] which they apply to each conditional $p(x_i | x_{-i})$, but couldn't this Theorem be applied directly to the joint target directly $p(x)$? The authors also mention that the joint target has a multi-dimensional input, while the conditionals have a one-dimensional input: still, that is not a satisfying explanation. Sampling from a one-dimensional distribution is not necessarily easier. Again, could the authors provide some high-level arguments on why sampling from the conditional distributions helps, versus sampling from the joint distribution? 

[1] Ma, Y. A., Chen, Y., Jin, C., Flammarion, N., & Jordan, M. I. (2019). Sampling can be faster than optimization. Proceedings of the National Academy of Sciences, 116(42), 20881-20885.