Score matching is a method introduced by Hyvärinen (2005) for learning the parameters of a distribution from data, useful for parametric families in which evaluating the likelihood is intractable. An illustrative example are energy-based models (EBMs), parametric families of the form p θ (x) ∝ exp(E θ (x)), for which evaluating and optimizing the likelihood is comutationally hard due to the partition function Z θ = x exp(E θ (x)) (Pabbaraju et al., 2023). Score matching obviates evaluating the partition function by instead fitting the score of the distribution (i.e., ∇ x log p(x)). While there is algorithmic gain, the statistical cost can be substantial. In recent work, Koehler et al. (2022) show that score matching is statistically much less efficient (i.e., the estimation error, given the same number of samples is much bigger) than maximum likelihood when the distribution being estimated has poor isoperimetric properties (i.e. a large Poincaré constant). However, even very simple multimodal distributions like a mixture of two Gaussians with far away meansÐhave a very large Poincaré constant. As many distributions of interest (e.g., images) are multimodal in nature, the score matching estimator is likely to be statistically untenable.
In the generative models literature, the seminal paper by Song and Ermon (2019) proposes a way to deal with multimodality and manifold structure in the data by annealing: namely, estimating the scores of convolutions of the data distribution with different levels of Gaussian noise. The intuitive explanation they propose is that the distribution smoothed with more Gaussian noise is easier to estimate (as there are no parts of the distribution that have low coverage by the training data), which should help estimate the score at lower levels of Gaussian noise. However, making this quantitative or formal seems very challenging.
In this paper, we show that there is a deep connection between the mixing time of broad classes of continuous, time-homogeneous Markov processes with stationary distribution p and generator L, and the statistical efficiency of an appropriately chosen generalized score matching loss (Lyu, 2012) that tries to match Op p . In the case that L is the generator of Langevin diffusion, and O = ∇ x , we recover the results of Koehler et al. (2022). This ªdictionaryº allows us to design score losses with better statistical behavior, by ªtranslatingº techniques for speeding up Markov chain convergence Ð e.g. preconditioning a diffusion and lifting the chain by introducing additional variables.
Our contributions are as follows:
1. A general framework for designing generalized score matching losses with good asymptotic sample complexity from fast-mixing diffusions. Precisely, for a broad class of diffusions with generator L and Poincaré constant C P , we can choose a linear operator O, such that the generalized score matching loss
1 2 E p Op p -Op θ p θ 2 2
has statistical complexity that is a factor C 2 P worse than that of maximum likelihood. (Recall, C P characterizes the mixing time of the Markov process with generator L in chi-squared distance.) In particular, for diffusions that look like ªpreconditionedº Langevin, this results in ªappropriately preconditionedº score loss.
2. We analyze a lifted diffusion, which introduces a new variable for temperature, and we provably show statistical benefits of annealing for score matching. Precisely, we exhibit continuouslytempered Langevin (CTLD), a Markov process which mixes in time poly(D, d, 1/λ min , λ max ) for finite mixtures of Gaussians in ambient dimension d with identical covariances whose smallest and largest eigenvalues are lower and upper bounded by λ min and λ max respectively, and means lying in a ball of radius D. (Note, the bound has no dependence on the number of components.) Moreover, the corresponding generalized score matching loss is a form of the annealed score matching loss introduced in (Song and Ermon, 2019;Song et al., 2020), with a particular choice of weighing for the different amounts of Gaussian convolution. This is the first result formally showing the statistical benefits of annealing for score matching. Technically, this result involves bounding the mixing time of a ªcontinuously temperedº version of Langevin dynamics for mixtures, using functional decomposition theorems.
On a conceptual level, our work draws on and brings together theoretical developments in understanding score matching, as well as designing and analyzing faster-mixing Markov chains based on strategies in annealing. An in-depth review of related prior work is included in Appendix I.
The conventional score-matching objective (Hyvärinen, 2005) is defined as
Note, the expression is asymmetric: p is the data distribution, p θ is the distribution that is being fit.
Written like this, it is not clear how to minimize this loss, when we only have access to data samples from p. The main observation of Hyvärinen (2005) is that the objective can be rewritten (under suitable decay conditions using integration by parts) in a form that is easy to fit given samples:
where K p is some constant independent of q. To turn this into an algorithm given samples, one simply solves
where p denotes the uniform distribution over the samples from p. This objective can be calculated efficiently given samples from p, so long as the gradient and Hessian of the log-pdf of p θ can be efficiently calculated. 1 Generalized Score Matching, first introduced in Lyu (2012), generalizes ∇ x to an arbitrary linear operator O:
Definition 1 The Generalized Score Matching (GSM) loss with a linear operator O acting on density functions is defined as
Again, for the special case of the family of operators O in (3), the integration by parts form of the objective can be easily written down explicitly (the proof is provided in Appendix B):
Lemma 2 (Integration by parts for the GSM in (3)) The generalized score matching objective in (3) satisfies the equality
In this section, we introduce the key definitions related to continuous-time Markov chains and diffusion processes:
Definition 2 (Markov semigroup) We say that a family of functions {P t (x, y)} t≥0 on a state space Ω is a Markov semigroup if P t (x, •) is a distribution on Ω and P t+s (x, dy) = Ω P t (x, dz)P s (z, dy) for all x, y ∈ Ω and s, t ≥ 0.
Definition 3 (Time-homogeneous Markov processes) A time-homogeneous Markov process (X t ) t≥0 on state space Ω is defined by a Markov semigroup {P t (x, y)} t≥0 as follows: for any measurable A ⊆ Ω, Pr(X s+t ∈ A|X s = x) = P t (x, A) = A P t (x, dy). Moreover, P t can be thought of as acting on a function g as
A particularly important class of time-homogeneous Markov processes is given by Itô diffusions, namely stochastic differential equations of the form dX t = b(X t )dt + σ(X t )dB t for a drift function b, and a diffusion coefficient function. In fact, a classical result due to Dynkin (Rogers and Williams (2000), Theorem 13.3) states that any ªsufficiently regularº time-homogeneous Markov process (specifically, a process whose semigroup is Feller-Dynkin) can be written in the above form.
We will be interested in Itô diffusions whose stationary distribution is a given distribution p(x) ∝ exp(-f (x)). Perhaps the most well-known example of such a diffusion is Langevin diffusion, namely dX t = -∇f (X t )dt + √ 2dB t , where B t is standard Brownian motion in R d . In fact, a completeness result due to Ma et al. (2015) states that we can characterize all Itô diffusions whose stationary distribution is p(x) ∝ exp(-f (x)) as ªpreconditionedº Langevin diffusion. Precisely:
Theorem 3 (Characterization of Itô diffusions w/ given stationary distribution, Ma et al. (2015)) Any Itô diffusion with stationary distribution p(x) ∝ exp(-f (x)) can be written in the form:
Intuitively, D(x) can be viewed as ªreshapingº the diffusion, whereas Q and Γ are ªcorrection termsº to the drift so that the stationary distribution is preserved.
Definition 4 The generator L corresponding to Markov semigroup is Lg = lim t→0 Ptg-g t . Moreover, if p is the unique stationary distribution, the Dirichlet form and the variance are
respectively. We will use the shorthand E(g) := E(g, g). By Itô's Lemma, the generator of diffusions of the form (5) have the form:
The Dirichlet form for diffusions of the form (5) also has a very convenient form:
Lemma 4 (Dirichlet form of continuous Markov Process) An Itô diffusion of the form (5) has a Dirichlet form
Notably, for Langevin diffusion, the Dirichlet form is just the l 2 norm of ∇g:
For a general diffusion of the form (5), we can think of D(x) as a (point-specific) preconditioner, specifying the norm with respect to which to measure ∇g. The proof of this lemma is given in Appendix B. Finally, we define the Poincaré constant:
Definition 5 (Poincaré inequality) A continuous-time Markov process satisfies a Poincaré inequality with constant C if for all functions g such that E(g) is defined (finite), we havefoot_0 E(g) ≥ 1 C Var p (g). We will abuse notation, and for a Markov process with stationary distribution p, denote by C P the Poincaré constant of p, the smallest C such that above Poincaré inequality is satisfied.
The Poincaré inequality implies exponential ergodicity for the χ 2 -divergence, namely χ 2 (p t , p) ≤ e -2t/C P χ 2 (p 0 , p), where p is the stationary distribution of the chain and p t is the distribution after running the Markov process for time t, starting at p 0 .
We will analyze mixing times using a decomposition technique similar to the ones employed in Ge et al. (2018); Moitra and Risteski (2020). Intuitively, these results ªdecomposeº the Markov chain by partitioning the state space into sets, such that: (1) the mixing time of the Markov chain inside the sets is good; (2) the ªprojectedº chain, which transitions between sets with probability equal to the probability flow between sets, also mixes fast. An example of such a result is Theorem 6.1 from Ge et al. (2018):
Theorem 5 (Decomposition of Markov Chains, Theorem 6.1 in Ge et al. (2018)) Let M = (Ω, L) be a continuous-time Markov chain with stationary distribution p and Dirichlet form E(g, g) = -⟨g, Lg⟩ p . Suppose the following hold.
1. The Dirichlet form for L decomposes as ⟨f, Lg⟩ p = m j=1 w j ⟨f, L j g⟩ p j , where p = m j=1 w j p j and L j is the generator for some Markov chain M j on Ω with stationary distribution p j .
2. (Mixing for each M j ) The Dirichlet form E j (f, g) = -⟨f, Lg⟩ p j satisfies the Poincaré inequality Var p j (g) ≤ CE j (g, g). [ḡ(k)ḡ(j)] P (j, k), where P (j, k) = w k max{χ 2 (p j , p k ), χ 2 (p k , p j ), 1} .
Let p be the stationary distribution of M . Suppose M satisfies the Poincaré inequality Var p(ḡ) ≤ C Ē(g, g).
Then M satisfies the Poincaré inequality:
We will also use some classical results about asymptotic normality of M-estimators, under standard identifiability and differentiability conditions. Namely, we will be considering estimators defined as
In general, we will denote by n the number of samples, and Ê will denote an empirical average, that is the expectation over the n training samples. Finally, θn will denote min θ∈Θ Ê[ℓ θ (x)] when the number of samples is n. In Section A.2 we recall sufficient conditions for θn to be asymptotically normal, and an expression for the asymptotic covariance.
The goal of this section is to provide a general framework that provides a bound on the sample complexity of a generalized score matching objective with operator O, under the assumption that some Markov process with generator L mixes fast. Precisely, we will show:
Theorem 6 (Main, sample complexity bound) Consider an Itô diffusion of the form (5) with stationary distribution p(x) ∝ exp(-f (x)) and Poincaré constant C P with respect to the generator of the Itô diffusion. Consider the generalized score matching loss with operator (Og)(x) :=
. Suppose we are optimizing this loss over a parametric family {p θ : θ ∈ Θ} that satisfies:
1. (Asymptotic normality) Let Θ * be the set of global minima of the generalized score matching loss D GSM , that is
Suppose the generalized score matching loss is asymptotically normal: namely, for every θ * ∈ Θ * , and every sufficiently small neighborhood S of θ * , there exists a sufficiently large n, such that there is a unique minimizer θn of Êl θ (x) in S, where:
Furthermore, assume θn satisfies
2. (Realizibility) At any θ * ∈ Θ * , we have p θ * = p.
Then, we have:
Remark 7 The two terms on the right hand sides qualitatively capture two intuitive properties necessary for a good sample complexity: the factor involving the covariances can be thought of as a smoothness term capturing how regular the score is as we change the parameters in the family we are fitting; the C P term captures how the error compounds as we ªextrapolateº the score into a probability density function.
Remark 8 This theorem generalizes Theorem 2 in Koehler et al. (2022), who show the above only in the case of L being the generator of Langevin diffusion and O = ∇ x , i.e. when D GSM is the standard score matching loss. Furthermore, they only consider the case of p θ being an exponential family, i.e. p θ (x) ∝ exp(⟨θ, T (x)⟩) for some sufficient statistics T (x). Finally, just as in Koehler et al. (2022), we can get a tighter bound by replacing C P by the restricted Poincaré constant, which is the Poincaré constant when considering only the functions of the form ⟨w, ∇ θ log p θ (x) | θ=θ * ⟩. Remark 9 Note that if we know √ n( θnθ * ) d -→ N (0, Γ SM ), we can extract bounds on the expected ℓ 2 2 distance between θn and θ * . Namely, from Markov's inequality (see e.g., Remark 4 in Koehler et al. ( 2022)), we have for sufficiently large n, with probability at least 0.99 it holds that
Some conditions for asymptotic normality can be readily obtained by applying standard results from asymptotic statistics (e.g. Van der Vaart (2000), Theorem 5.23, reiterated as Lemma 19 for completeness).From that lemma, when an estimator θ = arg min Êl θ (x) is asymptotically normal, we
, where L(θ) = E θ l(x). Therefore, to bound the spectral norm of Γ SM , we need to bound the Hessian and covariance terms in the expression above. The latter is a fairly straightforward calculation, which results in the following Lemma, proven in Appendix C.
Lemma 10 (Bound on smoothness) For l θ (x) defined in Theorem 6,
The bound on the Hessian is where the connection to the Poincaré constant manifests. Namely, we show:
Lemma 11 (Bounding Hessian) Let the operators O, L be such that for every vector w, the function g
Proof To reduce notational clutter, we will drop | θ=θ * since all the functions of θ are evaluated at θ * . Consider an arbitrary direction w. We have:
1 follows from a straightforward calculation (in Lemma 5), 2 follows from the definition of Poincaré inequality of a diffusion process with Dirichlet form derived in Lemma 4, applied to the function ⟨w, ∇ θ log p θ ⟩, and 3 follows since
Fisher matrix ( Van der Vaart, 2000)). Since this holds for every vector w, we have
. By monotonicity of the matrix inverse operator (Toda, 2011), the claim of the lemma follows.
In this section, we instantiate the framework from the previous section to the specific case of a Markov process, Continuously Tempered Langevin Dynamics, which is a close relative of simulated tempering (Marinari and Parisi, 1992), where the number of ªtemperaturesº is infinite, and we temper by convolving with Gaussian noise. We show that the generalized score matching loss corresponding to this Markov process mixes in time poly(D, d) for a mixture of K Gaussians (with identical covariance) in d dimensions, and means in a ball of radius D. More precisely, in this section, we will consider the following family of distributions: Assumption 1 Let p 0 be a d-dimensional Gaussian distribution with mean 0 and covariance Σ. We will assume the data distribution p is a K-Gaussian mixture, namely p = K i=1 w i p i , where p i (x) = p 0 (xµ i ), i.e. a shift of the distribution p 0 so its mean is µ i . We will assume the means µ i lie within a ball with diameter D. We will denote the min and max eigenvalues of covariance with λ min (Σ) = λ min and λ max (Σ) = λ max . We will denote the min and max mixture proportion with min i w i = w min and max i w i = w max . Let Σ β = Σ + βλ min I d be the shorthand notation of the covariance of individual Gaussian at temperature β.
Mixtures of Gaussians are one of the most classical distributions in statisticsÐand they have very rich modeling properties. For instance, they are universal approximators in the sense that any distribution can be approximated (to any desired accuracy), if we consider a mixture with sufficiently many components (Alspach and Sorenson, 1972). A mixture of K Gaussians is also the prototypical example of a distribution with K modes Ð the shape of which is determined by the covariance of the components.
Note at this point we are just saying that the data distribution p can be described as a mixture of Gaussians, we are not saying anything about the parametric family we are fitting when optimizing the score matching lossÐwe need not necessarily fit the natural unknown parameters (the means, covariances and weights). The primary reason this family of distributions is convenient for technical analysis is a closure property under convolutions: a convolution of a Gaussian mixture with a Gaussian produces another Gaussian mixture. Namely, the distributivity property of the convolution operator implies:
Proposition 1 (Convolution with Gaussian) Under Assumption 1, the distribution p * N (x; 0, σfoot_1 I) satisfies p * N (x; 0, σ 2 I) = i w i p 0 (xµ i ) * N (x; 0, σ 2 I) and (p 0 (xµ i ) * N (x; 0, σ 2 I)) is a multivariate Gaussian with mean µ i and covariance Σ + σ 2 I.
The Markov process we will be analyzing (and the corresponding score matching loss) is a continuous-time analog of the Simulated Tempering Langevin chain introduced in Ge et al. ( 2018):
Definition 6 (Continuously Tempered Langevin Dynamics (CTLD)) We will consider an SDE over a temperature-augmented state space, that is a random variable 1+β) and
λ min -1. Let p β := p * N (0, βλ min I d ) denotes the distribution p convolved with a Gaussian of covariance βλ min I d . Furthermore, L(dt) is a measure supported on the boundary of the interval [0, β max ] and ν t is the unit normal at the endpoints of the interval, such that the stationary distribution of this SDE is p(x, β) = r(β)p β (x) (Saisho, 1987) foot_2 .
If we ignore the boundary reflection term, the updates for CTLD are simply Langevin dynamics applied to the distribution p(x, β), and r(β) specifies the distribution over the different levels of noise. CTLD can be readily seen as a ªcontinuous-timeº analogue of the usual simulated tempering chain. Namely, in the usual (discrete-time) simulated tempering (Lee et al., 2018;Ge et al., 2018), the tempering chain has two types of moves: one which evolves the position in the current temperature, and one which tries to change the temperatures, followed by a Metropolis Hastings filtering step.
We point out several similarities and crucial differences with the chain proposed in Ge et al. (2018). The chain in Ge et al. (2018) has a finite number of temperatures and the distribution in each temperature is defined as scaling the log-pdf, rather than convolution with a GaussianÐthis is because the mode of access in Ge et al. (2018) is the gradient of the log-pdf, whereas in score matching, we have samples from the distribution. The distributions in Ge et al. (2018) are geometrically spaced outÐso β being distributed as exp(-Θ(β)) in our case can be thought of as a natural continuous analogue.
Since CTLD amounts to performing (reflected) Langevin dynamics on the appropriate joint distribution p(x, β), the corresponding generator L for CTLD is also readily written down (Proposition 6). The operator O that corresponds to the CTLD is also easy to derive:
Proof The first equality follows as a special case of Langevin on the lifted distribution. The second equality follows by writing ∇ β log p(x|β) and ∇ β log p θ (x|β) through the Fokker-Planck equation for p(x|β) (see Lemma 28).
This loss was derived from first principles from the Markov Chain-based framework in Section 3, however, it is readily seen that this loss is a ªsecond-orderº version of the annealed losses in Song and Ermon (2019); Song et al. (2020) Ð the weights being given by the distribution r(β).
With this setup in mind, we can proceed to the main technical results of this section.
Theorem 12 (Poincaré constant of CTLD) Under Assumption 1, the Poincaré constant of CTLD C P enjoys the upper bound
Note that perhaps surprisingly, the above result has no dependence on the number of components, or on the smallest component weight w min Ðonly on the diameter D, the ambient dimension d, and λ min and λ max . The results in Ge et al. (2018) have a dependence on w min , but in their model, it's not possible to query convolutions of the pdf with a Gaussian. This result can be seen as a ªtime-homogenousº analogue of recent results (Lee et al., 2023) that the reverse SDE (which is timeinhomogenous) converges to the data distribution in polynomial time under minimal assumptions (Lipschitzness of the score). This result is of independent technical interest as it illustrates the power of having an oracle for convolutions of the target distribution.
To get a complete bound on the asymptotic sample complexity of generalized score matching, according to the framework from Lemma 3, we also need to bound the smoothness terms as in Lemma 10. These terms of course depend on the choice of parametrization for the family of distributions we are fitting Ð in particular, there is no ªcanonicalº parametrization for multimodal distributions. To get a quantitative sense for how the smoothness might scale, we will consider one natural parametrization for a mixture:
Assumption 2 Consider the case of learning unknown meansfoot_3 , such that the parameters to be learned are a vector θ = (µ 1 , µ 2 , . . . , µ K ) ∈ R dK .
With this parametrization, the smoothness term can be bounded as follows:
Theorem 13 (Smoothness under the natural parameterization) Under Assumptions 1 and 2, the smoothness defined in Theorem 6 enjoys the upper bound
Finally, we show that the generalized score matching loss is asymptotically normal. The proof of this is in Appendix F, and proceeds by verifying the conditions of Lemma 19. Putting this together with the Poincaré inequality bound Theorem 12 and Theorem 6, we get a complete bound on the sample complexity of the generalized score matching loss with O: Theorem 14 (Main, Polynomial Sample Complexity Bound of CTLD)
Let the data distribution p satisfy Assumption 1. Then, the generalized score matching loss defined in Proposition 7 with parametrization as in Assumption 2 satisfies:
1. The set of optima
satisfies:
2. Let θ * ∈ Θ * and let C be any compact set containing θ * . Denote
Finally, let D be any closed subset of C not intersecting C 0 . Then, we have:
3. For every θ * ∈ Θ * and every sufficiently small neighborhood S of θ * , there exists a sufficiently large n, such that there is a unique minimizer θn of Êl θ (x) in S. Furthermore, θn satisfies:
We provide some brief comments the theorem:
• The goal of this result is not to provide a new algorithm for learning mixtures of Gaussians, but to provide a paradigm case in which annealing improves the statistical behavior of score matching. In fact, the generalized score matching loss is not necessarily convex, and our result is only statistical in nature. Moreover, from the point of view of score matching, the statistical complexity relative to MLE worsens when the modes are well-separated Ð which is exactly the case when many algorithms for learning mixtures of Gaussians tend to work. (Appendix I).
• Condition (1) is the standard identifiability condition (Yakowitz and Spragins, 1968) for mixtures of Gaussians: the means are identifiable up to ªrenamingº the components. This is inevitable if some of the weights are equal; if all the weights are distinct, Θ * would in fact only consist of one point, s.t. ∀i ∈ [K], µ i = µ * i . Condition (2) says that asymptotically, the empirical minimizers of D GSM are the points in Θ * . It can be viewed as (and follows from) a uniform law of large numbers.
• Condition (3) characterizes the sample complexity of minimizers in the neighborhood of each of the points in Θ * , and is a consequence of the CTLD Poincaré inequality estimate (Theorem 12) and the smoothness estimate (Theorem 13). Note that in fact the RHS of point 3 has no dependence on the number of components. This makes the result extremely general: the loss compared to MLE is very mild even for distributions with a large number of modes.foot_4
In this section, we will sketch the proof of Theorem 12. By slight abuse of notation, we will define the distribution of the ªindividual componentsº of the mixture at a particular temperature, namely for
Correspondingly, we will denote the conditional distribution for the i-th component by p(x, β|i) ∝ r(β)N (x; µ i , Σ + βλ min I d ).
The proof will proceed by applying the decomposition Theorem 5 to CTLD. Towards that, we denote by E i the Dirichlet form corresponding to Langevin with stationary distribution p(x, β|i). By Propositions 6, it's easy to see that the generator for CTLD satisfies E = i w i E i . This verifies condition (1) in Theorem 5. To verify condition (2), we will show Langevin for each of the distributions p(x, β|i) mixes fast (i.e. the Poincaré constant is bounded). To verify condition (3), we will show the projected chain ªbetweenº the components (as defined in Theorem 5) mixes fast. We will expand on each of these parts in turn. Fast mixing within a component: The first claim we will show is that we have fast mixing ªinsideº each of the components of the mixture. Formally, we show:
The proof of this lemma proceeds via another (continuous) decomposition theorem. Intuitively, what we show is that for every β, p(x|β, i) has a good Poincaré constant; moreover, the marginal distribution of β, which is r(β), is log-concave and supported over a convex set (an interval), so has a good Poincaré constant. Putting these two facts together via a continuous decomposition theorem (Theorem D.3 in Ge et al. ( 2018)), we get the claim of the lemma. The details are in Appendix E.1.
Fast mixing between components: Next, we show the ªprojectedº chain between the components mixes fast (full proofs are in Appendix E.2):
Lemma 16 (Poincaré constant of projected chain) Define the projected chain M over [K] with transition probability
where χ 2 max (p, q) = max{χ 2 (p, q), χ 2 (q, p)}. If j̸ =i T (i, j) < 1, the remaining mass is assigned to the self-loop T (i, i). The stationary distribution p of this chain satisfies p(i) = w i . Furthermore, the projected chain has Poincaré constant C ≲ D 2 λ -1 min .
The intuition for this claim is that the transition probability graph is complete, i.e. T (i, j) ̸ = 0 for every pair i, j ∈ [K]. Moreover, the transition probabilities are lower bounded, since the χ 2 distances between any pair of ªannealedº distributions p(x, β|i) and p(x, β|j) can be upper bounded. The reason for this is that at large β, the Gaussians with mean µ i and µ j are smoothed enough so that they have substantial overlap; moreover, the distribution r(β) is set up so that enough mass is placed on the large β. The precise lemma bounding the χ 2 divergence between the components is:
Lemma 17 For every i, j ∈ [K], we have χ 2 (p(x, β|i), p(x, β|j)) ≤ 14D 2 λ -1 min .
To obtain Theorem 13, we note ∥Cov (O∇ θ log p θ ) ∥ OP and ∥Cov ((O + O)∇ θ log p θ ) ∥ OP can be completely characterized by bounds on the higher-order derivatives with respect to x and µ i of the log-pdf, since derivatives with respect to β can be related to derivatives with respect to x via the Fokker-Planck equation (Lemma 28). The main technical tools involved are: (1) the convexity of the perspective map to relate derivatives of the mixture to derivatives of the components (Lemma 37);
(2) bounds on derivatives of the components via Hermite polynomial machinery (Lemma 38); (3) bounds on logarithmic derivatives via higher-order versions of the Faá di Bruno formula (Constantine and Savits, 1996). The complete proofs are in Appendix G.
In this paper, we provide a general framework about designing statistically efficient generalized score matching losses from fast-mixing Markov Chains. As a demonstration of the power of the framework, we provide the first formal analysis of the statistical benefits of annealing for score matching for multimodal distributions. A core technical result of this part is bounding the mixing time for a continuously tempered version of Langevin diffusion. The framework can be likely used to analyze other common continuous and discrete Markov Chains (and corresponding generalized score losses), like underdamped Langevin dynamics and Gibbs samplers. Contents 1 Introduction 2 Preliminaries 2.1 Generalized Score Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Continous-time Markov Processes . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Dirichlet forms and Poincaré inequalities . . . . . . . . . . . . . . . . . . . . . . . 2.4 Asymptotic efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Main Results: A Framework for Analyzing Generalized Score Matching 4 Main Results: Benefits of Annealing 4.1 Proof Sketch: Bounding the Poincaré Constant of CTLD . . . . . . . . . . . . . . 4.2 Proof Sketch: Bounding the Smoothness Terms . . . . . . . . . . . . . . . . . . . 5 Conclusion A Preliminaries A.1 Continuous Markov Chain Decomposition . . . . . . . . . . . . . . . . . . . . . . A.2 Asymptotic normality of M-estimators . . . . . . . . . . . . . . . . . . . . . . . . A.3 Hermite Polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 Logarithmic derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Moments of mixtures and the perspective map . . . . . . . . . . . . . . . . . . . . B Derivations of generators and score losses for diffusions C A Framework for Analyzing Generalized Score Matching D Overview of Continuously Tempered Langevin Dynamics E Polynomial mixing time bound: proof of Theorem 12 E.1 Mixing inside components: proof of Lemma 15 . . . . . . . . . . . . . . . . . . . E.2 Mixing between components: proof of Lemma 16 . . . . . . . . . . . . . . . . . F Asymptotic normality of generalized score matching for CTLD G Polynomial smoothness bound: proof of Theorem 13 G.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G.2 Detailed proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . H Technical Lemmas H.1 Moments of a chi-squared random variable . . . . . . . . . . . . . . . . . . . . . I Related work A.3. Hermite Polynomials
To obtain polynomial bounds on the moments of derivatives of Gaussians, we will use the known results on multivariate Hermite polynomials.
Definition 7 (Hermite polynomial, (Holmquist, 1996)) The multivariate Hermite polynomial of order k corresponding to a Gaussian with mean 0 and covariance Σ is given by the Rodrigues formula:
where ϕ(x; Σ) is the pdf of a d-variate Gaussian with mean 0 and covariance Σ, and ⊗ denotes the Kronecker product.
Note that ∇ ⊗k x can be viewed as a formal Kronecker product, so that ∇ ⊗k x f (x), where f : R d → R is a C k -smooth function gives a d k -dimensional vector consisting of all partial derivatives of f of order up to k.
Proposition 3 (Integral representation of Hermite polynomial, (Holmquist, 1996)) The Hermite polynomial H k defined in Definition 7 satisfies the integral formula:
where ϕ(x; Σ) is the pdf of a d-variate Gaussian with mean 0 and covariance Σ.
Note, the Hermite polynomials are either even functions or odd functions, depending on whether k is even or odd:
This property can be observed from the Rodrigues formula, the fact that ϕ(•; Σ) is symmetric around 0, and the fact that ∇ -x = -∇ x .
We establish the following relationship between Hermite polynomial and (potentially mixed) derivatives in x and µ, which we will use to bound several smoothness terms appearing in Section G.
Lemma 20 If ϕ(x; Σ) is the pdf of a d-variate Gaussian with mean 0 and covariance Σ, we have:
where the left-hand-side is understood to be shaped as a vector of dimension R d k 1 +k 2 .
Proof Using the fact that ∇ x-µ = ∇ x in Definition 7, we get:
Since the Kronecker product satisfies the property (A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD), we have (Σ∇ x ) ⊗k = Σ ⊗k ∇ ⊗k x . Thus, we have:
Since ϕ(µx; Σ) is symmetric in µ and x, taking derivatives with respect to µ we get:
Rearranging again and using (7), we get:
Combining ( 8) and ( 9), we get:
Applying the integral formula from Proposition 3, we have:
as we needed.
Now we are ready to obtain an explicit polynomial bound for the mixed derivatives for a multivariate Gaussian with mean µ and covariance Σ. We have the following bounds:
Lemma 21 (Lemma 38 restated) If ϕ(x; Σ) is the pdf of a d-variate Gaussian with mean 0 and covariance Σ, we have:
where the left-hand-side is understood to be shaped as a vector of dimension R d k 1 +k 2 .
Proof We start with Lemma 20 and use the convexity of the norm
Bounding the right-hand side, we have:
Applying Lemma 46 yields the desired result.
Similarly, we can bound mixed derivatives involving a Laplacian in x:
Lemma 22 If ϕ(x; Σ) is the pdf of a d-variate Gaussian with mean 0 and covariance Σ, we have:
Proof By the definition of a Laplacian, and the AM-GM inequality, we have, for any function
Thus, we have
Applying Lemma 21, the result follows.
Finally, we will need similar bounds for logarithic derivativesÐthat is, derivatives of log p(x), where p is a multivariate Gaussian.
We recall the following result, which is a consequence of the multivariate extension of the Faá di Bruno formula: Proposition 4 (Constantine and Savits (1996), Corollary 2.10) Consider a function f : R d → R, s.t. f is N times differentiable in an open neighborhood of x and f (x) ̸ = 0. Then, for any multi-index I ∈ N d , s.t. |I| ≤ N , we have:
The ≺ ordering on multi-indices is defined as follows: (a 1 , a 2 , . . . As a straightforward corollary, we have the following:
Corollary 23 (Restatement of Lemma 39) For any multi-index I ∈ N d , s.t. |I| is a constant, we have
|I| where J ∈ N d is a multi-index, and J ≤ I iff ∀i ∈ d, J i ≤ I i .
The main strategy in bounding moments of quantities involving a mixture will be to leverage the relationship between the expectation of the score function and the so-called perspective map. In particular, this allows us to bound the moments of derivatives of the mixture score in terms of those of the individual component scores, which are easier to bound using the machinery of Hermite polynomials in the prior section. Note in this section all derivatives are calculated at θ = θ * and therefore p(x, β) = p θ (x, β).
Lemma 24 (Convexity of perspective, Boyd and Vandenberghe ( 2004)) Let f be a convex function.
Then, its corresponding perspective map g(u, v)
We will apply the following lemma many times, with appropriate choice of differentiation operator D and power k.
Lemma 25 (Restatement of Lemma 37) Let D : F 1 → F m be a linear operator that maps from the space of all scalar-valued functions to the space of m-variate functions of x ∈ R d and let θ be such that p = p θ . For k ∈ N, and any norm ∥ • ∥ of interest
Note that since any norm is convex by definition, so is g, by Lemma 24. Then, we proceed as follows:
where (10) follows by linearity of D, and ( 11) by convexity of the function g.
First, we derive the Dirichlet form of Itô diffusions of the form (5). Namely, we show:
Lemma 26 (Dirichlet form of continuous Markov Process) For an Itô diffusion of the form (5), its Dirichlet form is:
Proof By Itô's Lemma, the generator L of the Itô diffusion (5) is:
The Dirichlet form is given by
Expanding and using the definition of Γ, term I can be written as:
We will simplify term II via a sequence of integration by parts:
The term 16 cancels out with term 12, i,j
The term 17 cancels out with the term 14.
For term 13,
Combining term 13 and term 15,
In the end, we are only left with term 18:
We also calculate the integration by parts version of the generalized score matching loss for (3).
Lemma 27 (Integration by parts for the GSM in (3)) The generalized score matching objective in (3) satisfies the equality
Proof Expanding the squares in (3), we have:
The cross-term can be rewritten using integration by parts (under suitable decay at infinity):
Appendix C. A Framework for Analyzing Generalized Score Matching Proposition 5 (Hessian of GSM loss) The Hessian of D GSM satisfies
Proof By a straightforward calculation, we have:
, the second term vanishes at θ = θ * .
By Lemma 2 in Koehler et al. (2022), we also have
which completes the proof.
Proposition 6 (Dirichlet form for CTLD) The Dirichlet form corresponding to CTLD has the form
where E β is the Dirichlet form corresponding to the Langevin diffusion (Lemma 4) with stationary distribution p(x|β).
Proof Equation 19 follows from the fact that CTLD is just a (reflected) Langevin diffusion with stationary distribution p(x, β). Equation 20 follows from the tower rule of expectation and the definition of the Dirichlet form for Langevin from Proposition 4.
Lemma 28 (β derivatives via Fokker Planck) For any distribution p β such that p β = p * N (0, λ min βI) for some p, we have the following PDE for its log-density:
As a consequence, both p(x|β, i) and p(x|β) follow the above PDE.
Proof Consider the SDE dX t = √ 2λ min dB t . Let q t be the law of X t . Then, q t = q 0 * N (0, λ min tI). On the other hand, by the Fokker-Planck equation, d dt q t (x) = λ min ∆ x q t (x). From this, it follows that
Hence, by the chain rule,
Furthermore, by a straightforward calculation, we have
Plugging this in ( 21), we have
as we needed.
Proposition 7 (Integration-by-part Generalized Score Matching Loss for CTLD) The loss D GSM in the integration by parts form (Lemma 1) as:
Moreover, all the terms in the definition of l 1 θ (x, β) and l 2 θ (x, β) can be written as a sum of powers of partial derivatives of ∇ x log p θ (x|β).
r(β)∇ β log p θ (x|β) + ∆ β log p θ (x|β)] + C By Lemma 28, ∇ β log p θ (x|β) is a function of partial derivatives of the score ∇ x log p θ (x|β). Similarly, ∇ 2 β log p θ (x|β) can be shown to be a function of partial derivatives of the score ∇ x log p θ (x|β) as well:
Appendix E. Polynomial mixing time bound: proof of Theorem 12
Proof The proof will follow by applying Theorem 5. Towards that, we need to verify the three conditions of the theorem:
1. (Decomposition of Dirichlet form) The Dirichlet energy of CTLD for p(x, β), by the tower rule of expectation, decomposes into a linear combination of the Dirichlet forms of Langevin with stationary distribution p(x, β|i). Precisely, we have
Proof The proof will follow by an application of a continuous decomposition result (Theorem D.3 in Ge et al. (2018), repeated as Theorem 18) , which requires three bounds:
1. A bound on the Poincaré constants of the distributions p(β|i): since β is independent of i, we have p(β|i) = r(β). Since r(β) is a log-concave distribution over a convex set (an interval), we can bound its Poincaré constant by standard results (Bebendorf, 2003). The details are in Lemma 29, C β ≤ 14D 2 πλ min . 2. A bound on the Poincaré constant C x|β,i of the conditional distribution p(x|β, i): We claim C x|β,i ≤ λ max + βλ min . This follows from standard results on Poincaré inequalities for strongly log-concave distributions. Namely, by the Bakry-Emery criterion, an α-strongly log-concave distribution has Poincaré constant 1 α (Bakry and Émery, 2006). Since p(x|β, i) is a Gaussian whose covariance matrix has smallest eigenvalue lower bounded by λ max
3. A bound on the ªrate of changeº of the density p(x|β, i), i.e.
L ∞ : This is done via an explicit calculation, the details of which are in Lemma 30. By Theorem D.3 in Ge et al. (2018), the Poincaré constant C x,β|i of p(x, β|i) enjoys the upper bound:
which completes the proof.
Lemma 29 (Bound on the Poincaré constant of r(β)) Let C β be the Poincaré constant of the distribution r(β) with respect to reflected Langevin diffusion. Then,
Proof We first show that r(β) is a log-concave distribution. By a direct calculation, the second derivative in β satisfies:
Since the interval is a convex set, with diameter β max , by Bebendorf (2003) we have
from which the Lemma immediately follows.
Lemma 30 (Bound on ªrate of change" of the density p(x|β, i))
We can apply Lemma 28 to derive explicit expressions for the right-hand side:
In 1 , we use (a + b) 2 ≤ 2(a 2 + b 2 ) for a, b ≥ 0; in 2 we apply the moment bound for the Chi-Squared distribution of degree-of-freedom d in Lemma 46; and in 3 we plug in the bound on β max .
Proof The stationary distribution follows from the detailed balance condition w i T (i, j) = w j T (j, i).
We upper bound the Poincaré constant using the method of canonical paths (Diaconis and Stroock, 1991). For all i, j ∈ [K], we set γ ij = {(i, j)} to be the canonical path. Define the weighted length of the path
where the inequality comes from Lemma 17 which provides an upper bound for the chi-squared divergence. Since D is an upper bound and λ min is a lower bound, we may assume without loss of generality that χ 2 max ≥ 1. Finally, we can upper bound the Poincaré constant using Proposition 1 in Diaconis and Stroock (1991)
Next, we will prove a bound on the chi-square distance between the joint distributions p(x, β|i) and p(x, β|j). Intuitively, this bound is proven by showing bounds on the chi-square distances between p(x|β, i) and p(x|β, j) (Lemma 32) Ð which can be explicitly calculated since they are Gaussian, along with tracking how much weight r(β) places on each of the β. Moreover, the Gaussians are flatter for larger β, so they overlap more Ð making the chi-square distance smaller.
Lemma 31 (χ 2 -divergence between joint ªannealedº Gaussians)
Proof Expanding the definition of χ 2 -divergence, we have:
where in Line 22, we apply our Lemma 32 to bound the χ 2 -divergence between two Gaussians with identical covariance. By a change of variable β :=
we can rewrite the integral as:
Since D is an upper bound and λ min is a lower bound, we can assume D 2 λ min ≥ 1 without loss of generality. Plugging in
Finally, we get the desired bound
The next lemma bounds the χ 2 -divergence between two Gaussians with the same covariance.
Lemma 32 (χ 2 -divergence between Gaussians with same covariance)
Proof Plugging in the definition of χ 2 -distance for Gaussians, we have:
In Equation 23, we apply Lemma G.7 from Ge et al. (2018) for the chi-square divergence between two Gaussian distributions. In Equation 24, we use the fact that Σ -1 β is PSD. 3. For every θ * ∈ Θ * and every sufficiently small neighborhood S of θ * , there exists a sufficiently large n, such that there is a unique minimizer θn of Êl θ (x) in S. Furthermore, θn satisfies:
for some matrix Γ SM .
Part 1 is shown in Lemma 34: the claim roughly follows by classic results on the identifiability of the parameters of a mixture (up to permutations of the components) (Yakowitz and Spragins, 1968).
Part 2 is shown in Lemma 36: it follows from a uniform law of large numbers. Finally, Part 3 follows from an application of Lemma 19Ðso we verify the conditions of the lemma are satisfied. The gradient bounds on l θ are verified Lemma 35Ðand it largely follows by moment bounds on gradients of the score derived in Section G. Uniform law of large numbers is shown in Lemma 36, and the the existence of Hessian of L = D GSM is trivially verified.
For the sake of notational brevity, in this section, we will slightly abuse notation and denote D GSM (θ) := D GSM (p, p θ ).
Lemma 34 (Uniqueness of optima) Suppose for θ := (µ 1 , µ 2 , . . . , µ K ) there is no permutation
Proof For notational convenience, let D SM denote the standard score matching loss, and let us denote D SM (θ) := D SM (p, p θ ). For any distributions p θ , by Proposition 1 in Koehler et al. (2022), it holds that
where LSI(q) denotes the Log-Sobolev constant of the distribution q. If θ = (µ 1 , µ 2 , . . . , µ K ) is such that there is no permutation π :
, by Yakowitz and Spragins (1968) we have KL(p θ * , p θ ) > 0. Furthermore, the distribution p θ , by virtue of being a mixture of Gaussians, has a finite log-Sobolev constant (Theorem 1 in Chen et al. (2021)). Therefore, D SM (θ) > D SM (θ * ).
However, note that D GSM (p θ ) is a (weighted) average of D SM losses, treating the data distribution as p β θ * , a convolution of p θ * with a Gaussian with covariance βλ min I d ; and the distribution being fitted as p β θ . Thus, the above argument implies that if θ ̸ = θ * , we have D GSM (θ) > D GSM (θ * ), as we need.
We proceed to bound the right hand sides above. We have:
Where (25) follows by Lemma 25,and (26) follows by combining Corollaries 41 and 23. The same argument, along with Lemma 28, and the fact that max β (∇ β log r(β)) 4 ≲ D 8 λ -4 min by a direct calculation shows that
Lemma 36 (Uniform convergence) The generalized score matching loss satisfies a uniform law of large numbers:
Proof The proof will proceed by a fairly standard argument, using symmetrization and covering number bounds. Precisely, let T = {(x i , β i )} n i=1 be the training data. We will denote by ÊT the empirical expectation (i.e. the average over) a training set T .
We will first show that
from which the claim will follow. First, we will apply the symmetrization trick, by introducing a ªghost training setº
. Precisely, we have:
)
where ( 28) follows by noting the population expectation can be expressed as the expectation over a choice of a (fresh) training set T ′ , (29) follows by applying Jensen's inequality. Next, consider Rademacher variables {ε i } n i=1 . Since a Rademacher random variable is symmetric about 0, we have
For notational convenience, let us denote by
We will bound this supremum by a Dudley integral, along with covering number bounds. Considering T as fixed, with respect to the randomness in {ε i }, the process
In other words, we have
The proof of this is as follows: since ε i is 1-subgaussian, and
The Dudley entropy integral then gives
where N (ϵ, Θ, d) denotes the size of the smallest possible ϵ-cover of the set of parameters Θ in the metric d.
Note that the ϵ in the integral bigger than the diameter of Θ in the metric d does not contribute to the integral, so we may assume the integral has an upper limit
Plugging this estimate back in (31), we get
Taking expectations over the set T (keeping in mind that R is a function of T ), by Lemma 35 we get
This completes the proof of (27). By Markov's inequality, (27) implies that for every ϵ > 0,
Thus, for every ϵ > 0,
Lemma 39 For any multi-index I ∈ N d , s.t. |I| is a constant, we have
For more details, and proof, see Appendix A.4.
First, we need several easy consequences of the machinery developed in Section A.3, specialized to Gaussians appearing in CTLD.
Lemma 40 For all k ∈ N, we have:
where the last inequality follows by Lemma 46.
Combining this Lemma with Lemmas 21 and 22, we get the following corollary:
Finally, we will need the following simple technical lemma:
Lemma 42 Let X be a vector-valued random variable with finite Var(X). Then, we have
where (32) follows from the subadditivity of the spectral norm, (33) follows from the fact that
for any two vectors x, y, as well as the fact that by Jensen's inequality, ∥E
Given this lemma, it suffices to bound E∥(O∇ θ log p θ ) Op θ p θ ∥ 2 2 and E∥(O + O)∇ θ log p θ ∥ 2 2 , which are given by Lemma 43 and Lemma 44, respectively.
Lemma 43
where each µ i is a d-dimensional vector, and we are viewing θ as a dK-dimensional vector.
where the last step follows by Cauchy-Schwartz. To bound both factors above, we will essentially first use Lemma 25 to relate moments over the mixture, with moments over the components of the mixture. Subsequently, we will use estimates for a single Gaussian, i.e. Corollaries 41 and 23.
Proceeding to the first factor, we have:
where (34) follows from the fact that Of = (∇ x f, ∇ β f ) T , (35) follows from Lemma 25, and (36) follows by combining Corollaries 41 and 23 and Lemma 28.
The second factor is handled similarlyfoot_5 . We have:
where (37) follows from Lemma 25, and (38) follows by combining Corollaries 41 and 23 and Lemma 28, as well as the fact that max β (∇ β log r(β)) 4 ≲ D 8 λ -4 min by a direct calculation. Together the estimates ( 36) and ( 38) complete the proof of the lemma.
Lemma 44
where (39) follows by exchanging the order of derivatives, (40) since β is a scalar, so the Laplacian just equals to the Hessian, (41) by dropping the derivatives that are zero in the prior expression.
To bound both summands above, we will essentially first use Lemma 25 to relate moments over the mixture, with moments over the components of the mixture. Subsequently, we will use estimates for a single Gaussian, i.e. Corollaries 23 and 41. Precisely, we have:
of maximum likelihood only when the distribution satisfies a Poincaré inequality. In particular, even simple bimodal distributions in 1 dimension (like a mixture of 2 Gaussians) can significantly worsen the sample complexity of score matching (exponential with respect to mode separation).
For restricted parametric families (e.g. exponential families with sufficient statistics consisting of bounded-degree polynomials), recent work (Pabbaraju et al., 2023) showed that score matching can be comparably efficient to maximum likelihood, by leveraging the fact that a ªrestrictedº form of the Poincaré inequality suffices for good sample complexity. On the empirical side, Song and Ermon (2019) proposed an annealed version of score matching, in which they proposed fitting the scores of the distribution convolved with multiple levels of Gaussian noise. They also proposed using the learned scores to sample via annealed Langevin dynamics, which uses samples from Langevin at higher levels of Gaussian convolution as a warm start for running a Langevin at lower levels of Gaussian convolution. Subsequently, this line of work developed into score-based diffusion models (Song et al., 2020), which can be viewed as a ªcontinuously annealedº version of (Song and Ermon, 2019).
Theoretical understanding of annealed versions of score matching is still very impoverished. A recent line of work (Lee et al., 2022(Lee et al., , 2023;;Chen et al., 2022) explores how accurately one can sample using a learned (annealed) score, if the (population) score loss is successfully minimized. This line of work can be viewed as a kind of ªerror propagationº analysis: namely, how much larger the sampling error with a score learned up to some tolerance. It does not provide insight on when the score can be efficiently learned, either in terms of sample complexity or computational complexity.
Sampling by annealing There are a plethora of methods proposed in the literature that use temperature heuristics (Marinari and Parisi, 1992;Neal, 1996;Earl and Deem, 2005) to alleviate the slow mixing of various Markov Chains in the presence of multimodal structure or data lying close to a low-dimensional manifold. A precise understanding of when such strategies have provable benefits, however, is fairly nascent. Most related to our work, in Ge et al. (2018); Lee et al. (2018), the authors show that when a distribution is (close to) a mixture of K Gaussians with identical covariances, the classical simulated tempering chain (Marinari and Parisi, 1992) with temperature annealing (i.e. scaling the log-pdf of the distribution) mixes in time poly(K).
The mixing time bounds we prove for CTLD rely on decomposition techniques. At the level of the state space of a Markov Chain, these techniques ªdecomposeº the Markov chain by partitioning the state space into sets, such that: (1) the mixing time of the Markov chain inside the sets is good; (2) the ªprojectedº chain, which transitions between sets with probability equal to the probability flow between sets, also mixes fast. These techniques also can be thought of through the lens of functional inequalities, like Poincaré and Log-Sobolev inequalities. Namely, these inequalities relate the variance or entropy of functions to the Dirichlet energy of the Markov Chain: the decomposition can be thought of as decomposing the variance/entropy inside the sets of the partition, as well as between the sets.
Most related to our work are Ge et al. (2018); Moitra and Risteski (2020); Madras and Randall (2002), who largely focus on decomposition techniques for bounding the Poincaré constant. Related ªmultiscaleº techniques for bounding the log-Sobolev constant have also appeared in the literature Otto and Reznikoff (2007); Lelièvre (2009); Grunewald et al. (2009).
Learning mixtures of Gaussians Even though not the focus of our work, the annealed scorematching estimator with the natural parametrization (i.e. the unknown means) can be used to learn the parameters of a mixture from data. This is a rich line of work with a long history. Identifiability of the parameters from data has been known since the works of Teicher (1963); Yakowitz and Spragins (1968). Early work in the theoretical computer science community provided guarantees for clustering-based algorithms (Dasgupta, 1999;Sanjeev and Kannan, 2001); subsequent work provided polynomial-time algorithms down to the information theoretic threshold for identifiability based on the method of moments (Moitra and Valiant, 2010;Belkin and Sinha, 2010); even more recent work tackles robust algorithms for learning mixtures in the presence of outliers (Hopkins and Li, 2018;Bakshi et al., 2022); finally, there has been a lot of interest in understanding the success and failure modes of practical heuristics like expectation-maximization (Balakrishnan et al., 2017;Daskalakis et al., 2017).
Techniques to speed up mixing time of Markov chains SDEs with different choices of the drift and covariance term are common when designing faster mixing Markov chains. A lot of such schemas ªpreconditionº by a judiciously chosen D(x) in the formalism of equation (5). A particularly common choice is a Newton-like method, which amounts to preconditioning by the Fisher matrix (Girolami and Calderhead, 2011;Li et al., 2016;Simsekli et al., 2016), or some cheaper approximation thereof. More generally, non-reversible SDEs by judicious choice of D, Q have been shown to be quite helpful practically (Ma et al., 2015).
ªLiftingº the Markov chain by introducing new variables is also a very rich and useful paradigms. There are many related techniques for constructing Markov Chains by introducing an annealing parameter (typically called a ªtemperatureº). Our chain is augmented by a temperature random variable, akin to the simulated tempering chain proposed by Marinari and Parisi (1992). In parallel tempering (Swendsen and Wang, 1986;Hukushima and Nemoto, 1996), one maintains multiple particles (replicas), each evolving according to the Markov Chain at some particular temperature, along with allowing swapping moves. Sequential Monte Carlo (Yang and Dunson, 2013) is a related technique available when gradients of the log-likelihood can be evaluated.
Analyses of such techniques are few and far between. Most related to our work, Ge et al. ( 2018) analyze a variant of simulated tempering when the data distribution looks like a mixture of (unknown) Gaussians with identical covariance, and can be accessed via gradients to the log-pdf. We compare in more detail to this work in Section 4. In the discrete case (i.e. for Ising models), Woodard et al. (2009b,a) provide some cases in which simulated and parallel tempering provide some benefits to mixing time.
Another way to ªliftº the Markov chain is to introduce a velocity variable, and come up with ªmomentum-likeº variants of Langevin. The two most widely known ones are underdamped Langevin and Hamiltonian Monte Carlo. There are many recent results showing (both theoretically and practically) the benefit of such variants of Langevin, e.g. (Chen and Vempala, 2019;Cao et al., 2023). The proofs of convergence times of these chains is unfortunately more involved than merely a bound on a Poincaré constant (in fact, one can prove that they don't satisfy a Poincaré constant) Ð and it's not so clear how to ªtranslateº them into a statistical complexity analysis using the toolkit we provide in this paper. This is fertile ground for future work, as score losses including a velocity term have already shown useful in training score-based models (Dockhorn et al., 2021).
We will implicitly assume this condition whenever we discuss Poincaré inequalities.
2
The existence of the boundary measure is a standard result of reflecting diffusion processes via solutions to the Skorokhod problem(Saisho, 1987).
In this parametrization, we assume that the weights {wi} K i=1 and shared covariance matrix Σ are known, though the results can be straightforwardly generalized to the natural parametrization in which we are additionally fitting a vector {wi} K i=1 and matrix Σ, at the expense of some calculational complexity.
Of course, in the parametrization in Assumption 2, ∥ΓMLE∥OP itself will generally have dependence on K, which has to be the case since we are fitting Ω(K) parameters.
Note, ∇ β f (β) for f : R → R is a scalar, since β is scalar.