Near-Optimality of Contrastive Divergence Algorithms

Authors: pierre  glaser
kevin han huang
arthur  gretton
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 12:34:56

Abstract: We provide a non-asymptotic analysis of the contrastive divergence (CD) algorithm, a training method for unnormalized models. While prior work has established that (for exponential family distributions) the CD iterates asymptotically converge at an $O(n^{-1 / 3})$ rate to the true parameter of the data distribution, we show that CD can achieve the parametric rate $O(n^{-1 / 2})$. Our analysis provides results for various data batching schemes, including fully online and minibatch. We additionally show that CD is near-optimal, in the sense that its asymptotic variance is close to the Cramér-Rao lower bound.

Main paper: 
 introduction: 
 Describing data using probability distributions is a central task in multiple scientific and industrial disciplines [1, 2, 3] . Since the true distribution of the data is generally unknown, such a task requires finding an estimator of the true distribution among a model class that best describes the available data. An estimator can be characterized at multiple levels of granularity: at the highest level lies consistency [4] , a property which states that as the number of available data points increases, a given estimator will converge to the one best describing the data distribution. At a lower level, a consistent estimator can be further characterized by its convergence rate, a quantity upper-bounding its distance to the true distribution as a function of the number of samples. A convergence rate can be either asymptotic, e.g. hold only in the limit of an infinite sample size, or non-asymptotic, in which case the rate also holds for finite sample sizes. In their simplest form, convergence rates are provided in big-O notation, discarding finer grained information such as asymptotically dominated quantities as well as multiplicative constants. These constants play a role in the so-called asymptotic variance of the estimator, which is a precise descriptor of an estimator's statistical efficiency. Convergence rates and asymptotic variances have been the subject of extensive research in the statistical literature; in particular, well-known lower bounds exists regarding both the best possible (asymptotic) convergence rate of an estimator and its best possible asymptotic variance. These results set a clear frame of reference to interpret individual convergence rates, and are routinely present in the analysis of modern statistical algorithms such as noise-contrastive estimation [5, 6] or score matching [7, 8, 9] .In this work, we focus on cases where (1) the true data distribution admits a density with respect to some known base measure, and (2) the model class is parametrized by a finite-dimensional parameter. In this setting, provided that the true distribution belongs to the model class, a celebrated result in statistical estimation states that the model maximizing the average log-likelihood both achieves the best possible asymptotic convergence rate (called the parametric rate) and the best possible asymptotic variance, called the Cramér-Rao bound (see, e.g. [10] ). While this result shows that Maximum Likelihood Estimators (MLE) are asymptotically optimal, fitting them is complicated by computational hurdles when using models with intractable normalizing constants. Such unnormalized models are common in the Machine Learning literature due to their high flexibility [11, 12] ; their weakness however lies in the fact that expectations under these models have no unbiased approximation. For this reason, popular approximation algorithms such as unbiased gradient-based stochastic optimization of the empirical log-likelihood cannot a priori be used, as the gradient of the normalizing constant is given by an expectation under the model distribution.The Contrastive Divergence (CD) algorithm [13] is a popular approach that circumvents this issue by using a Markov Chain Monte Carlo (MCMC) algorithm to approximate the gradient of the loglikelihood. Unnormalized models trained with Contrastive Divergence have been shown to reach competitive performance in high-dimensional tasks such as image [14, 15, 16] , text [17] , and protein modeling [18, 19] , or neuroscience [20] . A consistency analysis of the Contrastive Divergence algorithm is delicate, however: indeed, the optimization error e.g. the difference between the estimate returned by CD and the MLE, is likely to be non-negligible as compared with statistical error -the distance between the MLE and the true distribution -and thus cannot be discarded, as often done when analyzing estimators that minimize tractable objectives [8, 5] . Recent work [21] elegantly established asymptotic O(n -1/3 )-consistency of the CD estimator for unnormalized exponential families when using only a finite number of MCMC steps. Key to their argument is the fact that the bias of the CD gradient estimate decreases as iterates approach the data distribution. However, as noted by the authors, their work left open the question of whether and under what conditions CD might achieve O(n -1/2 )-consistency.Contributions In this work, we answer this question by providing a non-asymptotic analysis of the CD algorithm for unnormalized exponential families. While existing convergence bounds [21] were derived for the "full batch" setting, where the CD gradient is estimated using the full dataset at each iteration, our analysis covers both the online setting (where data points are processed one at a time without replacement), and the offline setting with multiple data reuse strategies (including full batch).In the online case (Section 3), we show, under a restricted set of assumptions compared to Jiang et al. [21] , that the CD iterates can converge to the true distribution at the parametric O(n -1/2 ) rate. Our analysis reveals that CD contains two sources of approximation: a bias term, and a variance term. These sources are almost independent of each other, in the sense that decreasing the bias by increasing the number of MCMC steps will not decrease the variance. The impact of these two sources of approximation transparently propagates in our resulting bounds: in particular, as the bias of the CD algorithm goes to 0, our bounds recover well-known results in online stochastic optimization [22] . Finally, we study the asymptotic variance of an estimator obtained by averaging the CD iterates, a classic acceleration technique in stochastic optimization [23] . We show that provided that the number of steps m is sufficiently large, the asymptotic variance of this estimator matches (up to a factor 4) the Cramér-Rao bound.Next, we study the offline setting (Section 4), where the CD gradient is estimated by reusing (potentially random) subsets of a finite dataset. We show that a similar result to the online setup holds, up to an additional correlation term that arises from data reuse, and present several approaches to control this term. We improve over the results of [21] by showing a non-asymptotic and nearparametric rate at O((log n) 1/2 n -1/2 ) under their conditions, and also illustrate how different rates can be obtained under a variety of conditions. Our results also show an interesting tradeoff between the effect of initialization and the statistical error as a function of batch size.In summary, we establish the near-optimality of a variety of Contrastive Divergence algorithms for unnormalized exponential families in the so called "long-run" regime, where the number of MCMC steps is high enough to ensure that the CD gradient bias is sufficiently offset by the convexity of the negative log-likelihood. 
contrastive divergence in unnormalized exponential families: 
 Unnormalized Exponential Families Exponential families (EF) [24, 25] form a well-studied class of probability distributions, given by p ψ (dx) := e ψ ⊤ ϕ(x)-log Z(ψ) c(dx), Z(ψ) := X e ψ ⊤ ϕ(x) c(dx).(1)Here, X ∋ x is the data or sample space, which we set to be a subset of R d for some d ∈ N * , although our results are readily extendable to more general measurable spaces. c is a measure on X called the base or carrier measure. When X ⊆ R d , c is often set to be the corresponding Lebesgue measure. ψ ∈ Ψ ⊆ R p is a finite-dimensional parameter called the natural parameter, and ϕ : R d -→ R p is a function called the sufficient statistics, which, alongside with the base measure, fully describes an exponential family. Finally, log Z(ψ), the log-normalizing (or cumulant) function, is a quantity ensuring that p ψ integrates to 1 over X . Crucially, we will not assume that log Z(ψ) admits a closed form expression for all ψ. The latter fact provides the practitioner with a great deal of flexibility in designing the model class: indeed, the only requirement that should be satisfied prior to performing statistical estimation is to have Z(ψ) < +∞ for all ψ, something that can be readily verified and is often the case in practice. The drawback of unnormalized EFs is the fact that sampling (and thus approximating expectations under the model) cannot usually be performed in an unbiased manner. Instead, inference in unnormalized EFs is often performed using tools from the Bayesian Inference literature, such as MCMC [26] . Unnormalized EFs belong to the larger class of unnormalized models [27, 28, 7, 6] , of the form e -E ψ (x)-log Z(ψ) c(dx), Z(ψ) = ∫ e -E ψ (x) c(dx), for some parametrized function E ψ : R d -→ R referred to as the energy. Unnormalized models thus take the flexibility of unnormalized EFs one step further by allowing the (negative) unnormalized log-density to be an arbitrary function E ψ of x and ψ, instead of requiring a linear dependence on ψ as in Equation 1. We focus in this work on unnormalized EFs due to the multiple computational benefits they provide, as explained in the next section, but we believe that extending our analysis to more general unnormalized models is an interesting avenue for future work.Statistical Estimation in Unnormalized Exponential Families using Contrastive Divergence We now review the Contrastive Divergence algorithm, an algorithm used to fit unnormalized models, and our main object of study in this work. The general setting is the following: we assume access to n i.i.d. samples (X 1 , . . . , X n ) drawn from some unknown distribution p ⋆ , which we assume belongs to P ψ , e.g. p ⋆ = p ψ ⋆ for some ψ ⋆ ∈ Ψ. Given these samples, we aim to perform statistical estimation, e.g. find a parameter ψ n within Ψ that should approach ψ ⋆ as n grows.The starting point of the Contrastive Divergence algorithm is the unfortunate realization that Maximum Likelihood Estimation, which corresponds to minimizing the cross-entropy L(ψ) := -E pn log dp ψ /dc between the model p ψ and the empirical data distribution p n := 1/n n i=1 δ Xi , cannot be performed using exact (possibly stochastic) gradient-based optimization, as the gradient ∇ ψ L(ψ) of L with respect to the parameter ψ contains an expectation under the model distribution p ψ . Indeed, the cross entropy and its gradient are given byEQUATIONThe second line follows from the well known identity ∇ ψ log Z(ψ) := E p ψ ϕ; we refer to [25, Proposition 3.1] for a proof. The Contrastive Divergence algorithm circumvents this issue by running approximate stochastic gradient descent (SGD) on L, where the intractable expectation in ∇ ψ log Z is estimated using an MCMC algorithm initialized at the empirical data distribution. In more details, given a number of epochs T , a sequence of data batches B t,j of size B (e.g.B t,j ∈ [[1, n]] B , 1 ≤ t ≤ T, 1 ≤ j ≤ N ⌈n/B⌉), and a family of Markov kernels {k ψ , ψ ∈ Ψ} each with invariant distribution p ψ , at the j th minibatch of epoch t, ∇ ψ log Z(ψ t,j-1 ) is approximated by1 B i∈Bt,j ϕ( Xm i ),where Xm i is produced by running the recursion Xki ∼ k ψt ( Xk-1 i , •), X0 i = X i up to k = m.Throughout the paper, we will refer to the conditional distribution of Xm i given X i as k m ψ (X i , •). The resulting gradient estimate arising from combining this approximation with the other (tractable) sum over the data samples present in ∇ ψ L(ψ), which we refer to as the CD gradient and denote as h t , is thush t,j := 1 B i∈Bt,j ϕ(X i ) - 1 B i∈Bt,j ϕ( Xm i ) = Algorithm 1 Online CD Input: (X 1 , . . . , X n ) i.i.d. ∼ p ψ ⋆ Parameters: Model class {p ψ , ψ ∈ Ψ}, Markov kernels {k ψ , ψ ∈ Ψ}, number of MCMC steps m, learning rate schedule η t := Ct -β , β ∈ [0, 1], C > 0, initial parameter ψ 0 for t = 1, . . . , n do //Approx. sample from p ψt-1 Xm t ∼ k m ψt-1 (X t , •) h t := ϕ(X t ) -ϕ( Xm t ) ψ t ← ψ t-1 -η t h t ψ t ← Proj Ψ (ψ t ) end for return ψ n Algorithm 2 Offline CD Input: (X 1 , . . . , X n ) i.i.d.∼ p ψ ⋆ Parameters: Same as Algorithm 1, plus number of epochs T , batch size B, batching schedule B t,j , initial parameter ψ 0,0 for t = 1, . . . , T do for j = 1, . . . ⌈n/B⌉ do[ Xm t,i,j ∼ k m ψt-1 (X i , •) for i in B t,j ] h t,j := 1 B i∈Bt,j (ϕ(X i )-ϕ( Xm t,i,j )) ψ t,j ← ψ t,j-1 -η t h t,j ψ t,j ← Proj Ψ (ψ t,j )end for end for return ψ T 3 Non-asymptotic analysis of Online CD 
preliminaries and assumptions: 
 Recall that the chi-squared divergence between two probability measures p and q is defined as: χ 2 (p, q) := ( dp dq (x) -1) 2 q(dx) 1/2 if p ≪ q, and +∞ otherwise. Here, p ≪ q denotes that p is absolutely continuous with respect to q and dp/dq is the Radon-Nikodym derivative [29] of p with respect to q. Let L 2 p ψ (R d ) be the space of square-integrable functions with respect to p ψ . For aEQUATIONwhich is a measure of how quick a Markov chain with kernel k ψ mixes, relative to the function f [30] . With these definitions in hand, we now state the assumptions required by our analysis of online CD. These assumptions form a strict subset of the assumptions considered in prior work [21] , which required additional regularity and tail conditions on the Markov kernels k ψ .Assumption A1. P ψ is a subset of a regular and minimal [25, Section 3.2] exponential family with natural parameter domain D ⊆ R p , Ψ is a convex and compact subset of D, and ψ * lies in the interior of Ψ.Assumption A2. There exists a constantC χ > 0 such that χ 2 (p ψ ⋆ , p ψ ) ≤ C 2 χ ∥ψ -ψ ⋆ ∥ 2 Assumption A3. α := sup{α(f, ψ), f ∈ {ϕ i } p i=1 ∪ {ϕ i ϕ j } p i,j=1 , ψ ∈ Ψ} < 1, where ϕ i is the i-th component of the function ϕ, andϕ 2 i is the i-th component of the function x -→ ϕ(x) 2 .A well known property of EFs [25, Proposition 3.1] is that their negative cross-entropy (against any other measure) is C ∞ , convex, and strictly so if the exponential family is minimal (meaning that the set of sufficient statistic functions ϕ i are not linearly dependent). Leaving aside the issue of intractable expectations, this convexity suggests that L can be efficiently minimized using stochastic approximation algorithms [31, 22] . The compactness of Ψ provided by Assumption A1 thus ensures, by the extreme value theorem [32] , the existence of finite positive constants µ and L defined as:EQUATIONwhere ∇ 2 ψ L is the Hessian of L with respect to ψ. µ (called the strong convexity constant) and L (a bound controlling the smoothness of the problem) play a critical role in the analysis of convex optimization algorithms [31] . While it is possible to obtain convergence rates in non-smooth or non-strongly-convex settings, our analysis follows the spirit of [21] by leveraging the strong convexity of the problem to compensate for the bias introduced by using CD gradients instead of unbiased stochastic gradients. Assumption A2 allows to link variations in distribution space to variations in parameter space, and will be instrumental to control the bias of the CD gradient. Note that since χ 2 (p ψ ⋆ , p ψ ) = e log Z(2ψ-ψ ⋆ )-(2 log Z(ψ)-log Z(ψ ⋆ )) -1 provided that 2ψ -ψ ⋆ ∈ D (see [33, Lemma 1] ), we expect Assumption A2 to hold in many cases of interests. On the other hand, the possible exponential scaling of C χ w.r.t log Z suggests that this constant may be large in some instances.Assumption A3 is a restricted spectral gap condition: it guarantees that the time required by the MCMC algorithm to estimate expectations of ϕ and ϕ 2 under p ψ will be uniformly bounded. This assumption is weaker than the (unrestricted) uniform spectral gap condition of [21] , which requires that α controls the convergence rate of all functions in L 2 p ψ (R d ). Note that standard results in stochastic analysis [34] guarantee that α ≤ 1: thus, it only remains to ensure that α is strictly less than 1. Spectral gaps are strongly dependent on two properties of distribution: their tail behavior and their multimodality. While multimodality poses the risk of pushing the constant α close to 1, very heavy tails distributions may not verify the spectral gap condition at all.In this section, we show that under the assumptions stated in Section 3.1, the iterates ψ t produced by the online CD algorithm described in Algorithm 1 will converge to the true parameter ψ ⋆ at the parametric rate O(n -1/2 ). To do so, we follow a well known paradigm in convex optimization [22] by deriving a recursion on the quantity δ t := E ∥ψ t -ψ ⋆ ∥ 2 , which will allow, after unrolling, to obtain convergence rates for the iterates ψ t . We aim to characterize precisely the impact of performing CD as opposed to performing online SGD on L, which would consist of replacing the CD gradient h t of Algorithm 1 by the unbiased (stochastic) gradient, given by:EQUATIONwhich satisfies Eg t (ψ) = ∇ ψ L(ψ). The only stochasticity in g t comes from the sampling of a single data point x t from the true distribution, which is unavoidable in the online setting, and we haveE∥g t (ψ ⋆ )∥ 2 = Cov p ψ ⋆ ϕ =: σ 2 ⋆ . σ 2⋆ plays a key role in the analysis of Stochastic Gradient Descent [22] . We expect that replacing g t by h t will introduce two sources of approximation: a bias term coming from using a finite number of MCMC steps m, and an additional variance term, coming from using a single sample xm t to estimate ∇ ψ log Z(ψ t ). With that in mind, we derive a recursion on δ t in the following lemma. Lemma 3.1. Let (ψ t ) 0≤t≤n be the iterates from Algorithm 1. Denoteδ t = E∥ψ t -ψ ⋆ ∥ 2 , σ ⋆ = (E p ψ ⋆ ∥ϕ -E p ψ ⋆ ϕ∥ 2 ) 1/2, and σ t = (E p ψ t ∥ϕ -E p ψ t ϕ∥ 2 ) 1/2 . Then, under A1, A2 and A3, for all t ≥ 1,EQUATIONwhere ∥log Z∥ 3,∞ is a constant, μm,t := µ -α m σ t C χ , and σm,t := (σ 2 ⋆ + σ 2 t + 2σ 2 t α 2m ) 1/2 .Lemma 3.1 is proved in Appendix D.2, which details the form of ∥log Z∥ 3,∞ , a constant that we expect to scale roughly as dL. Loosely speaking, this recursion suggests that as the learning rate η t goes to 0, the two terms scaling in η 2 t will be negligible, in which case we will have:δ t ≤ (1 -2η t μm,t-1 )δ t-1 < δ t-1, yielding convergence of δ t to 0. We make these arguments formal in the next theorem. The reader familiar with the convex optimization literature will note the similarities between this recursion and the one derived in [22] , which would apply as is to online SGD on L using g t . The difference between the two recursions is that the roles of the strong convexity constant µ and the noise σ ⋆ are now played respectively byμm,t-1 = µ -α m σ t-1 C χ and σ2 m,t-1 = σ 2 ⋆ + σ 2 t + 2σ 2 t α 2mThese two modifications respectively characterize the impact of the bias and the additional variance introduced by the CD gradient. The last term in Equation 7, scaling in α m/2 η 2 t √ δ t , is a residual higher order mixed term coming from relating the variance of the Markov chain sample xm t to σ 2 t . This term can be easily controlled as done next, and disappears as m → ∞. Investigating the impact of m in the recursion, we notice that as m → ∞, μm,t → µ. As we will see later, this ensures that CD will converge for a sufficiently high m. On the other hand, in that same regime, σm,t does not converge to σ ⋆ , but rather to (σ 2 ⋆ + σ 2 t ) 1/2 , showing the irreducible impact of the variance term. While we precisely investigate the impact of the residual variance term in the next section, we now unify σ ⋆ and σ t by introducingEQUATION)σ is an upper bound on the noise induced both by the CD gradient and by the online setup, and was used in prior work [21] . Note that by the properties of log Z, σ 2 also equals sup ψ∈Ψ tr(∇ 2 ψ L(ψ)), where tr(A) is the trace of A ∈ R p×p , and thus finite by the extreme value theorem. The following theorem is obtained by invoking standard unrolling arguments in the convex optimization literature. In the next result, we use the function φ γ (t), defined as φ γ (t) = t γ -1 γ if γ ̸ = 0, and log t if γ = 0. Theorem 3.2. Fix n ≥ 1. Let (ψ t ) 0≤t≤n be the iterates produced by Algorithm 1, and defineδ t := E ∥ψ t -ψ ⋆ ∥ 2 . Moreover, assume that m > log(σCχ/µ) log |α|, i.e. μm := µ -α m σC χ > 0. Then under Assumptions A1, A2 and A3, for η t = Ct -β with C > 0, we have:δ n ≤    2 exp 4 LC 2 φ 1-2β (n) exp -μmC 4 n 1-β δ 0 + σ2 m L2 + 4C σ2 m μmn β , if 0 ≤ β < 1 exp(2 L2 C 2 ) n μC δ 0 + σ2 m L2 + 2σ 2 m C 2 φ μmC/2-1 (n) n μmC/2 , if β = 1 ,where σm = σ 2 (2 + 2α 2m ) + α m/2 ∥log Z∥ 2 3,∞ C 2 χ and L = (L 2 + α m/2 ) 1/2 . Consequently, if η n = Cn with an initial learning rate C > 2μ -1 m , we have√ δ n ≤ 2σ m C μmC μmC-2 1 √ n + o 1 n .Theorem 3.2 is proved in Appendix D. 3 . It shows that the iterates produced by online CD will converge to the true parameter ψ ⋆ at the rate O(n -1/2 ) provided that the number of steps m is sufficiently large, improving over the asymptotic O(n -1/3 ) rate of [21] , while imposing slightly weaker conditions on the number of steps m (see [21, Theorem 2.1] ). This proves that online CD can be asymptotically competitive with other methods for training unnormalized models, such as Noise Contrastive Estimation [6] , or Score Matching [7] . However, the asymptotic variance of ψ t (e.g. the multiplicative factor in front of the O(n -1/2 ) term) is likely to be suboptimal, e.g. much larger than the Crámer-Rao bound, given by the trace of the inverse of the Fisher information matrix [25] . Given the statistical optimality of MLE, and the fact that CD in an approximate MLE method, this motivates the further goal or obtaining a CD estimator with near-optimal statistical properties. In the next section, we achieve this goal by showing that averaging the iterates ψ t will produce a near statistically-optimal estimator, in a sense that we will make precise.Polyak-Ruppert averaging [23] is a simple yet surprisingly effective way to construct an asymptotically optimal estimator ψn := 1 n n t=1 ψ i from a sequence of iterates (ψ t ) 0≤t≤n obtained by running a standard online SGD algorithm [22] . As shown in [22] , when the objective is the cross-entropy of a model, and assuming the unbiased stochastic gradients are available, averaging yields an estimator ψ with the asymptotic variance tr(I(ψ ⋆ ) -1 )/n, where I(ψ) := Cov p ψ ⋆ ϕ is the Fisher information matrix of the data distribution p ψ ⋆ . I(ψ ⋆ ) -1 being the Cramér-Rao lower bound on asymptotic variances of statistical estimators [10] , this estimator ψn is asymptotically optimal. The following theorem shows conditions under which averaging CD iterates can give rise to a near-optimal estimator. Theorem 3.3 (Contrastive Divergence with Polyak-Ruppert averaging). Let (ψ t ) t≥0 the sequence of iterates obtained by running the CD algorithm with a learning rate η t = Ct -β for β ∈ ( 1 2 , 1). Define ψn := 1 n n i=1 ψ i . Then, under the same assumptions as Theorem 3.2, and assuming additionally that m := m(n) > (1-β) log n 2| log α| , we have, for all n ≥ 1,E ψ n -ψ ⋆ 2 1/2 ≤ 2 tr(I(ψ ⋆ ) -1 ) n + o(n -1/2 )Consequently, we have that lim sup n→∞ nE( ψ n -ψ ⋆ 2 ) ≤ 4tr(I(ψ ⋆ ) -1 ).Theorem 3.3, alongside with a statement which includes the asymptotic order of the residual term, is proved in Appendix D. 4 . It shows that at the cost of an increase in computational complexity of the entire algorithm from O(n) to O(n log n), ψn will be a near-optimal statistical estimator of ψ ⋆ . While this increase in complexity emerges from the bias of CD, the additional variance of CD results in an asymptotic variance inflated by a factor of 4 compared to the Cramér-Rao bound.Theorem 3.3 concludes our analysis of online CD. Despite their asymptotic near-optimality, the bounds provided for online CD and its averaged version have weaknesses: the online CD iterates are not robust to choices of C. On the other hand, as shown in Appendix D.4, the bound of the averaged iterates contain higher-order terms that could be large in intermediate sample regimes. Next, we show that offline CD, which processes data points multiple times, can alleviate these issues.4 Non-asymptotic analysis of offline CDIn practice, CD gradient approximation schemes are commonly used within an offline stochastic gradient descent (SGD) algorithm, where one is given the full size-n dataset upfront and each update uses some stochastic subset of the data. We study CD under offline SGD with replacement (SGDw), i.e. Algorithm 2 with batches B t,j being i.i.d. uniform draws of size-B subsets of [n], and include SGD without replacement in Appendix B.2. To do so, we follow the setting of prior work on offline CD [21] , which established its asymptotic O(n -1 3 ) consistency. We show that by slightly strengthening a moment assumption used in [21] , the offline CD iterates converge to the true parameter at a near-parametric O((log n)1 2 n -12 ) rate. Our proof proceeds by controlling a "tail probability" term specific to the offline setting which characterizes the strength of the correlations between the offline CD iterates and the training data. While, as we show, the assumptions of [21] provide a tail control sufficient to obtain a near-parametric rate, other strategies are possible to obtain convergence guarantees. In particular, we show that non-asymptotic convergence can be obtained by either (1) relaxing assumptions on the Markov kernel required by prior work, or (2) making a specific mixing assumption the Markov chain. 
background: asymptotic consistency of offline cd in subexponential settings: 
 Prior work [21] has established asymptotic O(n -13 ) consistency of the (averaged) offline CD iterates in the full-batch case. We summarize their results and assumptions below.Assumption A4. There exists ν ≥ 2 s.t. for all m ∈ N, there is κ ν;m < ∞ s.t.sup x∈X sup ψ∈Ψ E ϕ(K m ψ (x)) -E[ϕ(K m ψ (x))] ν 1/ν ≤ κ ν;m .Assumption A5. There exists some C m > 0 such that, for allψ 1 , ψ 2 ∈ Ψ, sup x∈X ∥E[ϕ(K m ψ1 (x))]- E[ϕ(K m ψ2 (x))]∥ ≤ C m ∥ψ 1 -ψ 2 ∥.Assumption A6. There exist some σ m , ζ m > 0 such that, for any z ∈ R p with ∥z∥ ≤ ζ m ,E[e z ⊤ (ϕ(K m ψ * (X1))-E[ϕ(K m ψ * (X1))])] ≤ e σ 2 m ∥z∥ 2 /2 . Theorem 4.1 (Theorem 2.1 of [21] ). Assume assumptions A1,A2, A3, A4 (for ν = 2), A5 and A6. Let ψ t,1 be the t-th iterate of offline CD with full-batch gradient descent and constant stepsize η t = C, i.e the iterates produced by Algorithm 2 using B t,1 = [ [1, n] ]. Then for any learning rate C and number of Markov kernel steps m satisfying µ -α m σC χ -C 2 (L + α m σC χ ) 2 > 0, we have, for some A m > 0,lim n→∞ P lim sup T →∞ 1 T T t=1 ψ SGDw t,1 -ψ ⋆ > A m n -1 3 = 0This result shows convergence of the averaged full-batch CD iterates to the true parameter in the large n and T limit. As discussed, this result is asymptotic both in n and T : the probability of the error exceeding A m n -1 3 goes to 0 as n → ∞ and T → ∞, but at an unknown rate. Moreover, the O(n -1 3 ) does not match the optimal O(n - for t ≥ 1) SGD with constant step sizes η t = C case, which is the setting of [21] . Analogue bounds holding for the other mentioned batching and step sizes schedules can be found in Appendix B. Theorem 4.2. Assume the setup of Theorem 4.1, except that Assumption A4 holds for some ν > 2, and that μm = µ -α m σC χ > 4CL 2 . Let δ SGDw t,j := E∥ψ SGDw t,j -ψ * ∥ 2 . Then, we have:EQUATIONwhereEQUATIONare functions decreasing exponentially in T , andC ′ (p, ν, m, Ψ) is a constant in n, T . Consequently, lim T →∞ δ SGDw T,1 ≤ e μmC 2 μmC C ′ (p, ν, m, Ψ, β) √ log n √ n + 1 √ n .The precise values of all the constants can be found in Theorem B.1 (for E T,1 1 , E T,1 2 ) and Lemma B.3 (for C ′ (p, ν, m, Ψ, β)), including their expressions for N > 1 and β ∈ [0, 1]. We comment on the main differences between our result and the one of [21] . First our bound holds for any epoch T and number of samples n. Second, fixing n but taking T → ∞, the final bound matches the parametric O( √ n) up to a log(n) factor, a significant improvement over the O(n -1 3 ) rate of [21] . Finally, we control an L 2 error, which is a stronger control than a high probability bound by Markov's inequality; we hypothesize this as the reason why a slightly stronger moment assumption is required for our setup, compared to the one used for the high probability bound in [21] .Inspecting Equation 9, we notice the presence of two transient terms, and a stationary term, reminiscent of the structure of upper bound of Theorem 3.2. The transient terms (i.e. the ones containing E T,11 and E T,12 ) vanish exponentially fast in the total number of CD updates T . However, unlike in online CD where the number of updates and the number of samples are tied (e.g. T = n), these two values are now decoupled, and these terms can be made arbitrarily small by increasing the number of gradient steps T without having to collect more samples n. The stationary term, which is the only one remaining in the limit of T → ∞, decreases with n at a rate that is independent of hyperparameters like the step size C or the learning rate schedule β (see Lemma B.3). In that sense, offline CD compares favorably to online CD, whose rate is sensitive to β and C, and averaged online CD, whose bound contains higher-order (in n) terms which can be large in the moderate n regime. On the other hand, the stationary term in offline CD is asymptotically suboptimal: its rate is larger (while only up to a log factor) than the best-case O( √ n) one achieved by online CD algorithms, and the leading constant does not match the optimal one.The high-level proof of Theorem 4.2 follows a similar strategy as the online one: first, derive a recursion for the quantity δ SGDw t,1:= E∥ψ SGDw t,1 -ψ * ∥ 2 , then unroll it explicitly to obtain a final bound on δ SGDw T,1. The main difference to online CD is the presence of an additional offline-specific correlation between the iterates and the data. We thus break down the proof into three steps: (1) deriving a controllable, uniform-in-time upper bound of the data-iterate correlations, (2) deriving and unrolling a recursion on δ SGDw t,1 containing this new term, and (3) controlling the tail probability term to obtain a final bound on δ SGDw T,1.Step 1:characterizing the data-iterate correlations in offline CD In offline CD, at each epoch t ≥ 1, the iterate ψ SGDw t-1,1 and the data samples X i are correlated: this is because these samples were used in previous epochs t ′ < t -1 to obtain ψ SGDw t ′ ,1 and the data samples X i , which themselves influenced ψ t-1,1 . With such correlations, we now have P(X i |ψ t-1,1 ) ̸ = P(X i ), preventing us from obtaining an unrollable recursion on δ SGDw t,1 by first marginalizing X i out to obtain an upper bound of E ψ SGDw t,1 -ψ * 2 |ψ SGDw t-1,1 that only depends on ψ SGDw t-1,1 -ψ ⋆ , and then marginalizing over ψ SGDw t-1,1 to obtain a recursion as in Lemma 3.1. As this problem would not have occurred had we used "fresh samples" (e.g. i.i.d copies of X i not present in the training data) to perform our update, the core of the proof lies in controlling the following quantity:∆(ψ SGDw t,1 ) := 1 n i≤n E ϕ K m i;ψ SGDw t,1 (X i ) ψ SGDw t,1 , X i -E ϕ K m i;ψ SGDw t,1 (X ′ 1 ) ψ SGDw t,1whereX ′ 1 is an i.i.d. copy of X 1 . ∆(ψ SGDw t,1) is the expected (over the data and iterates) error between a quantity that allows to obtain a recursion (the rightmost term) and the one actually used by offline CD (the leftmost term). To control it, we upper-bound it using a tail decomposition:EQUATIONWe invoke an additional assumption to ensure that (E[∆(ψ SGDw t,1) ν ]) 2/ν is finite; in the results of [21] , this is automatically implied by assumptions A4 and A6. For simplicity we assume the same bounding constant κ ν;m . Assumption A7. There exists ν ≥ 2 s.t. for all m ∈ N, κ ν;m from A4 moreover verifies thatsup ψ∈Ψ E ϕ(K m ψ (X 1 )) -E[ϕ(K m ψ (X 1 ))] ν 1/ν ≤ κ ν;m .Note the similarity of this assumption with assumption A4: the only difference is that X 1 is now random training point instead of an deterministic (arbitrary) one. Ensuring assumption A7 in addition to assumption A4 thus requires controlling a ν-th order moment, instead of all moments as implied by assumption A6.Step 2: Deriving and unrolling the recursion on δ SGDw t,1The right-hand side of Equation ( 10) does not depend on t, allowing for the derivation of an "unrollable" recursion on δ SGDw t,1 and its subsequent unrolling, which is performed in the following theorem. For simplicity, we again assume β = 0 and N = 1 and defer the general case to Theorem B.1 in appendix. Theorem 4.3 (Convergence up to a tail control). Assume A1, A2, A3, A4 and A7. Let η t = C for some C > 0, and assume that μm = µ -α m σC χ > 4CL 2 . Then for any ϵ > 0,δ SGDw T,1 ≤ E T,1 1 δ SGDw 0,0 + C ε SGDw n,m,T ;ν (ϵ) + 5σ + 5κν,m √ n e μm C 2 μmC + E T,1 2 L 2 C 2where ε SGDw n,m,T ;ν (ϵ) is defined in Equation (10) .Note that in the general, non-full batch B ≤ n case, Step 3: Controlling the tail probability term Theorem 4.3 is just one step away from the final bound of Theorem 4.2: it remains to control the tail term ε SGDw n,m,T ;ν (ϵ). Under the assumptions of [21] , minimizing ε SGDw n,m,T ;ν (ϵ) over ϵ yields the following result:Lemma 4.4. Assume the setup of Theorem 4.2. Let n ∈ N be sufficiently large s.t. log n n < σ 2 m ζ 2 m p+ν-2 .Denote r Ψ as the radius of the smallest sphere in R p that contains Ψ, which is finite under A1. Theninf ϵ>0 ε SGDw ν;n,m,T (ϵ) ≤ 3σm p((ν -2)p + 2ν) √ ν -2 + κ ν;m 2 ν-2 2ν (r Ψ ) (ν-2)p 2ν 1 + 2Cm(ν -2) 1/2 σmp 1/2 ((ν -2)p + 2ν) 1/2 (ν-2)p 2ν √ log n √ n.To obtain this result, we control the moment term(sup t E[∆(ψ SGDw t,1) ν ]) using A7, and we control the tail probability term sup t P(∆(ψ SGDw t,1) > ϵ) as in [21, Lemma 3.1] using an union bound, a covering argument and A6. Theorem 4.2 then follows by plugging Lemma 4.4 into Theorem 4.3.As discussed above, the general unrolling result of Theorem 4.3 holds without the subexponentiality assumptions of assumption A6; this assumption was only used in Lemma 4.4 to control ε SGDw n,m,T ;ν (ϵ). We now discuss two alternative ways to control this quantity without requiring subexponential tails. The first generalizes the idea of Jiang et al. [21] , while the second exploits mixing of the Markov chain K m ψ (x) as m → ∞. As before we only state partial results (full batch, β = 0) and defer the full explicit bounds to Appendix B.3.Our first alternative uses Markov Inequality to yield the following. Theorem 4.5. Assume the setup of Theorem 4.3 and additionally that A5 holds. Theninf ϵ>0 ε SGDw ν;n,m,T (ϵ) ≤ C(p, ν, m, Ψ) n - (ν-2)ν 2(ν 2 +(ν-2)p) , lim T →∞ δ SGDw T,1 ≤ C′ (p, ν, m, Ψ) n - (ν-2)ν 2(ν 2 +(ν-2)p) + 1 √ n ,where C and C′ are functions whose explicit expressions are given in Lemma B.4 in the appendix.In the case p = 1 and ν = 3, the sub-optimal error from Theorem 4.5 reads O(n -3/20 ). Theorems 4.2 and 4.5 reveal that, depending on the tail condition imposed on the noise introduced by the Markov kernel, the convergence rate of offline CD varies: A subexponential tail, as assumed in prior work, in fact leads to near-parametric rate. Meanwhile, consistency can be obtained without assuming subexponentiality, albeit at a sub-optimal rate.Control via Markov chain mixing. Alternatively, notice that E[∆(ψ SGDw t,1) 2 ] involves an average of E ϕ K m ψ SGDw t,1 (X i ) 2 X i , ψ SGDw t,1 -E ϕ K m ψ SGDw t,1 (X ′ 1 ) ψ SGDw t,1. When m → ∞, the effect of initialization vanishes, and one may expect the difference to converge to zero. We defer to Lemma B.5 in the appendix to show that, under a ϕ-discrepancy mixing condition ( [35] ) with a mixing coefficientα ∈ [0, 1), inf ϵ>0 ε SGDw ν;n,m,T (ϵ) = O(κ ν;m α (ν-2)m 3ν-2) and limT →∞ δ SGDw T,1 = O κ ν;m α (ν-2)m 3ν-2 + σ + κν;m √ n .As m → ∞, this recovers the parametric rate O(n -1/2 ). This alternative convergence guarantee comes at the cost of requiring m, the number of Markov chain steps, to grow with the sample size n.Remark (Examples). In our main results (Theorems 3.2, 3.3 and 4.3) and the tail condition for offline SGD (Theorem 4.2), we employed a weaker set of assumptions than those in [21] (except for the mild ν > 2 moment assumption in (A4)). Consequently, our results apply to all three examples studied in [21] : A bivariate Gaussian model with unknown mean and random-scan Gibbs sampler, a fully visible Boltzmann machine with random-scan Gibbs sampler, and an exponential-family random graph model with a Metropolis-Hastings sampler. 
related work: 
 Central to this paper is the prior work of Jiang et al. [21] , who provided a rigorous theoretical foundation to analyze the convergence of full-batch CD, and which we refine. The study of optimization with biased gradient descent has attracted a lot of attention in recent years [36, 37, 38, 39] . These works, while closely connected to ours, analyze algorithms with different implementation choices than the CD algorithm: i.i.d. noise setup [36] , or setup where a persistent Markov chain is maintained through the iterations [36, 37, 38, 39] . The latter is akin to a variant of the CD algorithm, called the persistent CD [40] . In contrast, our analysis focus on the CD algorithm that restarts a batch of Markov chains from the data distribution at every iteration. Finally, there is a rich body of work on convergence guarantees for offline multi-pass SGD [41, 42, 43, 44, 45, 46] . A notable difference of our analysis is that we are primarily concerned with statistical errors associated with convergence to the true parameter ψ * in number of samples n, and not the commonly studied convergence rate in number of epochs T . Consequently, most of our work for the offline setup goes into handling the correlations that accumulate by reusing data across epochs. 
discussion: 
 In this work, we provide a non-asymptotic analysis of the Contrastive Divergence algorithms, showing, in the online setting, their potential to converge at the parametric rate and to have near-optimal asymptotic variance, and proving a near-parametric rates in the offline setting, significantly extending prior results. Our results apply to unnormalized exponential families: despite their flexibility, these models only cover log-densities with linear relationships on the model parameters. We believe that extending our results to more general forms of unnormalized models is an important direction for future work.The supplementary material provides the proofs of the main results of the paper:Section B states full explicit bounds for the offline CD algorithm.Section C collects a list of useful tools for our proofs. These include the properties of φ γ introduced before Theorem 3.2 in the main text, as well as several contraction and integrability results.Section D provides the proofs for the online CD algorithm.Sections E, F and G contain the proofs about the offline CD algorithm and the tail control.Throughout the proofs, we will denote by P m ψ the following operator from L 2 p ψ (R) to itself:EQUATIONwhere k m (x, x ′ ) is the m-iterated kernel of the Markov chain with transition kernel k ψ , given by:EQUATIONProj Ψ : R p -→ Ψ denotes the projection operator onto the convex set Ψ, e.g.Proj Ψ (ψ) := argmin ψ ′ ∈Ψ ∥ψ -ψ ′ ∥ .We also frequently use the following function, used in standard convex optimization results [22] .φ γ (t) = t γ -1 γ if γ ̸ = 0 log t if γ = 0 which is defined on R + \ {0}.In this section, we provide the full statements on error bounds for SGD with replacement (SGDw), SGD with reshuffling (SGDo) and tail moment bounds, which complement the results in Section 4. Proofs are deferred to Appendix F, which make use of L 2 approximation by auxiliary gradient updates derived in Appendix E.n,m,T . Denote the SGDw iterates by (ψ SGDw t,j) t∈N,j≤N and let X ′ 1 be an i.i.d. copy of X 1 . Throughout the remaining of the appendix, we define given ϵ > 0 and n ∈ Nϑ SGDw n,m,T (ϵ):= sup t∈[T ] j∈[N ] P n i=1 E ϕ K m ψ SGDw t-1,j (X i ) X i , ψ SGDw t-1,j -E ϕ K m ψ SGDw t-1,j (X ′ 1 ) ψ SGDw t-1,j n > ϵ . = sup t,j P(∆(ψ SGDw t,j ) > ϵ)and ϑ SGDo n,m,T (ϵ) analogously. Using these notations, we can redefine the quantity ε SGDw n,m,T ;ν (ϵ) in the main as ε SGDw n,m,T ;ν (ϵ) := ϵ 2 + κ 2 ν;m ϑ SGDw n,m,T (ϵ) ν-2 ν, consider assumption A7In the result below, we write δ SGDw t,j:= E ψ SGDw t,j-ψ * 2 and, for a fixed ϵ > 0, the quantity is upper bounded byσ SGDw n,T = ε SGDw n,m,T ;ν (ϵ) + 5σ + 5κm √ B = ϵ 2 + κ 2 m ϑ SGDw n,m,T (ϵ) ν-2 ν + 5σ + 5κm √ B . Theorem B.                       E T,N 1 δ SGDw 0,0 + Cσ SGDw n,T 4e μm CN (T +1) 1/2 μmC + 2N (1 + μm C) N -1 φ 1 2 -L 2 C 2 N (T + 1) E T,N 2 for β = 1 2 , E T,N 1 δ SGDw 0,0 + Cσ SGDw n,T 4 μmC + 3N 1 + L 2 C 2 2 N -1 e 2L 2 C 2 N log(T + 1) (T + 1) ( μmCN )/2 for β = 1 , E T,N 1 δ SGDw 0,0 + Cσ SGDw n,T 2 2β+1 μmC e μm C 2(1-β) N (T +1) β + 3 β (1 + μmC) N -1 (T + 2) β L 2 C 2 E T,N 2 otherwise ,where E T,N1 and E T,N 2 are two decreasing functions in T defined byE T,N 1 := exp 1 -N μm Cφ 1-β (T + 1) + N L 2 C 2 2 φ 1-2β (T + 1) , E T,N 2 := exp - N μmC 2 φ 1-β (T + 1) + 2N L 2 C 2 φ 1-2β (T + 1) .We emphasize that the full result above holds for any β ∈ [0, 1], which in particular includes the constant step size β = 0 regime considered by [21] . When β = 0, for E T,N1 and E T,N 2 to decay to zero as T → ∞, we additionally need the conditionμm = µ -α m σC χ > -4CL 2 .This is almost identical to the condition they use ((2.5), Theorem 2.1, [21] ), except that they replace 4L 2 by (L + α m σC χ ) 2 . Notably this says that an additional step size condition is needed for our results to hold in the constant step size regime, but not necessary for a decreasing step size.SGD with reshuffling (SGDo, also called SGD without replacement) is an optimization scheme that is also widely used in practice compared to SGDo and online SGD. In the context of CD, it corresponds to Algorithm 2 with batches chosen as(B t,1 , . . . , B t,N ) = π({1, . . . , n}) ,where π is a uniform draw of the permutation group on n elements. We denote the iterates of SGDo (ψ So t,j ) t∈N,j∈ [N ] . Analogously to ϑ SGDw n,m,T , we define, for X ′1 an i.i.d. copy of X 1 , ϵ > 0 and n ∈ N, the tail probability term ϑ SGDo n,m,T (ϵ):= sup t∈[T ] j∈[N ] P n i=1 E ϕ K m ψ SGDo t-1,j (X i ) X i , ψ SGDo t-1,j -E ϕ K m ψ SGDo t-1,j (X ′ 1 ) ψ SGDo t-1,j n > ϵ . Also denote ε SGDo n,m,T ;ν (ϵ) = ϵ 2 + κ 2 m ϑ SGDo n,m,T (ϵ) ν-2 ν and σ SGDw n,T = ε SGDo n,m,T ;ν (ϵ) + 5σ + 5κm √ B .The following result says that ψ SGDo t,j enjoys exactly the same convergence guarantee as ψ SGDo t,j in Theorem 4.3. The statement is identical to that of Theorem B.1 and is stated in full for completeness; see Appendix F.2 for the proof, which is a slight adaptation of the proof for Theorem B.1. As before we write δ SGDo t,j:= E ψ SGDo t,j -ψ * 2 .Theorem B.2 (Convergence of CD-SGDo). Assume A1 (where Ψ may be non-compact), A2, A3, A4 and A7. Let η t = Ct -β for some β ∈ [0, 1] and C > 0, and assume that m > log(σCχ/µ) log |α| s.t. μm = µ -α m σC χ > 0 as in Theorem 3.2. Then for any ϵ > 0, δ SGDo T,N is upper bounded by                       E T,N 1 δ SGDo 0,0 + Cσ SGDo n,T 4e μm CN (T +1) 1/2 μmC + 2N (1 + μm C) N -1 φ 1 2 -L 2 C 2 N (T + 1) E T,N 2 for β = 1 2 , E T,N 1 δ SGDo 0,0 + Cσ SGDo n,T 4 μmC + 3N 1 + L 2 C 2 2 N -1 e 2L 2 C 2 N log(T + 1) (T + 1) ( μmCN )/2 for β = 1 , E T,N 1 δ SGDo 0,0 + Cσ SGDo n,T 2 2β+1 μmC e μm C 2(1-β) N (T +1) β + 3 β (1 + μmC) N -1 (T + 2) β L 2 C 2 E T,N 2 otherwise ,where E T,N1 and E T,N 2 are two decreasing functions in T defined byE T,N 1 := exp 1 -N μm Cφ 1-β (T + 1) + N L 2 C 2 2 φ 1-2β (T + 1) , E T,N 2 := exp - N μmC 2 φ 1-β (T + 1) + 2N L 2 C 2 φ 1-2β (T + 1) .Remark. We also remark that existing works [47, 48] show that the standard SGDo typically gives a faster convergence rate in T than SGDw. An analogous result for the CD setup would involve additional technical hurdles of jointly controlling the correlations across minibatches and from reusing data samples, and we defer this to future work.We now provide the full explicit tail control bounds. All results in this section hold directly for ε SGDo ν;n,m,T (ϵ) and δ SGDo T,N , and we omit them here. In the result below, we denote r Ψ as the radius of the smallest sphere in R p that contains Ψ, which is finite under A1. Lemma B.3. Assume A5 and A6. Let n ∈ N be sufficiently large s.t. log n n <σ 2 m ζ 2 m p+ν-2 . Then inf ϵ>0 ε SGDw ν;n,m,T (ϵ) ≤ 3σm p((ν -2)p + 2ν) √ ν -2 + κ ν;m 2 ν-2 2ν (r Ψ ) (ν-2)p 2ν 1 + 2Cm(ν -2) 1/2 σmp 1/2 ((ν -2)p + 2ν) 1/2 (ν-2)p 2ν √ log n √ n .In particular, if we additionally assume the conditions of Theorem B.1, we havelim T →∞ δ SGDw T,N ≤ C ′ (p, ν, m, Ψ, β) √ log n √ n + 1 √ BwhereC ′ (p, ν, m, Ψ, β) := 8(1 + 5σ + 5κm) μm × 3σm p((ν -2)p + 2ν) √ ν -2 + κ ν;m 2 ν-2 2ν (r Ψ ) (ν-2)p 2ν 1 + 2Cm(ν -2) 1/2 σmp 1/2 ((ν -2)p + 2ν) 1/2 (ν-2)p 2ν.Lemma B.4. Assume the conditions of Theorem 4.3 and additionally that A5 holds. Theninf ϵ>0 ε SGDw ν;n,m,T (ϵ) ≤ 3C m + κ ν/2 ν;m (r Ψ ) (ν-2)p 2ν C -ν-2 2 m 3 (ν-2)p 2ν n - (ν-2)ν 2(ν 2 +(ν-2)p) , lim T →∞ δ SGDw T,N 1/2 ≤ C′ (p, ν, m, Ψ, β) n - (ν-2)ν 2(ν 2 +(ν-2)p) + B -1/2 ,whereC′ (p, ν, m, Ψ) := 8(1 + 5σ + 5κm) μm 3C m + κ ν/2 ν;m (r Ψ ) (ν-2)p 2ν C -ν-2 2 m 3 (ν-2)p 2ν.The next result considers a ϕ-discrepancy mixing condition ( [35] ), which is a mixing assumption on K m ψ but with respect to a specific test function ϕ, and we impose it uniformly over ψ ∈ Ψ. We also recall that X ψ 1 ∼ p ψ . Lemma B.5. Assume that there exist α ∈ [0, 1) and CK > 0 such that, for all ψ ∈ Ψ and x ∈ X ,∥E[ϕ(K m ψ (x))] -E[ϕ(X ψ 1 )]∥ ≤ CK αm . Then inf ϵ>0 ε SGDw ν;n,m,T (ϵ) ≤ 1 + 2 ν-2 2ν κ ν;m ( CK ) ν-2 2ν α (ν-2)m 3ν-2.In particular, if we additionally assume the conditions of Theorem 4.3, we havelim T →∞ δ SGDw T,N ≤ 8 µ -α m σCχ 1 + 2 ν-2 2ν κ ν;m ( CK ) ν-2 2ν α (ν-2)m 3ν-2 + 5σ + 5κm √ B . C Auxiliary Tools C.1 Properties of φ γThe following lemma collects some identities used in [22] . Lemma C.1. φ γ satisfies the following properties:(i) φ γ is increasing on R + for all γ ;(ii) φ γ (t) ≤ t γ γ for γ > 0, andφ γ (t) ≤ -1 γ for γ < 0 ; (iii) φ 1-β (t) ≥ t 1-β for β ∈ (0, 1] ; (iv) φ γ (t) -φ γ ( t 2 ) ≥ 1 2 x γ for γ ∈ (0, 1] .The next lemma provides some additional results on φ γ . Lemma C.2. φ γ satisfies the following properties:(i) φ γ is positive on t > 0 and increasing for every γ ∈ R;(ii) For 1 ≤ t 1 ≤ t 2 and β ≥ 0, we haveφ 1-β (t 2 + 1) -φ 1-β (t 1 ) ≤ t2 t=t1 t -β ≤ 2 φ 1-β (t 2 + 1) -φ 1-β (t 1 ) . If instead β < 0, we have 1 2 φ 1-β (t 2 + 1) -φ 1-β (t 1 ) ≤ t2 t=t1 t -β ≤ φ 1-β (t 2 + 1) -φ 1-β (t 1 ) ; (iii) For 1 ≤ t 1 ≤ t 2 and γ ̸ = 0, we have t γ-1 1 ≤ φ γ (t 2 ) -φ γ (t 1 ) ≤ t γ-1 2 if γ ≥ 1 , t -(1-γ) 2 ≤ φ γ (t 2 ) -φ γ (t 1 ) ≤ t -(1-γ) 1 if γ ≤ 1 ; (iv) Let 1 ≤ t 1 < t 2 and κ, β ≥ 0. If κ ̸ = 1 and a > 0, we have t2 t=t1 (t + 1) -β exp a φ 1-κ (t -1) ≤ (t 2 + 1) max{κ-β,0} a exp a φ 1-κ (t 2 + 1) ,and if κ ̸ = 1 and a < 0, we havet2 t=t1 (t + 1) -β exp a φ 1-κ (t) ≤ (t 2 + 1) max{κ-β,0} (-a)exp a φ 1-κ (t 1 ) .Proof of Lemma C.2. (i) follows from checking γ > 0, γ = 0 and γ < 0 respectively. The first set of bounds in (ii) follow by noting that t → t -β is decreasing for β ≥ 0:t2 t=t1 t -β ≥ t2+1 t1 t -β dt = φ 1-β (t 2 + 1) -φ 1-β (t 1 ) , t2 t=t1 t -β ≤ 2 t2 t=t1 (t + 1) -β ≤ 2 t2 t1-1 (t + 1) -β dt = 2 φ 1-β (t 2 + 1) -φ 1-β (t 1 ) .The second set of bounds follows from noting that t → t -β is increasing for β < 0:t2 t=t1 t -β ≥ 1 2 t2 t=t1 (t + 1) -β ≥ 1 2 t2 t1-1 (t + 1) -β dt = 1 2 φ 1-β (t 2 + 1) -φ 1-β (t 1 ) , t2 t=t1 t -β ≤ t2+1 t1 t -β dt = φ 1-β (t 2 + 1) -φ 1-β (t 1 ) .For (iii), we note that for γ ̸ = 0,φ γ (t 2 ) -φ γ (t 1 ) = t γ 2 -t γ 1 γ, so by the mean value theorem,inf t1≤t≤t2 t γ-1 ≤ φ γ (t 2 ) -φ γ (t 1 ) ≤ sup t1≤t≤t2 t γ-1 .The desired bounds then follow from an explicit computation of the infimum and the maximum in each of the two cases γ ≥ 1 and γ ≤ 1.For (iv), we first consider the case κ ̸ = 1 and a > 0. Thent2 t=t1 (t + 1) -β exp a φ 1-κ (t -1) = e -a 1-κ t2 t=t1 (t + 1) -β exp at 1-κ 1 -κ ≤ e -a 1-κ max t1≤t≤t2 (t + 1) κ (t + 1) β t2 t=t1 (t + 1) -κ exp at 1-κ 1 -κ ≤ e -a 1-κ (t 2 + 1) max{κ-β,0} t2 t=t1 (t + 1) -κ exp at 1-κ 1 -κ . Since, for x ≥ 0, x → (x + 1) -κ is decreasing and x → exp(ax 1-κ /(1 -κ)) is increasing, we have that for t 1 ≤ t ≤ t 2 and x ∈ [t, t + 1], (t + 1) -κ ≤ x -κ and exp(at 1-κ /(1 -κ)) ≤ exp(ax 1-κ /(1 -κ)) .This implies that t2 t=t1 (t + 1) -β exp a φ 1-κ (t -1)≤ (t 2 + 1) max{κ-β,0} e -a 1-κ t2 t=t1 t+1 t x -κ exp ax 1-κ 1 -κ dx = (t 2 + 1) max{κ-β,0} e -a 1-κ t2+1 t1 x -κ exp ax 1-κ 1 -κ dx ≤ (t 2 + 1) max{κ-β,0} a e -a 1-κ e a(t 2 +1) 1-κ 1-κ = (t 2 + 1) max{κ-β,0} a exp a φ 1-κ (t 2 + 1) .The main difference in the case κ ̸ = 1 and a < 0 is that we now use x → exp(a(x + 1)1-κ /(1 -κ)) is decreasing to obtain, for t 1 ≤ t ≤ t 2 and x ∈ [t, t + 1], exp(a(t + 1) 1-κ /(1 -κ)) ≤ exp(ax 1-κ /(1 -κ)) .A similar argument then yields t2 t=t1 (t + 1) -β exp a φ 1-κ (t)≤ (t 2 + 1) max{κ-β,0} e -a 1-κ t2 t=t1 (t + 1) -κ exp a(t + 1) 1-κ 1 -κ ≤ (t 2 + 1) max{κ-β,0} e -a 1-κ t2+1 t1 x -κ exp ax 1-κ 1 -κ dx = (t 2 + 1) max{κ-β,0} a exp a φ 1-κ (t 2 + 1) -exp a φ 1-κ (t 1 ) ≤ (t 2 + 1) max{κ-β,0} (-a)exp a φ 1-κ (t 1 ) .We also need the following lemma, which is useful for controlling the accumulation of errors from the noise terms over iterations.Lemma C.3. For any a, b ≥ 0, T, N ∈ N and κ, β ≥ 0 such that bt -β -at -κ ≤ 1 for all 1 ≤ t ≤ T , we have thatT t=1 (1 -bt -β + at -κ ) N ≤ exp -bN φ 1-β (T + 1) + aN φ 1-κ (T + 1) .Moreover, for any ζ ≥ 0, we have thatT t=1 t -ζ N j=1 (1 -bt -β + at -κ ) T s=t+1 (1 -bs -β + as -κ ) N ≤ Q 1 + exp - bN 2 φ 1-β (T + 1) + 4aN φ 1-κ (T + 1) Q 2 ,whereEQUATIONQ 2 :=    3 ζ (1 + a) N -1 2a (T + 2) max{κ-ζ,0} if κ ̸ = 1 and a > 0 , 2N (1 + a) N -1 φ 1-ζ-2aN (T + 1) if κ = 1 or a = 0 .In the special case ζ = β = 1 < κ , we haveEQUATION.Proof of Lemma C.3. By assumption, bt -β -at -κ ≤ 1 for all 1 ≤ t ≤ T . Since 0 ≤ 1 -x ≤ e -x for all x ≤ 1, we have that for anyEQUATION)Applying this to the first quantity of interest followed by noting that a, b ≥ 0 and using Lemma C.2(ii), we obtain the first bound thatT t=1 (1 -bt -β + at -κ ) N ≤ exp -bN T t=1 t -β + aN T t=1 t -κ ≤ exp -bN φ 1-β (T + 1) + aN φ 1-κ (T + 1) .For the second bound, we definet 0 := sup t ≤ T b 2 ≤ at -(κ-β) .Then by noting that 1 -bt -β + at -κ ≥ 0 for all 1 ≤ t ≤ T again,we can bound the quantity of interest asT t=1 t -ζ N j=1 (1 -bt -β + at -κ ) j-1 T s=t+1 (1 -bs -β + as -κ ) N = T t=t0+1 t -ζ N j=1 (1 -bt -β + at -κ ) j-1 T s=t+1 (1 -bs -β + as -κ ) N + T s=t0+1 (1 -bs -β + as -κ ) × t0 t=1 t -ζ N j=1 (1 -bt -β + at -κ ) j-1 t0 s=t+1 (1 -bs -β + as -κ ) N ≤ T t=t0+1 t -ζ N j=1 1 - b 2 t -β j-1 T s=t+1 1 - b 2 s -β N + T s=t0+1 1 - b 2 s -β N t0 t=1 t -ζ N j=1 (1 + at -κ ) j-1 t0 s=t+1 (1 + as -κ ) N ≤ N × T t=t0+1 t -ζ T s=t+1 1 - b 2 s -β N =:S1 + N (1 + a) N -1 × T s=t0+1 1 - b 2 s -β N =:S3 × t0 t=1 t -ζ t0 s=t+1 (1 + as -κ ) N =:S2.In the last line, we have used that 0 ≤ 1 -b 2 t -β ≤ 1 for t ≥ t 0 + 1 and 1 + at -κ ≤ 1 + a. To control the three quantities, we first note that by ( 13) , we haveS 3 ≤ exp - bN 2 T s=1 s -β exp bN 2 t0 s=1 s -β (a) ≤ exp - bN 2 T s=1 s -β exp aN t0 s=1 s -κ (b) ≤ exp - bN 2 φ 1-β (T + 1) + 2aN φ 1-κ (T + 1) .In (a) above, we have noted that b 2 ≤ as -(κ-β) for s ≤ t 0 ; in (b), we have used t 0 ≤ T and Lemma C.2(ii) with a, b ≥ 0. In the special case β = 1 < κ, the above yieldsS 3 ≤ (T + 1) -bN 2 exp 2aN 1 -(T + 1) -(κ-1) κ -1 ≤ (T + 1) -bN 2 e 2aN κ-1 .We now control S 2 . By (13) again, we haveS 2 ≤ t0 t=1 t -ζ exp aN t0 s=t+1 s -κ ≤ T t=1 t -ζ exp aN T s=t+1 s -κ (c) ≤ exp 2aN φ 1-κ (T + 1) × T t=1 t -ζ exp -2aN φ 1-κ (t + 1) (d) ≤ 3 ζ exp 2aN φ 1-κ (T + 1) T t=1 (t + 2) -ζ exp -2aN φ 1-κ (t + 1) .In (c) above, we have applied Lemma C.2(ii); in (d), we have noted that sup t∈N (t + 2) β /t β = 3 β . If κ ̸ = 1 and a > 0, we can apply Lemma C.2(iv) to get thatS 2 ≤ 3 ζ 2aN (T + 2) max{κ-ζ,0} exp 2aN φ 1-κ (T + 1) = Q 2 N (1 + a + c) N -1 exp 2aN φ 1-κ (T + 1) .If κ = 1 or a = 0, the bound from (c) above readsS 2 ≤ exp 2aN φ 1-κ (T + 1) T t=1 t -ζ (t + 1) -2aN ≤ exp 2aN φ 1-κ (T + 1) T t=1 t -ζ-2aN ≤ 2 φ 1-ζ-2aN (T + 1) exp 2aN φ 1-κ (T + 1) = Q 2 N (1 + a) N -1 exp 2aN φ 1-κ (T + 1), where we have used Lemma C.2(ii) in the last line. Now consider the special case with ζ = 1 < κ, the bound from (d) becomesS 2 ≤ 3 exp 2aN φ 1-κ (T + 1) T t=1 (t + 2) -1 exp -2aN φ 1-κ (t + 1) ≤ 3 exp 2aN 1 -(T + 1) -(κ-1) κ -1 T t=1 (t + 2) -1 ≤ 3e 2aN κ-1 log(T + 1) .We are left with controlling S 1 , which follows from a similar strategy as controlling S 2 :EQUATION≤ 4 ζ exp - bN 2 φ 1-β (T + 1) T t=1 (t + 3) -ζ exp bN 2 φ 1-β (t + 1) .In (a) above, we used Lemma C.2(ii). For β ̸ = 1 and b ̸ = 0, we can apply Lemma C.2(iv) with b 2 > 0 to obtainS 1 ≤ 4 ζ exp - bN 2 φ 1-β (T + 1) (T + 3) max{β-ζ,0} bN/2 exp bN 2 φ 1-β (T + 3) = 2 2ζ+1 (T + 3) max{β-ζ,0} bN exp bN 2(1 -β) (T + 3) 1-β -(T + 1) 1-β (b) ≤ 2 2ζ+1 (T + 3) max{β-ζ,0} bN exp bN 2(1 -β)(T + 1) β = Q 1 N .In (b), we have used Lemma C.2(iii) with 1 -β ≤ 1. Meanwhile, if β = 1 or b = 0, we haveS 1 ≤ exp - bN 2 φ 1-β (T + 1) T t=1 t -ζ (t + 1) bN/2 ≤ exp - b 2 φ 1-β (T + 1) T t=1 t -ζ+bN/2 ≤ 2φ 1-ζ+bN/2 (T + 1) exp - bN 2 φ 1-β (T + 1) = Q 1 N .For the special case with ζ = β = 1, the bound in (14) becomesS 1 ≤ exp - bN 2 φ 0 (T + 1) T t=1 t -1 exp bN 2 φ 0 (t + 1) = (T + 1) -bN 2 T t=1 t -1 (t + 1) bN 2 ≤ (T + 1) -bN 2 T t=1 t -(1-bN 2 ) (c) ≤ 2(T + 1) -bN 2 φ bN 2 (T + 1) = 2(T + 1) -bN 2 (T + 1) bN/2 -1 bN/2 ≤ 4 bN.In (c), we have used Lemma C.2(ii) for both the case 1 -bN 2 ≤ 0 and 1 -bN 2 ≥ 0.Combining the bounds for the general cases, we obtain the first desired inequality thatT t=1 t -ζ N j=1 (1 -bt -β + at -κ ) j-1 T s=t+1 (1 -bs -β + as -κ ) N ≤ Q 1 + exp - b 2 φ 1-β (T + 1) + uφ 1-ξ (T + 3) + 4a φ 1-κ (T + 1) Q 2 .For the special case ζ = β = 1 < κ, γ, combining the earlier bounds givesT t=1 t -ζ N j=1 (1 -bt -β + at -κ ) j-1 T s=t+1 (1 -bs -β + as -κ ) N ≤ 4 b + 3N (1 + a) N -1 e 4aN κ-1 log(T + 1) (T + 1) bN 2.The next result is a standard result in convex analysis, needed to handle projections performed in Algorithms 1 and 2.Lemma C.4. Let Ψ a be convex subset of R p . Let ψ ⋆ ∈ Ψ. Then, for all ψ ∈ R p , we have:∥Proj Ψ (ψ) -ψ ⋆ ∥ ≤ ∥ψ -ψ ⋆ ∥Proof. We have:∥ψ -ψ ⋆ ∥ 2 = ∥ψ -Proj Ψ (ψ) + Proj Ψ (ψ) -ψ ⋆ ∥ 2 = ∥ψ -Proj Ψ (ψ)∥ 2 + 2 ⟨ψ -Proj Ψ (ψ), Proj Ψ (ψ) -ψ ⋆ ⟩ + ∥Proj Ψ (ψ) -ψ ⋆ ∥ 2Since by [49, Proposition 1.1.9], we have:⟨ψ -Proj Ψ (ψ), ψ ′ -Proj Ψ (ψ)⟩ ≤ 0for all ψ ′ ∈ Ψ, we can use this inequality at ψ ′ = ψ ⋆ ∈ Ψ to obtain:∥ψ -ψ ⋆ ∥ 2 ≥ ∥ψ -Proj Ψ (ψ)∥ 2and the result follows by taking the square root.We now state two lemmas that guarantee an amount of integrability sufficient to our analysis.Lemma C.5. Let p, q ∈ P(X ) such that dp dq exists, and such that χ 2 (p; q) < +∞. Assume thatf ∈ L 2 q (R) Then |E p f | < +∞.Proof. By assumption, f ∈ L 2 q (R d ). Moreover, χ 2 (p, q) < +∞, and thus we have dp dq -1 ∈ L 2 q (R). Thus, the inner product is finite, andf dp dq -1 dq = f dp -f dq = |E p f -E q f | := M < +∞ =⇒ M -|E q f | < E p f < M + |E q f | Lemma C.6.For all ψ ∈ Ψ, for all m ≥ 1, and for all k ≥ 1, we have:E p ψ ⋆ P m ψ ϕ k < +∞Proof. By analycity of the log partition function ψ -→ log Z(ψ), we have that ∥ϕ∥ k dp ψ (x) < +∞ for all ψ ∈ Ψ, and thus, the function x -→ ∥ϕ∥ k (x) ∈ L 2 p ψ (R) for all ψ. Consequently,P m ∥ϕ∥ k ∈ L 2 p ψ (R d ).We can apply lemma C.5 to P m ∥ϕ∥ k to obtain E ψ ⋆ P m ∥ϕ∥ k < +∞ for all k ≥ 1 and for all m ≥ 0. As a byproduct, we obtain that P m ∥ϕ∥ k ∈ L 2 p ψ ⋆ (R), and thus so it ∥P m ϕ∥ k .The following lemma is used multiple time in our analysis. Lemma C.7. For any function f ∈ L 2 p ψ ⋆ (R d ) ∩ L 2 p ψ (R d ), we have:∥E p ψ ⋆ P m ψ (f -E p ψ f )∥ ≤ α m C χ (E p ψ f -E p ψ f 2 ) 1/2 ∥ψ -ψ ⋆ ∥ Proof. Let us note first that ∥E p ψ ⋆ P m ψ (f -E p ψ f )∥ 2 = p i=1 P m ψ (f i -E p ψ f i )(x)(p ψ ⋆ (x) -p ψ (x))dx 2 = p i=1 P m ψ (f i -E p ψ f i )(x)( p ψ ⋆ p ψ -1)p ψ (x)dx 2 ≤ χ 2 (p ψ , p ψ ⋆ ) p i=1 P m ψ (f i -E p ψ f i )(x) 2 p ψ (x)dx ≤ χ 2 (p ψ , p ψ ⋆ ) p i=1 α m (f i -E p ψ f i )(x) 2 p ψ (x)dx ≤ α m C 2 χ ∥ψ -ψ ⋆ ∥ E p ψ f -E p ψ f 2 .Lemma C.8. Let f : Ψ → R p be a differentiable function in the interior of Ψ ⊆ R p . For ψ ∈ Ψ, define σ min (ψ) := inf θ∈Ψ,∥θ∥=1 θ ⊤ ∇f (ψ)θ and σ max (ψ) := sup θ∈Ψ,∥θ∥=1 θ ⊤ ∇f (ψ)θ with respect to the Jacobian matrix ∇f (ψ). Then for any ψ 1 , ψ 2 ∈ Ψ, we have thatinf ψ∈Ψ σ min (ψ) ≤ (ψ 1 -ψ 2 ) ⊤ f (ψ 1 ) -f (ψ 2 ) ≤ sup ψ∈Ψ σ max (ψ)Proof of Lemma C.8. By the mean value theorem, there exists some a ∈ (0, 1) such that(ψ 1 -ψ 2 ) ⊤ f (ψ 1 ) -f (ψ 2 ) = (ψ 1 -ψ 2 ) ⊤ f (1 × ψ 1 + 0 × ψ 2 ) -f (0 × ψ 1 + 0 × ψ 2 ) = (ψ 1 -ψ 2 ) ⊤ ∇f (aψ 1 + (1 -a)ψ 2 )(ψ 1 -ψ 2 ) = ∥ψ 1 -ψ 2 ∥ 2 (ψ 1 -ψ 2 ) ⊤ ∥ψ 1 -ψ 2 ∥ ∇f (aψ 1 + (1 -a)ψ 2 ) ψ 1 -ψ 2 ∥ψ 1 -ψ 2 ∥.Plugging in the definition of σ max gives the desired upper bound and similarly σ min implies the lower bound.We recall the following notations used in the next lemmas, namely σ ψ := E p ψ ∥ϕ -E p ψ ϕ∥ 2 , σ ⋆ := σ ψ ⋆ and σ := sup ψ∈Ψ σ ψ .We now provide two intermediary lemmas necessary to analyze the impact of variance in the CD gradient. The strategy is similar in both of them: we change the integration from p ψ ⋆ to p ψ to obtain contraction, at the cost of an additional term scaling with C χ ∥ψ -ψ ⋆ ∥. We obtain exact constants that we choose to describe in terms of the smoothness parameters of the problem, e.g. the k th derivatives of the log partition function log Z, which, for k ≥ 2, equals the k th derivative of the negative cross-entropy model w.r.t p ψ ⋆ .The following lemmas guarantee the second moment of a sample from k m ψ p ψ ⋆ approaches the second moment of a sample from the target distribution p ψ . Lemma D.1. Under A1, A2 and A3, for all ψ ∈ Ψ, we have:|E p ψ ⋆ P m ψ ∥ϕ∥ 2 -E p ψ ∥ϕ∥ 2 | ≤ α m C χ ∥ψ -ψ ⋆ ∥ ∥log Z∥ 1,∞ where ∥log Z∥ 1,∞ := sup ψ∈Ψ p i=1 (4∂ 1 i log Z(ψ) 2 ∂ 2 i log Z(ψ) + 2∂ 2 i log Z(ψ) 2 + 4∂ 1 i log Z(ψ)∂ 3 i log Z(ψ) + ∂ 4 i log Z(ψ)) 1/2 < +∞ Proof. Applying Lemma C.7 to each f i := ϕ 2 i , we have |E p ψ ⋆ P m ψ ϕ 2 i -E p ψ ϕ 2 i | = α m C χ ∥ψ -ψ ⋆ ∥ E p ψ ϕ 2 i -E p ψ ϕ 2 i 2 1/2 = α m C χ ∥ψ -ψ ⋆ ∥ E p ψ ϕ 4 i -(E p ψ ϕ 2 i ) 2 1/2We map the two moments to derivatives of log Z(ψ), since the k th derivative of log Z(ψ) is the k th cumulant. It can be shown, using the multivariate moment to cumulant mapping, thatEϕ 4 i = ∂ log Z ∂ψ i 4 + 6 ∂ log Z ∂ψ i 2 ∂ 2 log Z ∂ψ 2 i + 3 ∂ 2 log Z ∂ψ 2 i 2 + 4 ∂ log Z ∂ψ i ∂ 3 log Z ∂ψ 3 i + ∂ 4 log Z ∂ψ 4 i = ∂ 1 i log Z(ψ) 4 + 6∂ 1 i log Z(ψ) 2 ∂ 2 i log Z(ψ) + 3∂ 2 i log Z(ψ) 2 + 4∂ 1 i log Z(ψ)∂ 3 i log Z(ψ) + ∂ 4 i log Z(ψ) where ∂ k i log Z(ψ)denotes the k th derivative of log Z with respect to ψ i . On the other hand,Eϕ 2 i = ∂ 1 i log Z(ψ) 2 + ∂ 2 i log Z(ψ) =⇒ (Eϕ 2 i ) 2 = ∂ 1 i log Z(ψ) 4 + 2∂ 1 i log Z(ψ) 2 ∂ 2 i log Z(ψ) + ∂ 2 i log Z(ψ) 2 implying E p ψ ϕ 4 i -(E p 2 ψ ϕ 2 i ) 2 = 4∂ 1 i log Z(ψ) 2 ∂ 2 i log Z(ψ) + 2∂ 2 i log Z(ψ) 2 + 4∂ 1 i log Z(ψ)∂ 3 i log Z(ψ) + ∂ 4 i log Z(ψ)The result follows by summing over i, since:|E p ψ ⋆ P m ψ ∥ϕ∥ 2 -E p ψ ∥ϕ∥ 2 | ≤ d i=1 |E p ψ ⋆ P m ψ ϕ 2 i -E p ψ ϕ 2 i | Note that ∥log Z∥ 1,∞ is finite since Ψ is compact and log Z is analytic.The next lemma provides convergence (in squared absolute value) of the first moment of the m-iterated markov kernel k m ψ . Lemma D.2. Under A1, A2 and A3, for all ψ ∈ Ψ, we have|E p ψ ⋆ ∥P m ψ ϕ∥ 2 -E p ψ ϕ 2 | ≤ α m σ 2 ψ + C χ α m/4 ∥log Z∥ 2,∞ ∥ψ -ψ ⋆ ∥ where ∥log Z∥ 2,∞ = sup ψ∈Ψ p i=1 F (ψ)∂ 2 i log Z(ψ) 1/4 + 2 ∂ 1 i log Z(ψ)∂ 2 i log Z(ψ) 1/2 and F (ψ) := 15∂ 2 i log Z(ψ) 3 + 10∂ 3 i log Z(ψ) 2 + 15∂ 2 i log Z(ψ)∂ 4 i log Z(ψ) + ∂ 6 i log Z(ψ)Proof. We have:E p ψ ⋆ (P m ψ ϕ i ) 2 = E p ψ ⋆ (P m ψ ϕ i -E ψ ϕ i + E ψ ϕ i ) 2 = E p ψ ⋆ (P m ϕ i -E p ψ ϕ i ) 2 + 2E ψ ϕ i × E ψ ⋆ P m (ϕ i -E ψ ϕ i ) + (E ψ ϕ i ) 2 =⇒ |E p ψ ⋆ (P m ψ ϕ i ) 2 -(E ψ ϕ i ) 2 | ≤ |E p ψ ⋆ (P m (ϕ i -E p ψ ϕ i )) 2 | ∆1 +2|E ψ ϕ i E ψ ⋆ P m (ϕ i -E ψ ϕ i )| ∆2 where ∆ 1 = E p ψ (P m ψ (ϕ i -E ψ ϕ i ∆1,1)) 2 + (E p ψ ⋆ (P m ψ ϕ i -E ψ ϕ i ) 2 -E p ψ (P m ψ ϕ i -E ψ ϕ i ) 2 ∆1,2 ) ≤ α 2m E ψ (ϕ i -E ψ ϕ i ) 2 + C χ ∥ψ -ψ ⋆ ∥ E p ψ (P m (ϕ i -E p ψ ϕ i ) 4 1/2 ≤ α 2m E ψ (ϕ i -E ψ ϕ i ) 2 + C χ ∥ψ -ψ ⋆ ∥ E p ψ P m (ϕ i -E p ψ ϕ i ) 2 1/4 E p ψ P m (ϕ i -E p ψ ϕ i ) 6 1/4 ≤ α 2m E ψ (ϕ i -E ψ ϕ i ) 2 + α m/2 C χ ∥ψ -ψ ⋆ ∥ E p ψ (ϕ i -E p ψ ϕ i ) 2 1/4 E p ψ (ϕ i -E p ψ ϕ i ) 6 1/4where we used the fact that P m ψ is a contraction in L 6 p ψ (R d ) and the L 2 p ψ (R d )-α spectral gap of P m ψ . As an aside, note that a simpler result can be obtained by making regularity assumption on the mapping ψ -→ P m ψ . Assuming that ψ -→ P m ψ (x) is uniformly L m -Lipschitz across x ∈ X for instance (as done in [21, Assumption 5] ), the second term ∆ 1,2 of ∆ 1 could have been handled using∆ 1,2 ≤ 2E p ψ ⋆ P m ψ ϕ -E ψ ⋆ ϕ 2 + 2 E p ψ ϕ -E p ψ ⋆ ϕ 2 ≤ 4(L m ∥ψ -ψ ⋆ ∥ + σ 2 ⋆ α 2m + 2 E p ψ ϕ -E p ψ ⋆ ϕ 2 ) ≤ 4(L m ∥ψ -ψ ⋆ ∥ + σ 2 ⋆ α 2m + 2L 2 ∥ψ -ψ ⋆ ∥)Although this result does not require possibly large constants related to sixth-order moments, it is less tight in the sense that it does not go to 0 as m → ∞. Back to the main proof, and ∆ 2 in particular. Applying Lemma C.7 to f := ϕ, we have∆ 2 ≤ α m E ψ ϕ i C χ ∥ψ -ψ ⋆ ∥ E p ψ (ϕ i -E p ψ ϕ i ) 2 1/2 ≤ α m C χ ∂ 1 i log Z(ψ)∂ 2 i log Z(ψ) 1/2 ∥ψ -ψ ⋆ ∥Putting everything together, we have:|E p ψ ⋆ (P m ψ ϕ i ) 2 -E p ψ ϕ 2 i | ≤ α 2m E p ψ (ϕ i -E p ψ ϕ i ) 2 + C χ ∥ψ -ψ ⋆ ∥ α m/2 (E p ψ (ϕ i -log ϕ i ) 6 ) 1/4 (E p ψ (ϕ i -log ϕ i ) 2 ) 1/4 + 2α m ∂ 1 i log Z(ψ)∂ 2 i log Z(ψ) 1/2 ≤ α 2m E p ψ (ϕ i -E p ψ ϕ i ) 2 + C χ α m/2 ∥ψ -ψ ⋆ ∥ (E p ψ (ϕ i -log ϕ i ) 6 ) 1/4 (E p ψ (ϕ i -log ϕ i ) 2 ) 1/4 + 2 ∂ 1 i log Z(ψ)∂ 2 i log Z(ψ) 1/2Summing over i, we obtain|E p ψ ⋆ ∥P m ψ ϕ∥ 2 -(E p ψ ϕ) 2 | ≤ p i=1 |E p ψ ⋆ (P m ψ ϕ i ) 2 -(E p ψ ϕ i ) 2 | ≤ α 2m p i=1 E p ψ (ϕ i -E p ψ ϕ i ) 2 + C χ α m/2 ∥log Z∥ 2,∞ ∥ψ -ψ ⋆ ∥ ≤ α 2m σ 2 ψ + C χ α m/2 ∥log Z∥ 2,∞ ∥ψ -ψ ⋆ ∥ where ∥log Z∥ 2,∞ = sup ψ∈Ψ p i=1 E p ψ (ϕ i -E p ψ ϕ i ) 6 1/4 E p ψ (ϕ i -E p ψ ϕ i ) 2 1/4 + 2 ∂ 1 i log Z(ψ)∂ 2 i log Z(ψ) 1/2Similarily to the previous lemma, one can upper bound E(ϕ i -E p ψ ϕ i ) 6 using the centered moment to cumulant formula:E p ψ (ϕ i -Eϕ i ) 6 = 15∂ 2 i log Z(ψ) 3 + 10∂ 3 i log Z(ψ) 2 + 15∂ 2 i log Z(ψ)∂ 4 i log Z(ψ) + ∂ 6 i log Z(ψ) =: F (ψ)To get a full description of ∥log Z∥ 2,∞ :∥log Z∥ 2,∞ = sup ψ∈Ψ p i=1 F (ψ)∂ 2 i log Z(ψ) 1/4 + 2 ∂ 1 i log Z(ψ)∂ 2 i log Z(ψ) 1/2 .We can now use the previous lemmas to obtain an expression on the second moment of the contrastive divergence gradient estimator, relating it to the one of the stochastic log-likelihood gradient estimator. Lemma D.3. Under A1, A2 and A3, we have:E ∥h t (ψ, X t )∥ 2 ≤ 2σ 2 ⋆ + 2σ 2 ψ + 2L 2 ∥ψ -ψ ⋆ ∥ 2 + 4(σ 2 ψ α 2m + α m/2 ∥log Z∥ 3,∞ C χ ∥ψ -ψ ⋆ ∥) where ∥log Z∥ 3,∞ := 2 max(∥log Z∥ 1,∞ , ∥log Z∥ 2,∞ ).Proof. We rely on the following decomposition:EQUATIONEQUATION(1a) + (1b) form the differentiable stochastic gradient g t of Equation 6. Note that (1a) is meanzero, and ( 2) is mean-zero conditionally on X t . Consequently, E ⟨(2), (3EQUATION⟩ = 0, and the only mixed-terms that remain to be controlled are E ⟨(1a), (3)⟩ and E ⟨(1b), (3)⟩. We first control the unmixed terms, and the simple ones first: we haveE ∥(1a)∥ 2 = σ 2 ⋆ , as well as E ∥(1b)∥ 2 ≤ L 2 ∥ψ -ψ ⋆ ∥ 2 .For (2), we have:E k m ψ (x,•) ∥(2)∥ 2 = E k m ψ (x,•) ∥ϕ(k m ψ (x, •))∥ 2 -∥E k m ψ (x,•) ϕ(k m ψ (x, •))∥ 2 = P m ψ ∥ϕ(x)∥ 2 -∥P m ψ ϕ(x)∥ 2We can invoke Lemmas D.2 and D.1, which guarantee|E p ψ ⋆ ∥P m ψ ϕ∥ 2 -E p ψ ϕ 2 | ≤ α 2m σ 2 ψ + C χ α m/2 ∥log Z∥ 2,∞ ∥ψ -ψ ⋆ ∥ |E p ψ ⋆ P m ψ ∥ϕ∥ 2 -E p ψ ∥ϕ∥ 2 | ≤ α m C χ ∥log Z∥ 1,∞ ∥ψ -ψ ⋆ ∥ to obtain E p ψ ⋆ E k m ψ (x,•) ∥(2)∥ 2 = E p ψ ∥ϕ∥ 2 -∥Eϕ∥ 2 + α 2m σ 2 ψ + C χ α m/2 ∥log Z∥ 1,∞ + ∥log Z∥ 2,∞ ∥ψ -ψ ⋆ ∥ = E ∥ϕ -Eϕ∥ 2 + α 2m σ 2 ψ + α m/2 ∥log Z∥ 3,∞ C χ ∥ψ -ψ ⋆ ∥ where ∥log Z∥ 3,∞ := 2 max(∥log Z∥ 1,∞ , ∥log Z∥ 2,∞ ).For (3), notice that (3) is precisely the term ∆ 1 in Lemma D.2, and we can thus bound it byE p ψ ⋆ ∥(3)∥ 2 ≤ σ 2 ψ α 2m + α m/2 ∥log Z∥ 2,∞ C χ ∥ψ -ψ ⋆ ∥ Finally, we simply bound 2E ⟨(1a), (3)⟩ by E ∥(1a)∥ 2 + E ∥(3)∥ 2 , and 2E ⟨(1b), (3)⟩ by E ∥(1b)∥ 2 + E ∥(3)∥ 2 .Putting everything together, we have:E ∥h t (ψ)∥ 2 = E ∥(1a)∥ 2 + E ∥(1b)∥ 2 + E ∥(2)∥ 2 + E ∥(3)∥ 2 + 2E ⟨(1a), (3)⟩ + 2E ⟨(1b), (3)⟩ ≤ 2E ∥(1a)∥ 2 + 2E ∥(1b)∥ 2 + E ∥(2)∥ 2 + 3E ∥(3)∥ 2 ≤ 2σ 2 ⋆ + 2L 2 ∥ψ -ψ ⋆ ∥ 2 + σ 2 ψ + 4(σ 2 ψ α 2m + α m/2 ∥log Z∥ 3,∞ C χ ∥ψ -ψ ⋆ ∥) ≤ 2σ 2 ⋆ + 2σ 2 ψ + 2L 2 ∥ψ -ψ ⋆ ∥ 2 + 4(σ 2 ψ α 2m + α m/2 ∥log Z∥ 3,∞ C χ ∥ψ -ψ ⋆ ∥) D.2 Proof of the SGD recursion (Lemma 3.1)We are now ready to provide a SGD-style recursion for the expected squared distance to the optimumE ∥ψ t -ψ ⋆ ∥ 2 .Lemma (Restatement of Lemma 3.1). Let ψ t be the iterates produced by Algorithm 1. Denoteδ t = E∥ψ t -ψ ⋆ ∥ 2 , σ ⋆ = (E p ψ ⋆ ∥ϕ -E p ψ ⋆ ϕ∥ 2 ) 1/2, and σ t = (E p ψ t ∥ϕ -E p ψ t ϕ∥ 2 ) 1/2 . Then, under Assumptions A1, A2 and A3, for all t ≥ 1, we have:δ t ≤ (1 -2η t μm,t-1 + 2η 2 t L 2 )δ t-1 + η 2 t σ2 m,t-1 + 4α m/2 η 2 t ∥log Z∥ 2,∞ C χ δ 1/2 t-1where ∥log Z∥ 3,∞ is a constant, μm,t := µ -α m σ t C χ , and σm,t := (σ 2 ⋆ + σ 2 t + 2σ 2 t α 2m ) 1/2 .Proof. In this proof, we note (F t ) t≥0 , the increasing family of σ-algebras generated by the random variables (X t ) t≥0 ∼ p ψ ⋆ and the markov chain samples Xm t |X t , ψ t ∼ k m ψt (X t , •). We decompose the integrand of δ t as follows:∥ψ t -ψ ⋆ n ∥ 2 = ∥Proj Ψ (ψ t-1 -η t h t (ψ t-1 )) -ψ ⋆ n ∥ 2 ≤ ∥ψ t-1 -η t h t (ψ t-1 ) -ψ ⋆ n ∥ 2 (By Lemma C.4) = ∥ψ t-1 -ψ ⋆ ∥ 2 -2η t ⟨h t (ψ t-1 ), ψ t-1 -ψ ⋆ n ⟩ + η 2 t ∥h t (ψ t-1 )∥ 2The first term is (up to an averaging operation) the previous iterate. The middle term will ensure (provided m is large enough) contraction of the expected distance to the optimum. Finally, the third term can be described by Lemma D.3, and anssentially behaves like the second moment of a log-likelihood stochastic gradient. Indeed, noting g(ψ) := -E p ψ ⋆ ϕ + E p ψ ϕ the expectation of g t w.r.t x t (which is the gradient of the negative cross-entropy between p ψ and p ψ ⋆ ), we have:⟨h t (ψ n,t-1 ), ψ t-1 -ψ ⋆ n ⟩ = ⟨g(ψ n,t-1 ), ψ t-1 -ψ ⋆ n ⟩ + ⟨(h t (ψ n,t-1 ) -g(ψ n,t-1 )), ψ t-1 -ψ ⋆ n ⟩ ∆and applying Lemma C.7, we get thath t (ψ t-1 ) -g(ψ t-1 ) = ϕ(k m (X t , •)) -E ψ ϕ =⇒ E ψ ⋆ E k m ψ (Xt,•) [h t (ψ t-1 ) -g(ψ t-1 )|F t-1 ] = P m ψ ϕ -E ψ ϕ , meaning |E [∆|F t-1 ]| ≤ E p ψ ⋆ P m ψ (ϕ -E ψ ϕ) ∥ψ t-1 -ψ ⋆ ∥ ≤ α m C χ (E p ψ t-1 f -E p ψ t-1 f 2 ) 1/2 ∥ψ -ψ ⋆ ∥ 2 ≤ α m σ t-1 C χ ∥ψ -ψ ⋆ ∥ 2On the other hand, by applying Lemma C.8 to g, we have:⟨g(ψ t-1 ), ψ t-1 -ψ ⋆ ⟩ = ⟨g(ψ t-1 ) -g(ψ ⋆ ), ψ t-1 -ψ ⋆ ⟩ ≥ µ ∥ψ t-1 -ψ ⋆ ∥ 2Combining the above results, we obtain:E p ψ ⋆ E k m ψ (x,•) [∥ψ t -ψ ⋆ ∥ |F t-1 ] ≤ (1 -2η t (µ -α m σ t-1 C χ ))∥ψ t-1 -ψ ⋆ ∥ 2 + η 2 t (2σ 2 ⋆ + 2σ 2 t-1 + 2L 2 ∥ψ t-1 -ψ ⋆ ∥ 2 + 4(σ 2 t-1 α 2m + α m/2 ∥log Z∥ 3,∞ C χ ∥ψ t-1 -ψ ⋆ ∥)And the result follows by integrating over F t-1 .We now prove Theorem 3.2. The recursion of Lemma 3.1 is (almost) identifiable, up to a crossterm of second order, with the one of an SGD algorithm as presented in the setting of [22] . To make the identification exact, we use the bound4α m/2 η 2 t ∥log Z∥ 3,∞ C χ δ 1/2 t-1 ≤ 2α m/2 η 2 t δ t + 2α m/2 η 2 t ∥log Z∥ 2 3,∞ C 2 χ, yielding the following recursion:δ t ≤ (1 -2η t (µ -α m σC χ ) + 2η 2 t (L 2 + α m/2 ))δ t-1 + (σ 2 (2 + 2α 2m ) + α m/2 ∥log Z∥ 2 3,∞ C 2 χ )η 2 twhere we used the fact that σ m,t ≤ σ. We can use arguments existing in prior work [22, Theorem 1] , noting that the only condition required to proceed is that µ -α m σC χ < (L 2 + 1) 1/2 , which automatically holds since µ < L. We then obtain, notingμm = µ -α m σC χ , σm = σ 2 (2 + 2α 2m ) + α m/2 ∥log Z∥ 2 3,∞ C 2 χ , and L = (L 2 + α m/2 ) 1/2 , δ n ⩽    2 exp 4 LC 2 φ 1-2β (n) exp -μmC 4 n 1-β δ 0 + σ 2 L 2 + 4C σ2 m µn β , if 0 ⩽ β < 1 exp(2 L2 C 2 ) n μC δ 0 + σ2 m L2 + 2σ 2 m C 2 φ μmC/2-1 (n) n μmC/2 , if β = 1.We start by three intermediate lemmas.Lemma D.4. Under Assumptions A1, A2, A3, the online CD iterates produced by Algorithm 1 usingη t = Ct -β for β ∈ ( 1 2 , 1) verify 1 n n i=1 E ∥ψ i -ψ ⋆ ∥ 1/2 1/2 ≤ P 3 (n) where P 3 (n) is a term of order n -1 2 -β 4 .Proof. Let us note δ n = E ∥ψ n -ψ ⋆ ∥ 2 For the second term, specific to our setting, we have, recalling that δ n is bounded thanks to Theorem 3.2, Recalling that δ n verifies the bound of Theorem 3.2, we have:1 n 2 n i=1 δ 1/2 i ≤ 2C 1/2 σ 2µ 1/2 n 2 n i=1 i -β/2 + √ 2 n 2 δ 0 + σ2 m L2 n i=1 e 2 L2 C 2 φ 1-2β (i) e -μm C 8 i 1-β A3 ≤ 2C 1/2 σm n 2 μ1/2 m φ 1-β/2 (n) + √ 2 n 2 δ 0 + σ2 m L 2 A 3 =⇒ 1 n n i=1 δ 1/2 i 1/2 ≤ 1 n 2C 1/2 σm μ1/2 φ 1-β/2 (n) 1/2 + 1 n √ 2(δ 0 + σ 2 L2 )A 3 1/2where A 3 is finite if β < 1, and A 3 = O(n) otherwise [22] .Lemma D.5. Under Assumptions A1, A2, A3, the online CD iterates produced by Algorithm 1 usingη t = Ct -β for β ∈ ( 1 2 , 1) verify E 1 n n i=1 h i (ψ i-1 ) 2 1/2 ≤ P 2 (n)for some term P 2 of order at most n β 2 -1 (an order smaller than n -1 2 ).Proof. The result follows from the fact that h t verifiesh t (ψ t-1 ) = 1 ηt (ψ t-1 -ψ t ) .A similar quantity was handled in the case of standard SGD [22, Theorem 3] , and the only condition needed to reuse their steps is that (ψ t ) t≤n satisfies an upper bound of the same form as the one [22, Theorem 1] derived. This is precisely the nature of our bound of ψ t estblished in Theorem 3.2, with μm , σm , L. Borrowing on their result, we have:E 1 n n i=1 h i (ψ i-1 ) 2 1/2 ≤ 4σmβ C 1/2 nμm φ β/2 (n) + 4β Cnμ 1/2 m δ 0 + σ2 m L2 1/2 A 2 + 1 nμ 1/2 m 1 C + 2 L δ 1/2 0 + 2 L nμ 1/2 m 2C 1/2 σm μ1/2 m φ 1-β (n) 1/2 + 4 L nμ 1/2 m δ 0 + σ2 m L2 1/2 A 1/2 2where μm , σm and L are defined in 3.2, andA 2 = n k=1 e -μC 16 k 1-β +16 L4 1 C 4 φ 1-2β (k)Lemma D.6. Under Assumptions A1, A2, A3, the online CD iterates produced by Algorithm 1 usingη t = Ct -β for β ∈ ( 1 2 , 1) verify 1 n n i=1 (E ∥ψ i -ψ ⋆ ∥ 4 ) 1/2 ≤ P 1 (n)where P 1 (n) is a term in n of order O(n -β ).Theorem (Contrastive Divergence with Polyak-Ruppert averaging). Let t ) t≥0 the sequence of obtained by running the CD algorithm with a learning rate η t = Ct -β for β ∈ (0, 1). Define ψn := 1 n n i=1 ψ i . Then, under the same assumptions as Theorem 3.2 we have, for all n ≥ 1, and assuming m >(1-β) log n 2| log α| E ψ n -ψ ⋆ 2 1/2 ≤ 2tr(I(ψ) -1 ) √ n + P (n)I(ψ ⋆ ) := Cov p ψ ⋆ ϕ is the Fisher information matrix of the data distribution, and the termP (n) is of order O(n -max(-β, β 2 -1,-( β 2 +m |log α| log n )) ).Proof. Let h n be the standard online (biased) CD gradient defined in Equation 3:h n (ϕ) = -ϕ(x n ) + ϕ(k m ψn-1 (•, x n )), x n ∼ p ψ ⋆ And let us note h(ϕ) = -E p ψ ⋆ ϕ(x n ) + E ψ ⋆ E k m ψ n-1 ϕ(k m ψn-1 (x n , •)), which corresponds to a CD gradient the negative cross entropy between the model and the true (and not empirical) data distribution, defined as:f (ψ) = E p ψ ⋆ ϕ(x i ) ⊤ ψ + log Z(ψ) (15) It holds that: f ′′ (ψ ⋆ )(ψ n-1 -ψ ⋆ ) = f ′ (ψ n-1 ) -f ′ (ψ ⋆ ) + (f ′′ (ψ ⋆ )(ψ n-1 -ψ ⋆ ) -f ′ (ψ n-1 ) + f ′ (ψ ⋆ )) = h n (ψ n-1 ) -f ′ (ψ ⋆ ) + (f ′′ (ψ ⋆ )(ψ n-1 -ψ ⋆ ) -f ′ (ψ n-1 ) + f ′ (ψ ⋆ )) + (f ′ (ψ n-1 ) -h(ψ n-1 )) + (h(ψ n-1 ) -h n (ψ n-1 ))and thus, we have that1 n n i=1 ψ i -ψ ⋆ = 1 n n i=1 f ′′ (ψ ⋆ ) -1 h i (ψ i-1 ) (i) + 1 n n i=1 f ′′ (ψ ⋆ ) -1 (f ′′ (ψ ⋆ )(ψ i-1 -ψ ⋆ ) -f ′ (ψ i-1 ) (ii) + 1 n n i=1 f ′′ (ψ ⋆ ) -1 (f ′ (ψ i-1 ) -h(ψ i-1 )) (iii) + 1 n n i=1 f ′′ (ψ ⋆ ) -1 (h(ψ i-1 ) -h i (ψ i-1 )) (iv).We seek to obtain a bound on E1 n n i=1 ψ i -ψ ⋆ 2 1/2. (i) and (ii) will be bounded using techniques from prior work [22] . (iii) captures the bias of the CD algorithm, while (iv) captures the variance.Bounding (iii) Note that the first difference (the bias term) is deterministic. Moreover, using lemma C.7, we have:f ′ (ψ n-1 ) -h(ψ n-1 ) = E p ψ n-1 ϕ -E p ψ ⋆ E k m ψ n-1 ϕ ≤ α m (E p ψ n-1 ϕ -E p ψ n-1 ϕ 2 ) 1/2 C χ ∥ψ n-1 -ψ ⋆ ∥ ≤ α m σC χ ∥ψ n-1 -ψ ⋆ ∥By Minkowski's inequality, we have:(E(iii) 2 ) 1/2 ≤ α m Cχ n n i=1 (E ∥ψ i-1 -ψ ⋆ ∥ 2 ) 1/2 = α m Cχ n n i=1 δ 1/2 i-1Recalling that δ i satisfies theorem 3.2, we have that δ1 2i is dominated by n -β/2 , and we thus thus haven i=1 δ 1/2 i which is dominated by n 1-β 2 . Overall, we obtain that (E(iii) 2 ) 1/2 is dominated by α m n -β/2 . Noting that α m n -β/2 = n -( β 2 +m |log α| log n ) ,this term is dominated by n -1/2 since, by assumption, we have m > (1-β) log n 2| log α| . To track the final order of the bound, we write (E(iii) 2 ) 1/2 = P 4 (n), whereP 4 (n) = O(n -( β 2 +m |log α| log n ) ).Bounding (iv) For (iv), we start by further breaking down the We havef (ψ ⋆ ) -1 (h(ψ ) -h n (ψ n-1 )) = f ′′ (ψ ⋆ ) -1 (ϕ(x n ) -E p ψ ⋆ ϕ) ∆1,n +f ⋆ (ψ ⋆ ) -1 (ϕ(k m ψn-1 (x n )) -E p ψ ⋆ E k m ψ n-1 ϕ(k m ψn-1 (x n , •))) ∆2,n = f ′′ (ψ ⋆ ) -1 ∆ 1,n + f ′′ (ψ ⋆ ) -1 ∆ 2,nand by Minkowski's inequality, we haveEQUATIONNote that this step was made possible because we are looking at the square-root of the variance, which is unlike the recursion in Lemma 3.1. This allows to separate the terms and use fewer intermediaries than in the proof of Lemma 3.1.Since both ∆ 1,n and ∆ 2,n are martingale differences with respect to the filtration F n-1 , the covariance terms vanish, and we have(E(iv) 2 ) 1/2 ≤ 1 n n i=1 E f ′′ (ψ ⋆ ) -1 ∆ 1,n 2 1/2 + 1 n n i=1 E f ′′ (ψ ⋆ ) -1 ∆ 2,n 2 1/2 ≤ 1 n n i=1 tr(f ′′ (ψ ⋆ ) -1 E∆ 1,n ∆ ⊤ 1,n f ′′ (ψ ⋆ ) -1 ) 1/2 + 1 n n i=1 tr(f ′′ (ψ ⋆ ) -1 E∆ 2,n ∆ ⊤ 2,n f ′′ (ψ ⋆ ) -1 ) 1/2which is key to obtain the desired acceleration. Since E∆ 1,n ∆ ⊤ 1,n = I(ψ ⋆ ), the first sum equalstr(I(ψ ⋆ ) -1 ) √ n. For the second term, since ∆ 2,n is mean-zero conditional on F n-1 , we have∆ 2,n = ϕ(k m ψn-1 (x n , •)) -E p ψ ⋆ E k m ψ n-1 ϕ(k m ψn-1 (x n , •)) =⇒ E ∆ 2,n ∆ ⊤ 2,n |F t-1 = E p ψ ⋆ E k m ψ n-1 ϕ(k m ψn-1 (x n , •))ϕ(k m ψn-1 (x n , •)) ⊤ |F t-1 -(E p ψ ⋆ E k m ψ n-1 [ϕ(k m ψn-1 (•, x n ))|F t-1 ])(E p ψ ⋆ E k m ψ n-1 ϕ(k m ψn-1 (•, x n ))|F t-1 ) ⊤ = E p ψ ⋆ P m ψn-1 ϕϕ ⊤ |F t-1 -(E p ψ ⋆ [P m ψn-1 ϕ|F t-1 ])(E p ψ ⋆ [P m ψn-1 ϕ ⊤ |F t-1 ])where we abused the notation for P m ψn-1 to handle functions mapping to R d and R d×d instead of only R. We now show that ∆ 2,n ∆ ⊤ 2,n approximately equals I(ψ n-1 ), by inspecting the second moment E p ψ ⋆ P m ψn-1 ϕϕ ⊤ and the squared first moment (E p ψ ⋆ P m ψn-1 ϕ)(E p ψ ⋆ P m ψn-1 ϕ) ⊤ . Looking at the second moment first, we have, omitting the conditioning on F t-1 for brevity:E p ψ ⋆ P m ψn-1 ϕϕ ⊤ = (E p ψ ⋆ P m ψn-1 ϕϕ ⊤ -E p ψ n-1 ϕϕ ⊤ ) + E p ψ n-1 ϕϕ ⊤ = E p ψ ⋆ P m ψn-1 (ϕϕ ⊤ -E p ψ n-1 ϕϕ ⊤ ) ∆3,n +E p ψ n-1 ϕϕ ⊤ since E p ψ n-1 P m ψn-1 (ϕϕ ⊤ -E p ψ n-1 ϕϕ ⊤ ) = 0, we have ∆ 3,n = E p ψ ⋆ P m ψn-1 (ϕϕ ⊤ -E p ψ n-1 ϕϕ ⊤ ) -E p ψ n-1 P m ψn-1 (ϕϕ ⊤ -E p ψ n-1 ϕϕ ⊤ ) =⇒ ∥∆ 3,n ∥ F ≤ α m C χ (E p ψ n-1 ϕϕ ⊤ -E p ψ n-1 ϕϕ ⊤ 2 F ) 1/2 ∥ψ n-1 -ψ ⋆ ∥which follows from applying lemma C.7 to the function ϕ i ϕ j -E p ψ n-1 ϕ i ϕ j , alongside with assumption A2. Noting σ2 := sup ψ∈Ψ E p ψ ϕϕ ⊤ -E p ψ ϕϕ ⊤ 2 F , which is finite by analyticity of the log-partition function, we have∥∆ 3,n ∥ F ≤ σα m C χ ∥ψ n-1 -ψ ⋆ ∥We now investigate the first moment. We have:E p ψ P m ψn-1 ϕ = (E ψ ⋆ P m ψn-1 ϕ -E p ψ n-1 ϕ) ∆4a,n +E p ψn ϕ =⇒ (E p ψ ⋆ P m ψn-1 ϕ)(E p ψ ⋆ P m ψn-1 ϕ) ⊤ -E p ψ n-1 ϕE p ψ n-1 ϕ ⊤ ∆4,n = ∆ 4a,n E p ψ n-1 ϕ ⊤ + E p ψ n-1 ϕ ⊤ ∆ 4a,nand thus, applying Lemma C.7 on ∆ 4a,n , we have:∥∆ 4,n ∥ F ≤ ∆ 4a,n E p ψ n-1 ϕ ⊤ + E p ψ n-1 ϕ ⊤ ∆ 4a,n ≤ ∆ 4a,n E p ψ n-1 ϕ ⊤ + E p ψ n-1 ϕ ⊤ ∆ 4a,n ≤ 2 ∥∆ 4a,n ∥ ∥E p ψ n-1 ϕ∥ ≤ 2 ∥log Z∥ 3,∞ α m C χ σ ∥ψ n-1 -ψ ⋆ ∥We can now combine our two matrix moment bounds to obtain: Consequently,E p ψ ⋆ E k m ψ n-1 ∆ 2,n ∆ ⊤ 2,n = E p ψn (ϕ -E p ψn ϕ)(ϕ -E p ψn ϕ) ⊤ + ∆ 3,n + ∆ 4,n Note that E p ψ n-1 (ϕ -E p ψ n-1 ϕ)(ϕ -E p ψ n-1 ϕ) ⊤ = Cov p ψ n-1 ϕ = ∇ 2 ψ L(ψ n-1 ). By analycity of L, there exists a constantM := sup ψ∈Ψ ∇ 3 log Z(ψ) op(∥•∥ F ,∥•∥ F ) such that ∥∇ 2 ψ L(ψ n-1 ) -∇ 2 ψ L(ψ ⋆ )∥ F := ∥∆ 5,n ∥ F ≤ M ∥ψ n -ψ ⋆ ∥ . Thus, E p ψ ⋆ E k m ψ n-1 ∆ 2,n ∆ ⊤ 2,n = I(ψ ⋆ ) + ∆ 3,n + ∆ 4,n + ∆ 5,n with ∥∆ 3,n + ∆ 4,n + ∆ 5,n ∥ F ≤ ∥∆ 3,n ∥ F + ∥∆ 4,n ∥ F + ∥∆ 5,n ∥ F ≤ (1 + M + ∥log Z∥ 3,∞ ) ∥ψ n-1 -ψ ⋆ ∥ .We are now ready to characterize the variance term alltogether. We have: Using the subadditivity of the square root ∆ 2,n sum, we have:1 n n i=1 tr(f ′′ (ψ ⋆ ) -1 E ψn-1 E p ψ ⋆ ∆ 2,i ∆ ⊤ 2,i f ′′ (ψ ⋆ ) -1 ) 1/2 ≤ 1 n n i=1 tr(I(ψ ⋆ )) -1 1/2 + 1 nµ n i=1 (1 + M + ∥log Z∥ 3,∞ )E ∥ψ n-1 -ψ ⋆ ∥ 1/2 ≤ tr(I(ψ ⋆ ) -1 ) 1/2 √ n + ∥log Z∥ 4,∞ µn n i=1 δ 1/2 i-1 1/2where ∥log Z∥ 4,∞ := (1 + M + ∥log Z∥ 3,∞ ) 1/2 . In these inequalities, we used the cyclicity of the trace, the fact that tr(A ⊤ B) ≤ ∥AB∥ F ≤ ∥A∥ op ∥B∥ F , the fact that ∥I(ψ ⋆ ) -1 ∥ op ≤ 1 µ , and the identityE ∥ψ i-1 -ψ ⋆ ∥ ≤ (E∥ψ i-1 -ψ ⋆ ∥ 2 ) 1/2 = δ 1/2 i-1 .We can combine this inequality with the rest of the term to obtain our final upper bound:E 1 n n i=1 ψ i -ψ ⋆ 2 1/2 ≤ 2tr(I(ψ ⋆ ) -1 ) √ n + ∥log Z∥ 4,∞ µn n i=1 δ 1/2 i 1/2 + E 1 n n i=1 g i (ψ i-1 ) 2 1/2 + M 2nµ n i=1 E ∥ψ i -ψ ⋆ ∥ 4 1/2We will use the readily available bounds of [22] for the last 2 terms. Before concluding the proof, we need to boundE ∥ψ i -ψ ⋆ ∥ 4 . Bounding E∥ψ i -ψ ⋆ ∥ 4 We use Lemma D.6 to obtain that M 2nµ n i=1 E ∥ψ i -ψ ⋆ ∥ 4 ≤ P 1 (n) with P 1 (n) of order n -β . Bounding 1/n n i=1 h i (ψ i-1 ) We use D.5 to obtain that (E 1 n n g i (ψ i-1 ) 2 ) 1/2 with P 2 (n) of order n β 2 -1 . Bounding ∥log Z∥ 4,∞ µn n i=1 δ 1/2 i 1/2 We use Lemma D.4 to obtain that ∥log Z∥ 4,∞ µn n i=1 δ 1/2 i 1/2 ≤ P 3 (n) with P 3 (n) of order n -1 2 + β 4 .Putting everything together, we have that:E ψ n -ψ ⋆ 2 1/2 ≤ 2tr(I(ψ) -1 ) √ n + P 1 (n) + P 2 (n) + P 3 (n) + P 4 (n) P 1 + P 2 + P 3 + P 4 is dominated by max -( 1 2 + β 4 ), -β, β 2 -1, β 2 + m |log α| log n = max -β, β 2 -1, β 2 + m |log α| log nwhich concludes the proof.In this section, we consider different gradient update schemes starting from some random initialization θ init , and control the L 2 distance between the different updates and the deterministic target ψ * ∈ Ψ.Notation. Recall the notation that X 1 , . . . , X n i.i.d.EQUATION∼ k m ψ (x, • ) Let θ init be some Ψ-valued random initialization that is possibly correlated with X 1 , . . . , X n . We capture the effect of correlation through the following quantities: For ϵ > 0 and ν > 2, letϑ init n,m (ϵ) := P n i=1 E ϕ K m θ init (X i ) X i , θ init -E ϕ K m θ init (X ′ i ) θ init n > ϵ , and ε init n,m;ν (ϵ) := ϵ 2 + κ 2 ν;m ϑ init n,m (ϵ) 2 ν-2 .We also consider the i.i.d. samples, drawn independently of X 1 , . . . , X n and on a given ψ ∈ Ψ, asX ψ 1 , . . . , X ψ n i.i.d. ∼ p ψ .For notational clarity, we shall use θ m,B to denote parameters arising from an one-step update, where the subscripts m, B represent performing the one-step update with length-m Markov chains and with batch size B. This is to be distinguished from ψ t elsewhere in the text, which denotes the parameter from the actual multi-step CD algorithm and the subscript t denotes the t-th CD iterate.Gradient update schemes. We consider five different updates. Let X ′ 1 be an i.i.d. copy of X 1 drawn independently of all other random variables. The SGD-with-replacement update is given byθ SGDw m,B := F SGDw m,B (θ init ) , where F SGDw m,B (ψ) := ψ - η B i∈S w ϕ(X i ) -ϕ K m i;ψ (X i )and S w is a uniformly drawn size-B subset of [n]. The SGD-without-replacement update, after renormalizing the learning rate, is given by the N -fold function compositionθ SGDo m,B := F SGDo m,B;N • . . . • F SGDo m,B;1 (θ init ) ,whereF SGDo m,B;j (ψ) := ψ - η N B i∈S o j ϕ(X i ) -ϕ K m i;ψ (X i ) for each j ∈ [N ] ,and S o j 's are disjoint size-B random subsets of [n], defined by (S o 1 , . . . , S o N ) = π([n]) for a uniformly drawn element π of the permutation group on n objects. The full-batch gradient update is given byθ GD m := F GD m (θ init ) , where F GD m (ψ) := ψ - η n i≤n ϕ(X i ) -ϕ K m i;ψ (X i ) .The full-batch gradient update with an infinite-length Markov chain is given byθ GD ∞ := F GD ∞ (θ init ) , where F GD ∞ (ψ) := ψ - η n i≤n ϕ(X i ) -ϕ(X ψ 1 ).The population gradient with an infinite-length Markov chain is given byθ pop := f pop (θ init ) , where f pop (ψ) := ψ -η E ϕ(X 1 ) -ϕ(X ψ 1 ) ,where we use the lowercase f to emphasize that f pop is a deterministic function.The forthcoming results are summarized below:θ SGDw m,B Lemma E.1 ≈ θ GD m Lemma E.2 ≈ θ GD ∞ Lemma E.3 ≈ θ pop Lemma E.4 ≈ ψ * .Lemma E.1. Let F n be the sigma algebra generated by{X i , K m i;θ init (X i ) | 1 ≤ i ≤ n}. Then E θ SGDw m,B -θ GD m θ init , F n = 0 almost surely .Moreover, under A1 and A7, we haveE θ SGDw m,B -θ GD m 2 ≤ 4η 2 (σ 2 + κ 2 ν;m ) B I {B<n} . Proof of Lemma E.1. Write A := (A 1 , . . . , A n ),whereA i := ϕ(X i ) -ϕ K m i;θ init (X i ) -E ϕ(X 1 ) -ϕ K m i;θ init (X 1 ) .Since S w is uniformly drawn from all size-B subsets of [n] and independently of all other variables, we have that almost surelyE θ SGDw m,B -θ GD m θ init , F n = E η B i∈S w A i - η n i≤n A i A = 0 .To prove the remaining bound, we note that the above relation implies θ SGDw m,B -θ GD m is zero-mean. By the law of total variance, we haveE θ SGDw m,B -θ GD m 2 = Tr Cov θ SGDw m,B -θ GD m = Tr E Cov θ SGDw m,B -θ GD m θ init , F n = η 2 Tr E E 1 B i∈S w A i 1 B i∈S w A i ⊤ A - 1 n i≤n A i 1 n i≤n A i ⊤ .To compute the covariance, recall that S w is a uniformly drawn size-B subset of [n] with B = n/N . Let P N ([n]) be the collection of all partitions of [n] into N size-B subsets. We can generate S w by the following two-step process:(i) Uniformly draw a partition P ′ = (P ′ 1 , . . . , P ′ N ) from P N ([n]); (ii) Uniformly sample an index K from [N ] and set S w = P ′ K . Then we have, almost surelyE 1 B i∈S w A i 1 B i∈S w A i ⊤ A - 1 n i≤n A i 1 n i≤n A i ⊤ = 1 |P N ([n])| P ′ ∈P N ([n]) 1 N k≤N 1 B 2 i,j∈P ′ k A i A ⊤ j - 1 n 2 i,j≤n A i A ⊤ j = 1 |P N ([n])| P ′ ∈P N ([n]) 1 N B 2 k≤N i,j∈P ′ k A i A ⊤ j - 1 N 2 B 2 k,l≤N i∈P ′ k ,j∈P ′ l A i A ⊤ j = 1 |P N ([n])| P ′ ∈P N ([n]) N -1 N 2 B 2 k≤N i,j∈P ′ k A i A ⊤ j - 1 N 2 B 2 k̸ =l i∈P ′ k ,j∈P ′ l A i A ⊤ j .By noting that A i 's are exchangeable, we obtainE θ SGDw m,B -θ GD m 2 = η 2 Tr E N -1 N B A 1 A ⊤ 1 + (N -1)(B -1) N B A 1 A ⊤ 2 - N -1 N A 1 A ⊤ 2 = η 2 Tr E N -1 N B A 1 A ⊤ 1 - N -1 N B A 1 A ⊤ 2 = η 2 (N -1) N B E∥A 1 ∥ 2 -E⟨A 1 , A 2 ⟩ (a) ≤ 2η 2 B E∥A 1 ∥ 2 = 2η 2 B E ϕ(X 1 ) -E[ϕ(X 1 )] -ϕ K m i;θ init (X 1 ) -E ϕ K m i;θ init (X 1 ) 2 ≤ 4η 2 B Tr Cov[ϕ(X 1 )] + Tr Cov ϕ K m i;θ init (X 1 ) (b) ≤ 4η 2 (σ 2 + κ 2 ν;m ) B .In (a), we have used a Cauchy-Schwarz inequality; in (b), we have used A1 and A7. Finally we note that if B = n, θ SGDw m,B = θ GD m almost surely, which implies the desired bound.Lemma E.2. Denote A n as the sigma algebra generated by θ init , X 1 , . . . , X n . Under A1, A2, A3 and A7, we have that for any ϵ > 0 and ν > 2,E E θ GD m -θ GD ∞ A n 2 ≤ η 2 α m σC χ E∥θ init -ψ * ∥ 2 + init n,m;ν (ϵ) 2 , E ∥θ GD m -θ GD ∞ ∥ 2 ≤ η 2 α m σC χ E∥θ init -ψ * ∥ 2 + ε init n,m;ν (ϵ) 2 + κ 2 ν;m + σ 2 n.Proof of Lemma E.2. The main challenge arises from the possible correlation between θ init and X 1 , . . . , X n . First note that for any ϵ > 0, ν > 2 and a real-valued random variable Y , by Hölder's inequality, we haveEQUATIONAlso note the useful inequality that for two real-valued random vectors (possibly correlated) V 1 , V 2 , we haveEQUATION)Now to control the first quantity of interest, by using a triangle inequality, we haveE E θ GD m -θ GD ∞ A n 2 = E E η n i≤n ϕ K m i;θ init (X i ) -ϕ X θ init i A n 2 (17) ≤ η 2 E[∆ 2 1 ] + E[∆ 2 2 ] 2 ,where∆ 1 := E 1 n i≤n ϕ K m i;θ init (X i ) -E ϕ K m i;θ init (X ′ 1 ) θ init A n = 1 n i≤n E ϕ K m i;θ init (X i ) θ init , X i -E ϕ K m i;θ init (X ′ 1 ) θ init , ∆ 2 := E 1 n i≤n E ϕ K m i;θ init (X ′ 1 ) θ init -ϕ X θ init i A n = E ϕ K m 1;θ init (X ′ 1 ) -ϕ X θ init 1 θ init ,and X ′ 1 is an i.i.d. copy of X 1 and in particular independent of θ init . ∆ 1 is controlled via ( 16):EQUATIONIn (a), we have plugged in the definition of ∆ 1 and applied a Jensen's inequality with respect to the empirical average; in (b), we have used A7 to bound the ν-th moment term asE ϕ(K m 1;θ init (X 1 )) -E ϕ(K m 1;θ init (X ′ 1 )) θ init ν 1/ν = E E ϕ(K m 1;θ init (X 1 )) -E ϕ(K m 1;θ init (X ′ 1 )) θ init ν θ init 1/ν ≤ sup ψ∈Ψ E ϕ(K m 1;ψ (X 1 )) -E ϕ(K m 1;ψ (X ′ 1 )) ν 1/ν = sup ψ∈Ψ E ϕ(K m 1;ψ (X 1 )) -E ϕ(K m 1;ψ (X 1 )) ν 1/ν ≤ κ ν;mand recalled the definitions of ϑ init n,m and ε init n,m;ν . On the other hand,E[∆ 2 2 ] = E R d ϕ(x)(K m 1;θ init p ψ * )(x)dx -E ϕ(X θ init 1 ) θ init 2 (a) = E R d (K m 1;θ init ϕ)(x) p ψ * (x)dx -E p θ init [ ϕ ] 2 (b) = E R d K m 1;θ init ϕ -E p θ init [ ϕ ] (x) × p ψ * (x)dx 2 (c) = E R d K m 1;θ init ϕ -E p θ init [ ϕ ] (x) × (p ψ * (x) -p θ init (x))dx 2 (d) = d l=1 E R d K m t1;θ init ϕ -E p θ init [ ϕ ] (x) ⊤ e l × p θ init (x) × p ψ * (x) -p θ init (x) p θ init (x) dx 2 (e) ≤ d l=1 E R d K m t1;θ init ϕ -E p θ init [ ϕ ] (x) ⊤ e l 2 p θ init (x)dx × R d p ψ * (x) -p θ init (x) p θ init (x) 2 p θ init (x)dx 2 (f ) ≤ d l=1 E α 2m Tr Cov ϕ(X θ init 1 ) θ init χ 2 (p ψ * , p θ init ) 2 (g) ≤ α 2m σ 2 C 2 χ E∥θ init -ψ * ∥ 2 .In (a), we have used that (Kf )(x) = K(x, y)f (y)dy; in (b), we have used that the Markov operator leaves the constant function invariant; in (c), we used that K 1;θ init leaves p θ init invariant; in (d), we denoted (e l ) l≤d as the standard basis vectors of R d and multiplied and divided by p θ init (x); in (e), we have used a Cauchy-Schwarz inequality; in (f ), we have used the definition of the spectral gap α in A3; in (g), we have used A1 and A2. Combining the bounds gives the first inequality thatE E θ GD m -θ GD ∞ A n 2 ≤ η 2 α m σC χ E∥θ init -ψ * ∥ 2 + ε init n,m;ν (ϵ) 2 .Now we handle the second quantity by conditioning on θ init and perform a bias-variance decomposition:E ∥θ GD m -θ GD ∞ ∥ 2 = η 2 E E 1 n i≤n ϕ K m i;θ init (X i ) -ϕ(X θ init i ) 2 A n = η 2 (Q B + Q V ) ,whereQ B := E E 1 n i≤n ϕ K m i;θ init (X i ) -ϕ(X θ init i ) A n 2 , Q V := E Tr Cov 1 n i≤n ϕ K m i;θ init (X i ) A n + Cov 1 n i≤n ϕ(X θ init i ) θ init .Note that the covariance terms separate because X θ init i is independent of K m i;θ init (X i ) conditioning on θ init . η 2 Q B is exactly the quantity controlled above, so it suffices to bound the variance term Q V . By explicitly computing the second covariance term while noting that X θ init 1 , . . . , X θ init n are conditionally i.i.d. given θ init and K m i;θ init (X i )'s are conditionally independent across 1 ≤ i ≤ n given A n , we haveQ V = i≤n E[Tr Cov[ϕ(K m i;θ init (X i )) | θ init , X i ]] n 2 + E Tr Cov[ϕ(X θ init 1 )|θ init ] n (a) = E[Tr Cov[ϕ(K m 1;θ init (X 1 )) | θ init , X 1 ]] n + E Tr Cov[ϕ(X θ init 1 )|θ init ] n ≤ κ 2 ν;m + σ 2 n, where we have used A4, A7 and A1 in the last line. Combining the bounds, we obtain thatE ∥θ GD m -θ GD ∞ ∥ 2 = η 2 (Q B + Q V ) ≤ η 2 α m σC χ E∥θ init -ψ * ∥ 2 + ε init n,m;ν (ϵ) 2 + κ 2 ν;m + σ 2 n . Lemma E.3. Under A1, E ∥θ GD ∞ -θ pop ∥ 2 ≤ 4η 2 σ 2 n .Proof of Lemma E.3. Since both F GD ∞ and f pop involve infinite-length Markov chains, the initializations do not matter and we can decouple the stochasticity of X i and K i;θ init . In particular,E ∥θ GD ∞ -θ pop ∥ 2 = η 2 E 1 n i≤n ϕ(X i ) -ϕ X θ init i -E ϕ(X 1 ) -ϕ X θ init 1 2 ≤ η 2 Q ′ 1 + 2 Q ′ 1 Q ′ 2 + Q ′ 2, whereQ ′ 1 := E 1 n i≤n ϕ(X i ) -E[ϕ(X 1 )] 2 = Tr Cov[ϕ(X 1 )] n = Tr ∇ 2 θ log Z(ψ * ) n ≤ σ 2 n , Q ′ 2 := E 1 n i≤n ϕ(X θ init i ) -E ϕ X θ init 1 2 = E Tr Cov[ϕ(X θ init 1 )|θ init ] n = E Tr ∇ 2 θ log Z(θ init ) n ≤ σ 2 n .In the computations above, we have used the relation ∇ 2 θ log Z(θ) = Cov X∼p θ [ϕ(X)] and the assumption sup θ∈Ψ tr(∇ 2 θ log Z(θ)) = σ 2 from A1. This implies the desired bound.Lemma E.4. Under A1, E ∥θ pop -ψ * ∥ 2 ≤ 1 -2µη + L 2 η 2 E ∥θ init -ψ * ∥ 2 .Proof of Lemma E.4. Recall thatf pop (θ ′ ) = θ -η E ϕ(X 1 ) -ϕ(X θ ′ 1 ) = θ -η ∇ ψ log Z(ψ * ) -∇ ψ log Z(θ ′ ). By construction, f pop is deterministic and f pop (ψ * ) = ψ * . By plugging in the recursions and expanding the square, we get thatE θ pop -ψ * 2 = E f pop (θ init ) -f pop (ψ * ) 2 = E (θ init -ψ * ) -η(∇ ψ log Z(θ init ) -∇ ψ log Z(ψ * )) 2 = E θ init -ψ * 2 -2η E θ init -ψ * , ∇ ψ log Z(θ init ) -∇ ψ log Z(ψ * ) + η 2 E ∇ ψ log Z(θ init ) -∇ ψ log Z(ψ * ) 2 ≤ E θ init -ψ * 2 -2µη E θ init -ψ * 2 + L 2 η 2 E θ init -ψ * 2 .In the last line, we have recalled inf ψ∈Ψ λ min (∇ 2 ψ log Z(ψ)) = µ and sup θ∈Ψ λ max (∇ 2 ψ log Z(ψ)) = L by A1 and applied Lemma C.8. Combining the coefficients gives the desired statement.We prove Theorem B.1 (which directly implies Theorem 4.3) and Theorem B.2 in this section. The key ingredient of both proofs is Lemma F.1 below, which provides an iterative error bound for the SGD-with-replacement scheme by combining different approximation bounds in Appendix E. Throughout this section, we denote δ SGDw t,j:= E ψ SGDw t,j -ψ * 2 .Lemma F.1. Under A1, A2, A3, A4 and A7, we have that for 1 ≤ j ≤ N -1, By adding and subtracting the auxiliary gradient updates followed by expanding the square, we obtainδ SGDw t,j ≤ 1 -η t µ -α m σC χ - L 2 2 η t δ SGDw t,δ SGDw t,j = E θ SGDw m,B -θ GD m + θ GD m -θ GD ∞ + θ GD ∞ -θ pop + θ pop -ψ * 2 (a) = E θ SGDw m,B -θ GD m 2 + E θ GD m -θ GD ∞ 2 + E θ GD ∞ -θ pop 2 + E θ pop -ψ * 2 + 2 E θ SGDw m,B -θ GD m , θ GD m -θ GD ∞ + 2 E θ SGDw m,B -θ GD m , θ GD ∞ -θ pop + 2 E E θ SGDw m,B -θ GD m θ init , F n , θ pop -ψ * + 2 E θ GD m -θ GD ∞ , θ GD ∞ -θ pop + 2 E E θ GD m -θ GD ∞ A n , θ pop -ψ * + 2 E θ GD ∞ -θ pop , θ pop -ψ * (b) ≤ 4η 2 t (σ 2 + κ 2 ν;m ) √ B I {B<n} + η 2 t α m σC χ δ SGDw t,j-1 + ε SGDw n,m,t;ν (ϵ) 2 + κ 2 ν;m + σ 2 n + 4η 2 t σ 2 n + (1 -2µη t + L 2 η 2 t )δ SGDw t,j-1 + 2 2ηt σ 2 + κ 2 ν;m √ B I {B<n} η t α m σC χ δ SGDw t,j-1 + ε SGDw n,m,t;ν (ϵ) + κ 2 ν;m + σ 2 √ n + 2 2ηt σ 2 + κ 2 ν;m √ B I {B<n} 2ηtσ √ n + 0 + 2η t α m σC χ δ SGDw t,j-1 + ε SGDw n,m,t;ν (ϵ) + κ 2 ν;m + σ 2 √ n 2ηtσ √ n + 2η t α m σC χ δ SGDw t,j-1 + ε SGDw n,m,t;ν (ϵ) (1 -2µη t + L 2 η 2 t )δ SGDw t,j-1 + 2 2ηtσ √ n (1 -2µη t + L 2 η 2 t )δ SGDw t,j-1 (c) ≤ 1 -2µη t + Lη 2 t + η 2 t α 2m σ 2 C 2 χ + 2η t α m σC χ 1 -2µη t + L 2 η 2 t × δ SGDw t,j-1 + 2 η 2 t α m σC χ ε SGDw n,m,t;ν (ϵ) + 2η 2 t ( σ 2 + κ 2 ν;m + σ) α m σCχ √ B + η t ε SGDw n,m,t;ν (ϵ) 1 -2µη t + L 2 η 2 t + 2ηtσ 1 -2µηt + L 2 η 2 t √ n × δ SGDw t,j-1 + η 2 t 4(σ 2 + κ 2 ν;m ) B + ε SGDw n,m,t;ν (ϵ) 2 + 5κ 2 ν;m + 9σ 2 B + 4 σ 2 + κ 2 ν;m ε SGDw n,m,t;ν (ϵ) √ B + 8σ σ 2 + κ 2 ν;m B + 4σε SGDw n,m,t;ν (ϵ) √ B =: (A 1 ) × δ SGDw t,j-1 + (A 2 ) × δ SGDw t,j-1 + (A 3 ) .In (a), we have expanded the square, used F n defined in Lemma E. EQUATIONc), we have grouped the terms by powers of δ SGDw t,j-1 , bounded the indicator function from above by 1 and used B ≤ n to replace n in the denominator by B. While the computation above is complicated, we remark that the key steps are taking the conditional expectations for the crossterms in (b), which gives us a tighter bound than directly applying a triangle inequality of the formE∥Y 1 + Y 2 ∥ 2 ≤ ( E∥Y 1 ∥ 2 + E∥Y 2 ∥ 2 ).To further simplify the bounds, we seek to bound each coefficient by a square, which yields(A 1 ) = η t α m σC χ + 1 -2µη t + Lη 2 t 2 , (A 2 ) ≤ 2 η t α m σC χ + 1 -2µη t + Lη 2 t × η t ε SGDw n,m,t;ν (ϵ) + 4σ + 2κν;m √ B , (A 3 ) ≤ η 2 t ε SGDw n,m,t;ν (ϵ) 2 + 8σ 2 + 4κ 2 ν;m √ B ε SGDw n,m,t;ν (ϵ) + 21σ 2 + 17κ 2 ν;m B ≤ η 2 t ε SGDw n,m,t;ν (ϵ) + 5σ + 5κν;m √ B 2 .This impliesδ SGDw t,j ≤ (A 1 ) × δ SGDw t,j-1 + (A 2 ) × δ SGDw t,j-1 + (A 3 ) ≤ η t α m σC χ + 1 -2µη t + Lη 2 t δ SGDw t,j-1 + η t ε SGDw n,m,t;ν (ϵ) + 5σ + 5κν;m √ B 2 . Now note that since µ ≤ L by the definitions in A1, 2µη t -L 2 η 2 t ≤ 2Lη t -L 2 η 2 t ≤ 1. By using √ 1 -x ≤ 1 -x 2 for all x ≤ 1, we get that 1 -2µη t + L 2 η 2 t ≤ 1 -µη t + L 2 2 η 2 t .Substituting this into the earlier bound and taking a square-root, we obtain the desired bound thatδ SGDw t,j ≤ 1 -η t µ -α m σC χ - L 2 2 η t δ SGDw t,j-1 + η t ε SGDw n,m,t;ν (ϵ) + 5σ + 5κν;m √ B .Moreover, since L ≥ µ by definition from A1, we have1 -η t µ -α m σC χ - L 2 2 η t = L √ 2 η t -1 2 + ( √ 2 L -µ)η t + α m σC χ η t > 0 ,which finishes the proof.Under A1, A2, A3, A4 and A7, Lemma F.1 impliesδ SGDw t,j ≤ 1 -μm Ct -β + L 2 C 2 2 t -2β δ SGDw t,j-1 + Ct -β σ SGDw n,Tfor 1 ≤ t ≤ T and 1 ≤ j ≤ N , where we have used μm = µ -α m σC χ , η t = Ct -β andε SGDw n,m,t;ν (ϵ) + 5σ + 5κν;m √ B ≤ ε SGDw n,m,T ;ν (ϵ) + 5σ + 5κν;m √ B = σ SGDw n,T.By an induction on j = 1, . . . , N , we haveδ SGDw t,N ≤ 1 -μm Ct -β + L 2 C 2 2 t -2β N δ SGDw t,0 + Cσ SGDw n,T t -β N j=1 1 -μm Ct -β + L 2 C 2 2 t -2β j-1 .By noting that δ SGDw t,0 = δ SGDw t-1,N almost surely for t ≥ 2 and using another induction on t = 1, . . . , T ,EQUATIONwhereQ t := T s=t+1 1 -μm Cs -β + L 2 C 2 2 s -2β N and A t := N j=1 1 -μm Ct -β + L 2 C 2 2 t -2β j-1 .Note that by Lemma F.1, we also have that 1EQUATIONt -2β > 0. This allows us to apply Lemma C.3 to obtainκ 0 ≤ exp 1 -N μm Cφ 1-β (T + 1) + N L 2 C 2 2 φ 1-2β (T + 1) = E T,N 1 . To control T t=1 Q t A t t -β , recall that β ∈ [0, 1], μm > 0, L 2 C 2 2> 0, and thatE T,N 2 = exp - N μmC 2 φ 1-β (T + 1) + 2N L 2 C 2 φ 1-2β (T + 1) .We can now apply Lemma C.3: If β ̸ ∈ { 1 2 , 1}, then T t=1 Q t A t t -β ≤ 2 2β+1 μmC e μmC 2(1-β) N (T +1) β + 3 β (1 + μmC) N -1 (T + 2) β L 2 C 2 E T,N 2 . If β = 1 2 , i.e. 2β = 1, we have T t=1 Q t A t t -β ≤ 4 μmC e μmCN (T +1) 1/2 + 2N (1 + μm C) N -1 φ 1 2 -L 2 C 2 N (T + 1) E T,N 2 . If β = 1, we get that T t=1 Q t A t t -β ≤ 4 μmC + 3N 1 + L 2 C 2 2 N -1 e 2L 2 C 2 N log(T + 1) (T + 1) ( μmCN )/2Substituting the bounds into (19) , we get the desired bound that δ SGDw T,N is upper bounded byEQUATIONF. , we condition on S o t , which in particular fixes S o t,j , the last size-B subset of [n] chosen. We then identify θ init = ψ SGDo t,j-1 , η = η t and the dataset used as D o t,j := (X i : i ∈ S o t,j ), which allows us to identify ψ SGDo t,j = θ GD m , the full-batch gradient descent update using D o t,j . Meanwhile, we note that almost surely θ GD m = θ SGDw m,B , the SGD-with-replacement iterate that uses the full dataset D o t,j . Observe that the proof of Lemma F.1 holds with δ SGDo t,j-1 replaced by any random initialization θ init possibly correlated with X 1 , . . . , X n , which allows us to obtainδ SGDo t,j ≤ 1 -η t µ -α m σC χ - L 2 2 η t δ SGDo t,j-1 + η t ε SGDw n,m,t;ν (ϵ) + 5σ + 5κν;m √ B .Since the error recursion for δ SGDo We present the proofs for results in Appendix B.3 that control the tail probability terms ϑ SGDw ν;n,m,T and ε SGDw ν;n,m,T .Proof of Lemma B.3. For δ > 0, let N δ be the δ-covering number of Ψ, which satisfies N δ ≤ r Ψ (1 + 2/δ) p (Example 5.8, [50] ). Note also that by the Jensen's inequality applied to E[ • |X 1 ]and Assumption A6, there exist some σ m , ζ m > 0 such that, for any z∈ R p with ∥z∥ ≤ ζ m , E[e z ⊤ (E[ϕ(K m ψ * (X1))|X1]-E[ϕ(K m ψ * (X1))]) ] ≤ E[e z ⊤ (ϕ(K m ψ * (X1))-E[ϕ(K m ψ * (X1))]) ] ≤ e σ 2 m ∥z∥ 2 /2 .Meanwhile under A5, recycling the proof of Lemma 3.1 of [21] shows that ifC m δ/ √ p < σ 2 m ζ m , P sup ψ∈Ψ 1 n n i=1 E ϕ K m ψ (X i ) X i -E ϕ(K m ψ (X i )) > 3C m δ ≤ 2N δ p exp - nC 2 m δ 2 2pσ 2 m ≤ 2(r Ψ ) p exp p log(1 + 2δ -1 ) - nC 2 m δ 2 2pσ 2 m .Since the probability above is an upper bound to ϑ SGDw n,m,T (3C m δ), we get that if Proof of Lemma B.4. Let δ > 0 and N δ be defined as in the proof of Lemma B.3, with N δ ≤ r Ψ (1 + 2/δ) p . Let (ψ l ) N δ l=1 be the centers of the covering δ-balls. The covering-ball argument of the proof of Lemma 3.1 of [21] shows thatC m δ/ √ p < σ 2 m ζ m , ε SGDw ν;n,m,T (3C m δ) 2 = 9C 2 m δ 2 + κ 2 ν;m ϑ SGDw n,m,T Cmδ √ p ν-2 ν ≤ 9C 2 m δ 2 + κ 2 ν;m 2 2(ν-2) ν (r Ψ ) (ν-2)p ν exp (ν -2)p ν log(1 + 2δ -1 ) - n(ν -2)C 2 m δ 2 2νpσ 2 m . Recall that by assumption, log n n < σ 2 m ζ 2 m p+ν-2 . We now choose δ = σm Cm p p + 2ν ν -2 × log n n = σm Cm p × (ν -2)p + 2ν ν -2 × log n n , which implies inf ϵ>0 ε SGDw ν;n,m,T (ϵ) 2 ≤ ε SGDw ν;n,m,T (3C m δ) 2 ≤ 9σ 2 m p((ν -2)p + 2ν) log n (ν -2)n + κ 2 ν;m 2 ν-2 ν (r Ψ ) (ν-2)p ν (1 + 2δ -1 ) (ν-2)p ν exp - (ν -2)p + 2ν 2ν log n ≤ 9σ 2 m p((ν -2)p + 2ν) log n (ν -2)n + κ 2 ν;m 2 ν-2 ν (r Ψ ) (ν-2)p ν 1 + 2Cm(ν -2) 1/2 σmp 1/2 ((ν -2)p + 2ν) 1/2 √ n √ log n (ν-2)p ν n -(ν-2)p+2ν 2ν ≤ 9σ 2 m p((ν -2)p + 2ν) log n (ν -2)n + κ 2 ν;m 2 ν-2 ν (r Ψ ) (ν-2)p ν 1 + 2Cm(ν -2) 1/2 σmp 1/2 ((ν -2)p + 2ν) 1/2 (ν-2)p ν n (ν-2)p 2ν - (ν-2)p+2ν 2ν ≤ 9σ 2 m p((ν -2)p + 2ν) ν -2 + κ 2 ν;m 2 (ν-2) ν (r Ψ ) (ν-2)p ν 1 + 2Cm(ν -2) 1/2 σmp 1/2 ((ν -2)p + 2ν) 1/2ϑ SGDw n,m,T (3C m δ) ≤ P sup ψ∈Ψ 1 n n i=1 E ϕ K m ψ (X i ) X i -E ϕ(K m ψ (X i )) > 3C m δ ≤ N δ l=1 P 1 n n i=1 E ϕ K m ψ l (X i ) X i -E ϕ(K m ψ l (X i )) ≥ C m δ .By a Markov's inequality followed by the Burkholder's inequality applied to an average of i.i.d. summands (see e.g. [51] for ν > 2), there exists a constant C ν > 0 depending only on ν such that ϑ SGDw n,m,T (3C m δ) ≤ N δ l=1 E n i=1 E ϕ K m ψ l (X i ) X i -E ϕ(K m ψ l (X i )) ν n ν C ν m δ ν ≤ N δ l=1 E E ϕ K m ψ l (X 1 ) X i -E ϕ(K m ψ l (X 1 )) ν n ν/2 C ν m δ ν• 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: [NA] Justification: The paper does not include experiments.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?Justification: As a theory paper, we believe that our work does not breach the code of conduct. With this work, we aim to advance the understanding of statistically efficient algorithms, a domain which has the potential to improve the overall computational efficiency of machine learning algorithms.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: [NA] Justification: We believe that our work will not have the negative impacts mentioned (e.g., disinformation, surveillance, fairness, privacy, security considerations). We believe that advancing the domain of statistical effiency has the potential to improve the overall computational efficiency of machine learning algorithms, which, all else left equal, could constitute a positive societal impact.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Since our work does not include experiments nor data, and only analyzes existing models, we do not believe that our work requires any specific guidelines.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: This paper studies the maximum-likelihood estimator for an exponential family. It provides guarantees for the Contrastive Divergence algorithm, which uses stochastic gradient descent on the cross-entropy loss between the model estimate and empirical data.
The reviewers found the presentation clear and the contribution valuable for the NeurIPS community. I recommend accepting this paper

Reviews: review: paper_summary: paper_summary: This paper studies the problem of finding a maximum-likelihood estimator for an exponential family. The paper proves guarantees for the Contrastive Divergence (CD) algorithm, which iteratively performs stochastic gradient descent on the cross-entropy loss between the estimated model and empirical data distribution. However, each step of gradient descent is performed by estimating the gradient with an MCMC algorithm.

The paper proves:
* Thm 3.2: Online CD (i.e., using a fresh sample on each step), with enough MCMC steps, leads to $O(n^{-1/2})$ error rate with $n$ samples
* Thm 3.3: Online CD with Polyak-Ruppert averaging achieves (up to a constant factor of <= 2), an asymptotically optimal rate.
* Thm 4.2: Offline CD (i.e., using a batch of samples at each step) leads to $O((\log(n) / n)^{1/2} + 1/\sqrt{B})$ error rate with batch size B.
strength_weakness: strength_weakness: My main concern is on whether the subject of this paper will be of sufficient interest to a large enough portion of the NeurIPS audience, or whether a different venue/journal with a more statistical audience might be more appropriate. The 
results of this paper are technical bounds on how CD converges for maximum-likelihood estimation in exponential families; and the CD algorithm is not (currently) a focus of the machine learning community. Furthermore, the rate improvement proved is relatively modest -- from the previously known $O(n^{-1/3})$ rate in Jiang et al. [21] to the current best possible $O(n^{-1/2})$ rate, so to some extent the qualitative nature of the result was previously known.

Very minor:
* Assumption A.4, which is needed for Theorem 4.2 is stronger than what Jiang et al. [21] assumes (see Remark 4.1). But Theorem 4.2 is stronger than what is proved in Jiang et al. [21].
* The bound in Line 23 of Theorem 3.2 is quite elaborate and difficult to check/use. However, the most important case of $\\beta = 1$ is conveniently solved.
* On line 95, the idiom “silver lining” is used incorrectly. It should indicate something positive, rather than something negative.
questions: questions: * What is $\\tilde{\\mu}$ in Theorem 3.2? Do you mean $\\tilde{\\mu}_m$?
limitations: limitations: Yes
paper_summary: paper_summary: The paper delivers a detailed non-asymptotic analysis of the Contrastive Divergence (CD) algorithm applied to unnormalized exponential family distributions. It significantly advances the understanding of CD by showing that, under certain regularity conditions, the algorithm can achieve the parametric convergence rate of \(O(n^{-1/2})\). This finding represents a notable improvement over the previously established convergence rate of \(O(n^{-1/3})\). The paper's analysis is rooted in a meticulous exploration of how CD's bias and variance components affect its convergence behavior, providing deeper insights into the algorithm's performance under various conditions.
strength_weakness: strength_weakness: The paper focuses heavily on theoretical analysis with less emphasis on empirical validation. Including comprehensive experiments that benchmark CD against other estimation techniques across various real-world datasets could strengthen the claims.

The paper focuses heavily on theoretical analysis with less emphasis on empirical validation. Including comprehensive experiments that benchmark CD against other estimation techniques across various real-world datasets could strengthen the claims.
questions: questions: Could you provide more detailed explanations or examples of the "specific regularity conditions" under which CD achieves the improved convergence rate of $O(n^{-1/2})$?

Are there plans to complement the theoretical findings with empirical validations?
limitations: limitations: The authors listed the limitations and future directions in their discussion section
paper_summary: paper_summary: The paper studies parameter estimation problem in naturally-parametrized exponential families using the contrastive divergence methods, i.e., SGD where the stochastic gradients are estimated via MCMC. Both online and offline variants of the algorithms are considered. Building upon existing works on stochastic approximation, the paper establishes upper bounds on the last iterate and averaged iterate for online algorithms. In the offline setting, the paper establishes a near-parametric rate under certain tail assumptions.
strength_weakness: strength_weakness: On the technical side, the results are mostly applying existing techniques directly. The online results are mostly from Bach and Moulines (2011), with a careful analysis of the biases.

I'm not totally convinced by the usefulness of the offline algorithm. It seems that the analysis is much more involved and requries additional technical assumptions. On the other hand, the results are weaker than the online setting.
questions: questions: - Unlike Bach and Moulines (2011), the claim of matching Cram\'{e}r--Rao lower bound in Theorem 3.3 involves an additional factor of 2. I think this is unavoidable because of the additional variance from the Monte Carlo sample. This needs to be discussed more thoroughly.
- The writing of Section 4 needs some improvement. The notations are extremely heavy, making the main message difficult to follow. I would suggest simplifying the expression under optimal choice of parameters, and present the main result clearly. The heavy notations can be introduced in the proof sections.
Minor comments:
- Eq (5). Definition of $\mu$ is incorrect. Strong convexity parameter should be the minimum eigenvalue, instead of the operator norm.
limitations: limitations: The paper addressed limitations of the results adequately.
paper_summary: paper_summary: This paper analyzes the contrastive divergence algorithm for learning exponential family distributions.
Both online and offline versions are studied. The main contribution of this paper is to establish the parametric convergence rate $O(n^{-1/2})$ of the estimator, which improves the previously best result of $O(n^{-1/3})$.
strength_weakness: strength_weakness: - That said, the paper is extremely technical, loaded with heavy notation. Though I am aware of some related works on learning unnormalized models, this paper is hard to appreciate. It is hard to penetrate without extreme familiarity with the previous analysis in the literature and all the optimization analysis tools. I guess this paper could be certainly of some audience, but it'd be great if the authors can provide some more high-level explanations, pointing out some difficulties in the analysis, required techniques, and so on.
- In particular, assumptions are not with concrete examples. Can you provide a few examples of exponential family distributions that satisfy and don't satisfy each assumption?
questions: questions: - Is it proper to say "near-optimality", if this refers to the parametric rate $O(n^{-1/2})$, and it indeed achieves CRLB up to a factor of 4? It's still not efficient due to the multiplicative factor. Is this dependence usually called "nearly optimal"? For example, compared to the NCE estimators, can you argue that what's guaranteed with the CD estimator is provably better in a certain regime? This kind of insights would be very useful.
limitations: limitations: There is no specific limitation noted in the paper.
It'd be very helpful if the authors can provide a comparison of assumptions and attained rates from previous papers to the ones in the paper, side-by-side.
paper_summary: paper_summary: The paper provides a non-asymptotic analysis of the Contrastive Divergence (CD) algorithms, demonstrating their ability to achieve parametric convergence rates and near-optimal asymptotic variance in the online setting. It also shows near-parametric rates in the offline setting, extending previous results. The analysis is applicable to unnormalized exponential families, which, despite their flexibility, cover only log-densities with linear relationships on the model parameters. 
The analysis also covers various data batching schemes, including fully online and minibatch methods.
The study highlights the importance of extending these results to more general models in future work. Additionally, the assumption of $\chi^2$-regularity of the exponential family model is crucial for mapping Markov contractions back to the parameter space. Future research could consider alternative conditions on the Markov Chain Monte Carlo (MCMC) chains, such as geometric ergodicity, to avoid the presence of potentially large regularity constants.
strength_weakness: strength_weakness: While I acknowledge that I am not an expert in this specific area, I have some suggestions that could improve the paper towards its stated goals. Although I did not thoroughly check the proofs due to the extensive 30-page appendix, here are my detailed comments:

**Experimental Validation**: While the primary focus of this paper is theoretical, it would be beneficial to include an intuitive demonstration of the $O(n^{-1/2})$ rate, if possible. This could help bridge the gap between theory and practical understanding, making the results more accessible and relatable to a broader audience.

**Clarity and Structure**: The paper currently presents a substantial amount of theoretical content, which can be overwhelming, especially for readers who may not be deeply familiar with the notation and concepts used. Highlighting the new techniques that significantly improve upon previous work would make the contributions clearer.
Consider focusing on one main theorem with detailed lemmas and techniques that build towards it. This approach could provide a more structured and digestible presentation, making it easier for readers to follow the logical progression of the work.

**Assumptions and Practicality**: The paper requires a slightly stronger moment assumption (Assumption A4) to achieve better results. Providing more detailed examples and explanations that demonstrate the practical acceptability of this assumption would be valuable. This could help convince readers that the stronger assumption is reasonable and applicable in real-world scenarios.

I welcome further discussion on these points to ensure a comprehensive evaluation. Thanks!
questions: questions: See the weaknesses part.
limitations: limitations: Yes, the authors adequately addressed the limitations in conclusion part. The paper is theoretical, and I believe there is no need to confirm societal impact.

Comments: title: Response to reviewer's questions
Response
Response
Discussion of proofs, comparision with score matching
comment: Thank you for the response. I have read the rebuttal and decide to keep my score.
We thank the reviewer for their response. We answer their questions below.

**Regarding our comparison between online and offline CD, and its fairness**:

- First, the bounds provided in our paper are **the only currently existing nonasymptotic bounds for the CD algorithm** (and the contribution of our work). We do not believe that comparing these two bounds (and stating when offline CD is likely to be more efficient than online CD) as done in the paper and in the rebuttal is unfair to online CD. Note that we take into account all step size schedules of the form $Ct^{-\beta}$ in our comparison.
- Second, while adapting the guarantees of **other** online unbiased SGD algorithms to the CD setting is an interesting avenue for future work (and we are happy to mention that in the paper), the bounds are **not** expected to transfer without fundamental modifications from unbiased SGD to CD: for instance, the tail-averaging algorithm of [58] (Equation 6a) contains averaging **in its update rule**, which may introduce CD-specific correlations that are precisely the ones requiring a tail control as done in Section 4, and that lead to the additional $\log (n)$ factor in the rate. The rate of such variants may thus be **worse** than (averaged) online CD. Offline CD, as we show, provides a simple optimization algorithm with an optimal (up to a log factor) convergence rate (sharply improving over [21]) and that is robust to multiple hyper-parameters of the problem, unlike averaged online CD. Trying to extend the guarantees of other online SGD algorithms to the CD setting would be interesting for future work, but the absence of such results in our paper does not take away the value of our analysis of offline CD, and its benefits compared to online CD.

**Regarding the rates, leading order terms, and high-order terms of offline CD**: all such quantities are given in the current manuscript. The rates (depending on the setting) are given in Lemma 4.3 and Lemma 4.4 of the main body. The leading order terms are made explicit (to keep the main body simple) in Lemma B.3 and B.4 of the appendix. The higher order terms are given in Theorem B.2 of the appendix, but importantly vanish for fixed $n$ as $T \to \infty$. The rate is optimal up to a $\log(n)$ factor, while the multiplicative constant may be improvable. We will mention these points and will make clearer references to these results in the main body.

We thank the reviewer for engaging with us and our manuscript. We are happy to answer any remaining questions they may have.

[58] ROOT-SGD: Sharp Nonasymptotics and Asymptotic Efficiency in a Single Algorithm, Li et. Al
Thanks for addressing my concerns. I appreciate the authors's efforts on improving the presentation of Section 4. I have the following questions:
- The drawbacks of online CD (dependence on initial distance, exponential high-order terms) can be overcome by more careful choice of stepsize and tail averaging techniques. These are not discussed and the comparison to offline CD is unfair.
- For offline CD, it would be helpful to make it clear what the rate, leading-order term, and high-order terms are, and whether they're optimal.
I think there could be a publishable paper here, but the current version looks unclear. So I'll keep my score.
Thank you for addressing my concerns. I am still hesitant about NeurIPS being a good audience, but I will raise my score since I would not be upset if this paper gets accepted.
I really appreciate the authors' detailed response.
As most of the reviewers raised concerns of its extreme technicality, the paper will be better appreciated by a larger audience, after a through revision as explained and promised by the authors. 
I'd like to say that the added simulation is also nice, in that it confirms the theory by a simple example.
I believe that the authors will do the best in revising the manuscript.
I will keep the score as is.
Thank you for the response. I maintain my current evaluation of the paper, as it is highly technical and covers areas I am unfamiliar with. I believe the paper could be improved based on the current modifications shown in the rebuttal. Good luck to you.
In the following comment, we compare our bounds with existing ones for score matching, and provide a high level discussion of the proofs of the theoretical results in the manuscript.

**Regarding score matching** The asymptotic variance of score-matching was analyzed in [8], where it is shown that:
- Under a Poincaré Inequality, the asymptotic variance of score matching scales as the square of CRLB (times a smoothness term) (See [8, Theorem 2]).
- For some specific exponential families, the asymptotic variance of score matching can scale exponentially badly with the parameters of the problem (see [8, Corollary 1]).


**Regarding high-level explanation of the proofs** As suggested by the reviewer, to further improve the accessibility of the paper, we provide below a high-level discussion of each result of the main body, and the main proof ingredients used. The overall approach leading the online and offline bounds consists in using non-asymptotic arguments from the convex optimization literature, combined with a control of the bias (and the correlations in the offline setting) introduced by using the CD gradients. While the online setting was not studied in prior analyses of CD, the offline setting was studied in [21]. Our approach significantly departs from [21], which obtained its results by considering constant step sizes $\eta_t = \eta$, thanks to which the sequence of iterates $\psi_t$ becomes a Markov Chain. [21] then obtain convergence by showing the recurrence of this Markov chain in a ball containing the MLE. We now discuss the online and the offline results in more detail.

**High level discussion of the online result proofs** In the online setting, the fact that every sample $x_t$ is only used once makes $\psi_{t}$ and $x_{t}$ conditionally independent given past information. Thanks to this conditional independence, the final recursion is free of correlation terms, and we only need to control the bias and variance induced by using CD instead of standard SGD. The lack of correlations allows to obtain an optimal rate of convergence. A fine control of the bias and the variance, as done in Theorem 3.3, allows to show the that the averaged CD iterates are near-optimal estimators. This correlation-free setting also allows to get rid of multiple assumptions compared the offline analysis of [21].

- Lemma 3.1 is a recursion on the iterates describing the distance of of the estimated parameter to the true parameter during optimization. This recursion is obtained by combining standard arguments in the convex optimization litterature (See [22, Equation 16]) with a control of the approximation error due to using the CD gradient instead of the true likelihood gradient. This control uses arguments similar in spirit to [21, Lemma 3.2] although with considerable simplification thanks to the online setup.
- Theorem 3.2 provides a bound on the distance of the n-th iterate of the previous sequence, by "unrolling" the recursion. The form of the recursion allows to directly invoke [20] to unroll the recursion, allowing for a very short argument.
- Theorem 3.3 provides a bound on $\mathbb E ||1/n \sum_{t=1}^{n} \psi_t - \psi^\star ||^2$, the distance of the averaged CD iterates to the true parameter. The proof consists in adapting the argument of [20, Theorem 3], with multiple modifications to account for the fact that 
  - The CD gradient is an unbiased estimate of a *biased* gradient. Consequently, a different decomposition to the one in the proof of [20, Theorem 3] must be used (l. 1133 to l. 1139 of our manuscript)
  - The CD gradient contains an additional source of noise coming from from using one MCMC sample to approximate $\mathbb E_{p_{\psi^\star}} \phi$. However (and this is key to our final result), we relate this additional noise to the Fisher Information of the model at the true parameter, which is why the CRLB is multiplied by factor of 4.
  - While averaging the CD iterates yields an estimator with smaller variance, it induces the accumulation of a bias term, which we can control by scaling the number of MCMC steps with $\log n$.


**High-level discussion of the offline results proof** In the offline setting, already studied in [21], the conditional independence structure of the online setting disappears. Instead, the iterates $\psi_t$ are correlated with the data points $x_t$, and these correlations need to be controlled to obtain our final bounds. From a technical point of view, the control of these correlations constitutes the main difference between the online and the offline setting (see Lemma E.2), and requires arguments that are specific to our setting, and not standard in the SGD literature. The rest of the arguments used to obtain our bounds (Lemma E.1, E.3 and E.4) are either adaptations from classical results in the SGD literature, or similar in spirit to the one developed in the online section (with adaptations to the offline setting).