Fast Doubly-Adaptive MCMC to Estimate the Gibbs Partition Function with Weak Mixing Time Bounds
Year: 2022

Abstract:
We present a novel method for reducing the computational complexity of rigorously estimating the partition functions (normalizing constants) of Gibbs (Boltzmann) distributions, which arise ubiquitously in probabilistic graphical models. A major obstacle to practical applications of Gibbs distributions is the need to estimate their partition functions. The state of the art in addressing this problem is multi-stage algorithms, which consist of a cooling schedule, and a mean estimator in each step of the schedule. While the cooling schedule in these algorithms is adaptive, the mean estimation computations use MCMC as a black-box to draw approximate samples. We develop a doubly adaptive approach, combining the adaptive cooling schedule with an adaptive MCMC mean estimator, whose number of Markov chain steps adapts dynamically to the underlying chain. Through rigorous theoretical analysis, we prove that our method outperforms the state of the art algorithms in several factors: (1) The computational complexity of our method is smaller; (2) Our method is less sensitive to loose bounds on mixing times, an inherent component in these algorithms; and (3) The improvement obtained by our method is particularly significant in the most challenging regime of high-precision estimation. We demonstrate the advantage of our method in experiments run on classic factor graphs, such as voting models and Ising models.

heading: 1 Introduction
text: The Gibbs (Boltzmann) distribution is a family of probability distributions of exponential form. First introduced in the context of statistical mechanics [25], Gibbs distributions are now ubiquitous in a variety of other disciplines, such as chemistry [24, 31], economics [1, 54] and machine learning. Gibbs distributions are typically used to model the global state of a system as a function of a collection of interdependent random variables, each representing local states in the system. The dependencies ∗ r© indicates randomized ordering and equal contribution.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
in the system are modeled by a Hamiltonian function, and the probability distribution is inversely proportional to exponent of the Hamiltonian scaled by the temperature (see eq. (1) § 1.1).
Gibbs distributions provide potent statistical inference tools in many machine learning applications. They appear in probabilistic graphical models [41, 49, 51], including restricted Boltzmann machines [43, 66], Markov random fields [40, 47], and Bayes networks [32], and are applied in the analysis of images and graphical data [21, 23, 44, 65], topic modeling (LDA) [27, 53, 57, 62], and more [2, 13, 18, 19, 26, 30, 48, 56, 58, 68].
A major obstacle in applying the Gibbs distribution in practice is the need to compute, or estimate, its partition function (normalizing constant), henceforth written GPF. The partition function is defined over the Cartesian product of supports of a (typically large) number of variables, making exact computation intractable. Furthermore, due to interdependence of variables, exact sampling is not practically feasible, thus Markov-chain Monte-Carlo (MCMC) solutions for this problem have been extensively studied [6, 22, 29, 34, 37, 39, 42, 52, 63, 67].
Like other MCMC methods, here various heuristics are used. The most well-known heuristics are the annealed importance sampling [39, 52, 64] or convergence diagnostics methods [10, 11, 20, 61]. Unfortunately, these methods are often error-prone, as their correctness is only proven asymptotically, without rigorous mathematical analysis to bound their estimation error with finite samples. In fact, theoretical findings have shown that with no prior knowledge of relevant measures, such as the variance of importance weights in annealed importance sampling, or upper bounds on mixing or relaxation times for diagnostic methods, these methods are either unreliable or computationally intractable (see [52, section 4] or [8, 33]).
On the other hand, theoreticians study this problem by designing Fully Polynomial Randomized Approximation Schemes (FPRAS) (see problem 1). The state of the art FPRAS for estimating the GPF is a multi-stage algorithm involving a sequence of functions at various temperatures, such that the expectation of the product of these functions, or the product of the expectations of said functions, is the GPF. FPRAS’s are proven to produce (approximate) solutions w.h.p., but their performance guarantees rely on available upper-bounds on various measures such as variances of estimators or mixing times of Markov chains. In static algorithms, these upper-bounds are given a priori, and adaptive2 algorithms estimate them dynamically, while increasing the sample size until desired properties are mathematically guaranteed. Thus, adaptive algorithms are less sensitive to looseness of known upper-bounds, more robust, often faster, and more easily applied to various settings.
Most of the research on designing FPRAS’s for the GPF is focused on designing adaptive algorithms to produce sequences (cooling schedules) with minimum length while keeping the variances of estimators small (thus removing the need to have a-priori known bounds on variances). In contrast, the computation of the sequence of mean estimates, which dominates the total computation cost, is done by black-box MCMC estimators, with a priori known upper bounds on the mixing times of the chains. These upper bounds are often loose, and improving them for particular models is a challenging active area of research [4, 5, 9, 12, 28, 63]. In order to complement the adaptive cooling schedule and reduce dependence on a priori bounds on Markov chains’ mixing times, it seems necessary to design an adaptive procedure with theoretical guarantees for MCMC-mean estimation.
In this work we develop a doubly adaptive FPRAS, combining the adaptive cooling schedule with adaptive MCMC mean estimator that dynamically adapts the number of Markov chain steps to the observed underlying chain. Through rigorous theoretical analysis, we prove that our method outperforms the state of the art algorithms in several factors: (1) The computational complexity of our method is smaller; (2) Our method is less sensitive to loose bounds on mixing times, an inherent component in these algorithms; and (3) The improvement obtained by our method is particularly significant in the most challenging regime of high precision estimates. We demonstrate the advantage of our method in experiments run on classic factor graphs, such as voting and Ising models [5, 7, 15].

heading: 1.1 Preliminaries and Prior Work
text: Assume a sample space Ω, Hamiltonian function H : Ω→ {0} ∪ [1,∞), and inverse temperature parameter β ∈ R, referred to as inverse temperature. The Gibbs distribution on Ω, H(·), and β is
2The usage of the word “adaptive” here refers to algorithms which draw samples progressively and adapt their sample complexity based on empirical estimates until desired conditions are met, as it has been used in [34, 42] (see § 1.1), and should not confused with the work of [60].
then characterized by probability law
∀x ∈ Ω : πβ(x) . =
1 Z(β) exp ( −βH(x) ) . (1)
Here Z(β) is the normalizing constant or Gibbs partition function (GPF) of the distribution, with
Z(β) . = ∑ x∈Ω exp ( −βH(x) ) . (2)
Estimating the GPF Z(β), is computationally challenging, since typically the size of Ω is exponential in the number ofvariables, and the values of random terms in the sum have large variance (due to the exponential). The following problem has been extensively studied, and is the focus of this paper.
Problem 1. Given a domain Ω, a Hamiltonian function H , and a parameter β, design a Fully Polynomial Randomized Approximation Scheme (FPRAS) for estimating the partition function Z(β) . = ∑ x∈Ω exp ( −βH(x) ) . In other words, for user-supplied ε, the task is to produce an estimate Ẑ(β), such that with probability at least 1 − δ, we have (1 − ε)Z(β) ≤ Ẑ(β) ≤ (1 + ε)Z(β), in time polynomial in 1/ε, ln(1/δ), and all other problem parameters (e.g., the number of vertices in an Ising model, or neurons in an RBM).
All known scalable solutions to this problem rely on Monte-Carlo Markov-chain (MCMC) methods, and their execution cost is dominated by the total number of Markov chain steps they execute. We therefore follow past work, and analyze our algorithms in terms of number of the Markov chain steps.
TPA-Based Adaptive Cooling Schedules Building on extensive earlier work [6, 22, 67], the current state of the art is due to Huber and Schott [35], with Kolmogorov’s sharper analysis [42]. They introduce the paired product estimator (PPE), see definition 1.1, and apply the tootsie-pop algorithm (TPA) to adaptively compute a near-optimal cooling schedule, i.e., a sequence of inverse temperatures β0 < β1 < · · · < β`−1 < β` satisfying β` = β, and that Z(β0) is easy to compute, e.g., β0 = 0 is often convenient, since Z(0) = |Ω|. We thus define Q . = Z(β)/Z(β0) and estimate it using the paired product estimator.
Definition 1.1 (PPE [34]). Assume a cooling schedule β0, β1, . . . , β`. For each pair (βi, βi+1) in the schedule, we define two random variables, Xi ∼ πβi and Yi ∼ πβi+1 , all independent, and we then define fβi,βi+1 . = exp ( −βi+1−βi2 H(Xi) ) and gβi,βi+1 . = exp (βi+1−βi 2 H(Yi) ) . It is easy to verify that E[fβi,βi+1 ] = Z( βi+βi+1 2 )/Z(βi), and E[gβi,βi+1 ] = Z( βi+βi+1
2 )/Z(βi+1). We then define F .= ∏k i=1 fβi,βi+1 , G . = ∏k i=1 gβi,βi+1 . Letting µ̂ and ν̂ denote empirical estimates of E[F ] and E[G], respectively, the paired product estimator (PPE) is Q̂ .= µ̂/ν̂ .
Denote by Vrel[X] . = E[X
2]/E[X]2 − 1 = V[X]/E[X]2 the relative variance of a random variable X . The TPA schedule [35, 36] is generated by an adaptive algorithm, which, by a proper setting of parameters, outputs a cooling schedule guaranteeing constant Vrel[F ] and Vrel[G] (see alg. 3 in the supplementary material). Kolmogorov [42] presents a tighter analysis of Huber’s TPA method, and proves that with slight modifications (see alg. 4. in the Appendix) the schedule has a shorter length, while preserving constant relative variance for the paired product estimators (see thm. 1.1). In this paper, we use Kolmogorov’s algorithm, and we denote it by TPA(k, d). For completeness, both of Huber’s and Kolmogorov’s versions of TPA are presented in the Appendix.
We will use the following result in our analysis:
Theorem 1.1 ([42]). Let Hmax . = maxx∈ΩH(x), using TPA(k, d), k = Θ(logHmax) and d = 16 to generate cooling schedule (β0, β1, . . . , β`). W.h.p., we have ` = Θ(log(Q) log(Hmax)) and Vrel[F ]+1 = ∏` i=1(Vrel[fβi,βi+1 ]+1) = Θ(1) and Vrel[G]+1 = ∏` i=1(Vrel[gβi,βi+1 ]+1) = Θ(1).
Kolmogorov [42] nearly matches known lower bounds when given oracle access to near-independent samples, but leaves open the possibility of better use of the dependent sequence of samples generated by MCMC chains. This fertile ground is ill-explored, since if an approximate sampling oracle draws samples by running a chain for T steps, there is a factor T potential improvement.
MCMC Mean-Estimator Huber and Schott [35] assume unit-cost for exact sampling from each πβi , and Kolmogorov [42] extends their analysis to include the complexity of generating approximate samples with standard MCMC processes, assuming a priori upper-bounds on their mixing times. The main contribution of our paper is a specialized, adaptive, multiplicative MCMC-mean estimator for the TPA-based PPE. Our method is significantly more efficient than using standard black-box MCMC sampling for this problem, thus we improve the best-known method for estimating the GPF.
LetM be an ergodic Markov chain with state space S and stationary distribution π. Let τmix(ε) denote the ε-mixing time ofM, and define τmix . = τmix(1/4). Letting λ denote the second largest absolute eigenvalue ofM’s transition matrix, the relaxation time ofM is τrx . = (1− λ)−1, and it is
related to the mixing time τmix, by ( τrx(M)− 1 ) ln(2) ≤ τmix(M) ≤ ⌈ τrx(M) ln ( 2/√πmin )⌉ [45]. Let T be an upper bound on max{τrx(M), τmix(M)}. Consider any i.i.d. sampling concentration bound like Chebyshev’s, Hoeffding’s, or Bernstein’s inequalities [50], with, say, sample complexity mε. Using MCMC as a black-box sampling tool, we obtain the same precision estimation guarantees, with a computational cost of mε · τmix(ε/mε), which is equal to mε log(mε · ε−1) · T in the absence of exact values for τmix. Other concentration bounds compute the average over the entire trace of a Markov chain, and their complexity is dependent on known upper-bounds on the relaxation time [14, 38, 46, 50, 55], or function specific mixing time [59]. Note that since log( 12ε )(τrx − 1) ≤ τmix(ε) ≤ log( 1 επmin
)τrx, using these bounds is often more efficient, saving at least log(mε) steps.
Recently, Cousins et al. [16] introduce a novel Markov chain statistical measure, the inter-trace variance. The inter-trace variance depends on both the function being estimated and the dependency structure between nearby samples in the chain, and unlike the mixing time, it can be efficiently estimated from data. By using progressive sampling, Cousins et al. show an additive MCMC mean estimator whose complexity is proved in terms of inter-trace variance and they show it it less sensitive to prior knowledge of the input parameters, such as relaxation time and trace variance. Unfortunately due to a few technical problems, their result can not directly be used with the TPA method. Thus, in order to obtain a doubly adaptive algorithm for problem 1, we tailor their techniques to our setting, which requires developing new algorithms and analysis tools.

heading: 1.2 Our Main Contributions
text: • We present a specialized mean estimator method that significantly improves the state of the art computational complexity of computing the partition function of Gibbs distribution.
• While all rigorous MCMC-based estimates depend on some a priori knowledge of the Markov chain properties (such as bounds on its mixing or relaxation time), the complexity of our method is less dependent on these a priori bounds, and decays gracefully as they become looser.
• The improvement of our method is particularly significant in the more challenging high precision regime, where the goal is to compute estimates with very small multiplicative error.
• Our method improves the computational cost of prior work by replacing standard black-box MCMC mean estimators with an adaptive MCMC estimator, specially tailored to this problem.
• The analysis of our method relies on a novel notion of sample variance in a sequence of observations obtained by Markov chains runs, which we term the relative trace variance.
• We demonstrate the practicality of our method through experiments on Ising and voting models.

heading: 2 Algorithms
text: In this section, we develop two doubly-adaptive fully polynomial randomized approximation schemes providing more efficient algorithmic solutions to problem 1. The proof of all of the lemmas and theorems are presented fully in the supplementary material.
Notation and Setting Parameters We use the following notation throughout: We use capital letters to denote upper-bounds. e.g., T denotes an upper-bound on max(τmix, τrx), and Λ denotes a upperbound on the second absolute eigenvalue λ. We use GH,β to denote any Markov chain with Gibbs stationary distribution πβ , eq. (1). Having the Hamiltonian H , we denote its maximum and minimum values as Hmax and Hmin, i.e., Hmax . = maxx∈Ω{H(x)} and Hmin . = minx∈Ω{H(x)}. Having a
schedule (β0, β1, . . . , β`), the paired product estimators fβi,βi+1 , gβi,βi+1 , F = ⊗` i=1 fβi,βi+1 and
G = ⊗`
i=1 gβi,βi+1 are as in definition 1.1. When writing (β0, β1, . . . , β`) = TPA(k, d), we mean the cooling schedule is obtained from running alg. 4 in the Appendix, and we always set k = logHmax and d = 64, as these parameters are shown to produce a near-optimal schedule w.h.p. [42].
We first introduce a novel MCMC-based multiplicative mean estimation procedure RELMEANEST (see alg. 1), and analyze its computational complexity in terms of a new quantity, which we coin the relative trace variance (see definition 2.1). RELMEANEST receives as input a Markov chainM, a function f , and precision parameters ε and δ, and it outputs a multiplicative estimate of the expected value of the function w.r.t. the stationary distribution of the Markov chain. For simplicity, we may refer to it as RELMEANEST(M, f), leaving out the precision parameters. Letting (β0, β1, . . . , β`) = TPA(k, d), we first present PARALLELTRACEGIBBS, in which we invoke both RELMEANEST(GH,βi , fβi,βi+1) and RELMEANEST(GH,βi , gβi,βi+1) for each i = 1, 2, . . . , `− 1. We then present an often-more-efficient algorithm, SUPERCHAINTRACEGIBBS, which invokes RELMEANEST once each on F and G on a “super” product chain (see definition 2.2). We prove correctness of both PARALLELTRACEGIBBS and SUPERCHAINTRACEGIBBS, and bound their complexity in terms of the relative trace variance of the estimators. Furthermore, we prove SUPERCHAINTRACEGIBBS improves the computational complexity of the state of the art [42] (thm. 2.4 and corollary 2.6). Both of these algorithms have low dependence on tightness of mixing time: They receive as input an upper-bound on mixing or relaxation time T , but we show for ε ≥ ε0 their computation complexity is dominated by the true relaxation time τrel (of each Gibbs chain or the product chain).

heading: 2.1 Relative trace variance and RELMEANEST
text: In this section we introduce a new variance notion, the relative trace variance, which captures the computational complexity of MCMC-mean estimation with multiplicative precision guarantees. The relative trace variance depends on both the chainM and the function f , and it generalizes the relative variance, defined as Vrel[f ] . = V[f ]/E[f ]2, which depends only on f , and is used in i.i.d. regimes.
Definition 2.1 (Relative Trace Variance). For arbitrary τ , consider a trace of length τ of a Markov chainM, and a real-valued function f . OnM, we define the relative trace variance of f as
ReltrvτM[f ] . =
E[f̄( ~X1:τ ) 2 ]
(E[f̄( ~X1:τ )])2 − 1 ,
where ~X1:τ . = X1, X2, . . . , Xτ is a trace of length τ ofM, and f̄( ~X1:τ ) . = ( 1τ ) ∑τ i=1 f(Xi). We may drop the subscript when the chain is clear from the context.
The above definition is similar to what Cousins et al. coined as the inter-trace variance, denoted by trv(τ)(M, f), which they showed it captures MCMC-mean estimation with additive precision guarantees [16]. In fact, the two terms are related as
ReltrvτM[f ] = trv(τ)(M, f) (E[f̄( ~X1:τ )])2 .
Note that the two terms are not easily convertible without knowing the mean, E[f̄( ~X1:τ )]. Lemma 2.1. For any τ we have ReltrvτM[f ] ≤ Vrel[f ] . (3) Furthermore, for τ ≥ τrx(M) we have,
ReltrvτM[f ] = O
( τrx(M)
τ Reltrv
τrx(M) M [f ]
) . (4)
Lemma 2.1 enables us to compare the computational complexity of our algorithms with the state of the art [42]. In particular, using (3), we show our results improve the state of the art (which is in terms of Vrel), and using (4), we show that for high-precision estimations, the sample complexity of our algorithms only depends on τrx, which improves the state of the art (which is in term of T ).
The relative trace variance is a better analysis tool for estimating the GPF, because, unlike the intertrace variance, it leads directly to relative error bounds, rather than absolute error bounds.We now present some definitions which can also be found in standard MCMC textbooks, e.g., [45].
Definition 2.2 (Product Chain and Tensor Product Function). Consider k Markov chains {Mi}ki=1 each defined on state space Si and assume real valued functions {fi : Si → R}ki=1. The product chain M⊗1:k is defined on the Cartesian product of Si as follows: at any stepM ⊗ 1:k chooses i with probability
ωi (thus ∑k i=1 ωi = 1), and moves from (x1, x2, . . . , xi, . . . , xk) to (x1, x2, . . . , yi, . . . , xk), with the transition probability of moving from xi to yi inMi. The tensor product of {fi}ki=1, denoted by⊗ 1:k fi, is defined as (⊗ 1:k fi ) (x1, x2, . . . , xk) = ∏k i=1 fi(xi).
RELMEANEST Let T denote an upper bound on the relaxation time of a Markov chain M. RELMEANEST receives T , M, f and precision parameters ε and δ as input. Before it starts collecting samples, it runs the chain for a warm start (§ 2.1 of alg. 1). Starting from a minimum sample size m↓, it runsM for T ·m↓ steps, and collect samples X1, X2, . . . , XT ·m↓ . It then computes for j = 1, 2, . . . ,m↓, f̄j . = ∑j·T i=(j−1)·T+1 f(Xi); using them, it calculates an empirical estimate of the mean, µ̂, and an empirical estimation for the trace variance ofM and f , v̂. Based on these estimates, we derive an upper-bound on the current trace variance ui and relative error ε̂×i , and check whether is smaller than the user-specified error ε (lines 18-19). If so, we return the current mean estimate, otherwise we double the sample size and repeat.
Algorithm 1 RELMEANEST 1: procedure RELMEANEST 2: Input: Markov chain M, upper-bound on relaxation time T , real-valued function f with range [a, b], letting R = b − a,
multiplicative precision ε, error probability δ. 3: Output: Multiplicative approximation µ̂ of µ = Eπ[f ].
4: T ← ⌈
1+Λ 1−Λ ln
√ 2 ⌉ ; Λ′ ← ΛT . Choose T to be an upperbound on relaxation time
5: I ← 1 ∨ ⌊
log2 ( bR 2a2 · (1−ε)2 (1+ε)ε )⌋ ; α← (1+Λ ′)R ln 3Iδ (1+ε) (1−Λ′)bε ; m0 ← 0 . Initialize sampling schedule
6: Tunif ← ⌈ T · ln(1/πmin) ⌉ ; ( ~X0,1, ~X0,2)←MTunif (⊥) . Warm-start two chains for Tunif steps from arbitrary ⊥ ∈ Ω 7: for i ∈ 1, 2, . . . , I do 8: mi ← ⌈ αri ⌉
. Total sample count at iteration i; r is the geometric ratio (constant, usually 2) size 9: for j ∈ (mi−1 + 1), . . . ,mi do 10: ( ~Xj,1, ~Xj,2)← (T steps ofM starting at ~Xj−1,1, ~Xj−1,2) . Run two independent copies ofM for T steps
11: f̄( ~Xj,1)← 1
T T∑ t=1 f ( ~Xj,1(t) ) ; f̄( ~Xj,2)← 1 T T∑ t=1 f ( ~Xj,2(t) ) . Average f over T -traces
12: end for 13: µ̂i ← 1
2mi mi∑ i=1 ( f( ~Xj,1) + f( ~Xj,2) ) ; v̂i ← 1 2mi mi∑ i=1 ( (f( ~Xj,1)− f( ~Xj,2) )2 . Compute empirical mean; trace variance
14: ui ← v̂i + (11 +
√ 21)(1 + Λ ′ / √
21)R2 ln 3Iδ (1− Λ′)mi +
√ (1 + Λ′)R2v̂i ln 3I δ
(1− Λ′)mi . Variance upper bound
15: ε̂+i ← 10R ln 3Iδ
(1− Λ′)mi +
√ (1 + Λ′)ui ln 3I δ
(1− Λ′)mi . Apply Bernstein bound
16: µ̂×i ← (µ̂i − ε̂+i ) ∨ a+ (µ̂i + ε̂ + i ) ∧ b
2 . Optimal mean estimate
17: ε̂×i ← ((µ̂i + ε̂
+ i ) ∧ b− (µ̂i − ε̂ + i ) ∨ a
2µ̂×i . Empirical relative error bound
18: if (i = I) ∨ (ε̂×i ≤ ε) then . Terminate if accuracy guarantee is met 19: return µ̂×i 20: end if 21: end for 22: end procedure
The following theorem, shows the correctness of RELMEANEST and bounds its complexity.
Theorem 2.2 (Efficiency and Correctness of RELMEANEST). With probability at least 1 − δ, RELMEANEST will output µ̂ satisfying (1 − ε)µ̂ ≤ µ ≤ (1 + ε)µ̂. Furthermore, with probability at least 1− δ3I , the total Markov chain steps of RELMEANEST, m̂, obeys
m̂ ∈ O ln( ln baε δ )( T ·R µε + τrxReltrv τrx ε2 ) . (5)

heading: 2.2 Doubly adaptive algorithms: SUPERCHAINTRACEGIBBS and PARALLELTRACEGIBBS
text: Let (β0, β1, . . . , β`) = TPA(k, d), and consider a family of Gibbs chains GH,βi , each corresponding to some βi, and the paired product estimators F = ⊗` i=1 fβi,βi+1 G = ⊗` i=1 gβi,βi+1 . The TPA method is designed to ensure Vrel of the estimators are bounded, which can be employed by concentration bounds (e.g., Chebyshev’s bound) to guarantee the multiplicative error is bounded with high probability for a given sample size.
In order to generalize the same machinery for samples generated from a Markov chain using RELMEANEST, we need to bound the two terms appearing in eq. (5), which dominate the computational complexity of RELMEANEST. We refer to the first term, T ·R/µ, as the range term, and to the term τrxReltrvτrx as the trace variance term. Note that as ε becomes smaller, the trace variance term dominates the sample complexity of RELMEANEST, thus dependence on loose bounds T and R is dominated by dependence on true and a priori unknown values τrx and Reltrvτrx .
In order to ensure that the ranges of estimators are small, we prove that the length of each inversetemperature interval in the TPA schedule is w.h.p. small. Having a schedule (β0, β1, . . . , β`) we define and use the following notation: for 0 ≤ i ≤ `− 1, interval length ∆i . = βi+1− βi, maximum interval length ∆max . = maxi ∆i, and total length ∆ . = β` − β0 .
Lemma 2.3. Let z(β) .= ln ( Z(β) ) , and let βi, βi+1 be two consecutive points generated by TPA(k, d). For arbitrary ε ≥ 0, we have: 1. P(z(βi)− z(βi+1) ≤ ε) ≥ (1− exp(−εk/d))d. 2. P ( ∆i ≥ ε/E[H(x)] ) ≤ d exp(−εk/d), where E[H(x)] is taken w.r.t. x ∼ πβi+1 .
SUPERCHAINTRACEGIBBS Let G⊗ the product of GH,βis with uniform weights i.e., ωi = 1` ,∀i (see definition 2.2). SUPERCHAINTRACEGIBBS calls RELMEANEST(G⊗, F ) and RELMEANEST(G⊗, G), with appropriate parameters, and simply outputs the ratio of the two estimates (see alg. 2, left).
Algorithm 2 SUPERCHAINTRACEGIBBS and PARALLELTRACEGIBBS 1: procedure SUPERCHAINTRACEGIBBS(. . . ) 2: (β0, β1, . . . , β`)← TPA(k, d)a 3: ε′ ← ε
2+ε ; δ′ ← δ 2
4: for i ∈ 1, 2, . . . , ` do 5: fi(x) . = exp(−βi+1−βi
2 H(x))
6: gi(x) . = exp( βi−βi−1 2
H(x)) 7: end for 8: F .= ⊗` i=1 fi; G . = ⊗` i=1 gi
9: G⊗ ← ⊗`
i=1 GH,βi , with ωi = 1 ` ,∀i
10: Rf ← exp(−β−β02 Hmin)− exp(− β−β0 2 Hmax) 11: Rg ← exp(β−β02 Hmax)− exp( β−β0 2 Hmin) 12: µ̂← RELMEANEST(G⊗, Rf , T, F, ε′, δ′) 13: ν̂ ← RELMEANEST(G⊗, Rg, T,G, ε′, δ′) 14: return Ẑ ← ν̂
µ̂
15: end procedure
16: procedure PARALLELTRACEGIBBS(. . . ) 17: (β0, β1, . . . , β`) = TPA(k, d) 18: ε′ ←
√̀ 1+ε−1√̀ 1+ε+1 ; δ′ ← δ 2`
19: for i ∈ 1, 2, . . . ` do 20: fi(x) . = exp(−βi+1−βi
2 H(x))
21: gi−1(x) . = exp( βi−βi−1 2 H(x)) 22: Rf ← exp(−βi+1−βi2 Hmin)− exp(− βi+1−βi 2 Hmax) 23: Rg ← exp(βi+1−βi2 Hmax)− exp( βi+1−βi 2 Hmin) 24: µ̂i ← RELMEANEST(Gi, Rf , Ti, fi, ε′, δ′) 25: ν̂i ← RELMEANEST(Gi, Rg, Ti, gi, ε′, δ′) 26: end for 27: return Ẑ ← ∏` i=1 ν̂i µ̂i 28: end procedure
ak = Θ(logHmax) and d = 64 as in [42]
We now show the correctness and efficiency of SUPERCHAINTRACEGIBBS. Let τprx denote G⊗’s true (and unknown) relaxation time and T a known upper-bound on it (T ≥ τprx), ε and δ are user specified precision parameters. For simlicity of presentation we use the following notation to refer to relative ranges: relR = Range(F )/µ+ Range(G)/ν, where µ = E[F ] and ν = E[G] . Theorem 2.4. With probability at least 1− δ, it holds that the total number m̂ of Markov chain steps taken by SUPERCHAINTRACEGIBBS is upper-bounded by
Õ ( ln ( 1
δ
)( T · relR
ε + τprx ·
( Reltrv
τprx G⊗ (F ) + Reltrv τprx G⊗ (G) ) ε2 )) .
Lemma 2.5. Defining α1 = √
Z(β0) Z(β0−∆max) , we have:
Range(F ) µ ≤ α1
√ Q
exp(∆Hmin) and
Range(G) ν ≤ α1
√ exp(∆Hmax)
Q .
Using lemma 2.5 and thm. 2.4, we identify ε0 such that for ε ≤ ε0 the trace variance term in will the dominate computational complexity of SUPERCHAINTRACEGIBBS. In order to make a fair comparison with the state of the art [42] we employ eq. (3) of lemma 2.1. Finally we use thm. 1.1 and conclude: Corollary 2.6. Let α1 be as in lemma 2.5, τmax . = maxi τi and ε0
. = (τprx/T ) ·(√
exp(∆Hmin) Q +
√ Q
exp(∆Hmax)
) · α1. When ε ≤ ε0, the number of Markov chain steps of
SUPERCHAINTRACEGIBBS is dominated by Õ(`τmax) .
PARALLELTRACEGIBBS For i = 1, 2, . . . , ` − 1, PARALLELTRACEGIBBS (alg. 2, right) runs RELMEANEST(GH,βi , fβi,βi+1) and RELMEANEST(GH,βi , gβi,βi+1) independently. We show the computational complexity of PARALLELTRACEGIBBS in thm. 2.7.
For i = 1, 2, . . . , `, assume τi is the true (unknown) relaxation time of GH,βi and Ti is a known bound on it. For simplicity of presentation we use the following notations: relRi . = Range(fβi,βi+1)/µi + Range(gβi−1,βi)/νi, where µi = E(fβi,βi+1) and νi = E(gβi,βi+1). Theorem 2.7 (Efficiency of PARALLELTRACEGIBBS). With probability at least 1 − δ, it holds that the total number m̂ of Markov chain steps taken by PARALLELTRACEGIBBS is upper-bounded by
Õ ( log ( `
δ )∑̀ i=1 ( ` · Ti · relRi ε + `2 ε2 τi · ( ReltrvτiGH,βi (fβi,βi+1) + Reltrv τi GH,βi (gβi−1,βi) ))) .
Furthermore, for all 1 ≤ i ≤ `, Range(fβi,βi+1)/µi ≤ `1/ log(n) and Range(gβi−1,βi)/νi ≤ `α0(i)/ logn, where α0(i) = (Hmax/2E[H(x)])− 1, for x ∼ πβi .
PARALLELTRACEGIBBS and SUPERCHAINTRACEGIBBS make different computational complexity tradeoffs. PARALLELTRACEGIBBS is usually slower than SUPERCHAINTRACEGIBBS, because in each iteration i = 1, 2, . . . , `, the mean estimator must acquire a higher-precision estimate so that all estimators together achieve an ε-δ relative-error guarantee. Relaxation times (true values and their upper-bounds) appear in a sum in the complexity of PARALLELTRACEGIBBS, whereas they appear in a maximum in SUPERCHAINTRACEGIBBS ( ∑` i=1 τi vs. maxi=1,...,` τi). Furthermore, dominance of the trace variance terms in both of these algorithms occur at different values of ε. A comparison of the complexity of these algorithms, in the high-precision regime, with Kolmogorov’s TPA + PPE (which uses MCMC as a black box) is presented in table 1.

heading: 3 Experimental Results
text: In this section we report our experiment results, comparing the performance of the two versions of our doubly adaptive method (alg. 2), to the performance of the state of the art algorithm in [42].
Setup. We run the experiments using the single site Gibbs sampler (known also as the Glauber dynamics) on two different factor graph models:
(A) The Ising model on 2D lattices. Having a 2-dimension lattice of size n× n, the Hamiltonian is defined on n2 random variables having values ±1 and their dependency is represented by the Hamiltonian: H(x) = − ∑ (i,j)∈E 1(x(i) = x(j)). We run the algorithms on lattices of sizes 2× 2,
3 × 3, 4 × 4, and 6 × 6. For each lattice, the parameter β ≥ 0 is chosen below the critical inverse temperature at which it undergoes a phase transition. We use known mixing time bounds for high temperature Ising models [3] (see fig. 1 and A.6. of supplementary material).
(B) The logical voting model. For a parameter n , we have 2n + 1 random variables: the query variable Q ∈ {−1, 1}, and the voter variables T1, T2, . . . , Tn and F1, F2, . . . , Fn all in {0, 1}. The factors have 2n+ 1 weights, ω, ωTi , ωFi , i = 1, . . . , n. The Hamiltonian is:
H(Q,T, F ) = ωQmax i Ti − ωQmax i Fi + n∑ i=1 ωTiTi + n∑ i=1 ωFiFi , where ω, ωTi , ωFi ∈ [−1, 1]
The parameters are reported in fig. 2. We follow De Sa et al. [17] and use hierarchy width to derive upper bounds on mixing times. To make a fair comparison, we always run the TPA algorithms once, and with the parameters given in [42]. At each iteration of RELMEANEST, the sample size is extended with geometric ratio 1.1 (see alg. 1 line 8). All code is available at https://github.com/ zysophia/Doubly_Adaptive_MCMC.
Results: Our experiments demonstrate the practical advantages of our doubly adaptive method, validating our theoretical analysis.
(1) We first compare the complexity of our algorithms to Kolmogorov’s algorithm. Our experiments show the superiority of both versions of our methods on different models and various sets of parameters. Figure 1 demonstrates the superiority of our methods on the Ising model for various sets of parameters, and in figs. 2a and 2c for the voting model, when ε is fixed and Z(β) is varying (fig. 2c), and when Z is fixed and ε is varying (fig. 2a). All of these hold while the precision of our algorithms beats [42] as ε→ 0 (fig. 1c). (2) To demonstrate the advantage of using the relative trace variance, in contrast to the relative variance, we run both of our algorithms using a simpler mean estimator which only uses progressive sampling, and we compare the results. This is done by setting T ← 1 in line 4 of RELMEANEST. In Figure 2b, we show the effectiveness of trace averaging, since both SUPERCHAINTRACEGIBBS and PARALLELTRACEGIBBS beat their simplified versions (T ← 1) after 1/ε passes a certain threshold. This is consistent for different parameters of the voting model.
(3) Comparing the performance of SUPERCHAINTRACEGIBBS and PARALLELTRACEGIBBS, we observe that in all of our experiments SUPERCHAINTRACEGIBBS has better performance than PARALLELTRACEGIBBS. In fig. 2b, we show the trace variance term PARALLELTRACEGIBBS becomes dominant earlier as 1/ε grows, thus it performs better in this perspective. This is consistent with our theoretical findings, because the ranges of estimators in PARALLELTRACEGIBBS are smaller than the ranges used in SUPERCHAINTRACEGIBBS.

heading: 4 Conclusions: advantages and limitations of proposed algorithms
text: We develop a doubly-adaptive MCMC-based estimator for the partition function of Gibbs distributions, which resolves a major impediment of prior methods that use MCMC as a black-box sampler. We show, both theoretically and experimentally, that our method requires substantially fewer MCMC steps than the state-of-the-art method. The better performance is due to several factors, which all stem from the use of an adaptive MCMC mean estimator instead of a standard "black-box" MCMC estimate. The complexity of the adaptive MCMC process depends on the (smaller) trace, rather than stationary, relative variances, and on relaxation times instead of mixing times. It is also less sensitive to weak upper-bounds on mixing and relaxation times.
In particular, Kolmogorov’s method requires Θ( /̀ε2) approximately independent samples, where ` is the length of cooling schedule. This requires tight convergence (total variance distance of O(ε2/`) from stationary) for each sample, which adds a multiplicative ln `ε2 , with ` = Θ(lnQ lnHmax), to its complexity (see column 3 of table 1 and [42], theorem 9). In contrast, our doubly adaptive method only depends on relaxation times, which do not depend on ε.
Limitations. While significantly improving the state of the art, our methods suffer from a several limitations. In SUPERCHAINTRACEGIBBS, the major limitation is the dependence on the relative ranges of F and G, which can be large, especially when the Hamiltonian range is large. Another issue is that the product chain’s mixing time is dominated by `max{τi}`i=1, as opposed to ∑` i=1 τi. While PARALLELTRACEGIBBS circumvents these issues by estimating each factor of the telescoping product independently, it fails to beat SUPERCHAINTRACEGIBBS’s efficiency in general, due both to the union bound and the higher-precision guarantees required for each subproblem. Improving performance further will likely require new estimators with smaller ranges and relative trace variances.
Statement of Broader Impact. While probabilistic graphical models as other machine learning methods that rely on MCMC estimations continue to grow in importance and popularity. But running the MCMC to theoretical convergence guarantees is often prohibitively expensive, while running it to apparent convergence is methodologically unsound, particularly in the modern context, where public confidence in machine learning systems is continuously eroded by ethical, accuracy, and safety failures. Our work attempts to bridge the gap between the definite, elegant and theoretically sound analytic methods, and efficiency-focused practical utility, as we seek to reduce proof-burden, while maintaining theoretical guarantees of accuracy, with adaptive methods that bound efficiency in terms of (potentially unknown) convergence rate metrics and variances.
Acknowledgements. Shahrzad Haddadan is supported by NSF Award CCF-1740741. Cyrus Cousins and Eli Upfal are supported by NSF grant RI-1813444 and DARPA/AFRL grant FA8750. The authors are thankful to anonymous reviewers of NeurIPS 2021 for several valuable inputs.

Review summaries:
The authors propose a doubly "adaptive" MCMC algorithm to approximate the partition function in gibbs distributions. In related work, the partition function of Gibbs distributions are computed using paired product estimators [1]. Such estimators use a cooling schedule $[\beta_{(\cdot)}]$ to estimate two product functions $F$ and $G$ using stationary samples from the Gibbs distribution at each temperature. The ratio of these two products is then used to compute the partition function.  In the current work, the first "adaptive" refers to the adaptive TPA cooling schedule to decide $[\beta_{(\cdot)}]$. The second refers to the automatic number of steps used in the MCMC MeanEstimator which is used to estimate the two products mentioned above $F$ and $G$.  [1] - Huber, Mark. "Approximation algorithms for the normalizing constant of Gibbs distributions." The Annals of Applied Probability 25.2 (2015): 974-985.

The objective of this work is to define a procedure for calculating Gibbs distribution normalization constants. The calculation of normalisation constants is a central problem in many fields and most known methods depend on many hyper-parameters that are difficult to optimise. The most efficient methods (such as "paired product estimator" extending "multiple  sampling" of Valleau and Card) are based on annealing schemes. The normalisation constant is obtained in this case as a product, each term of the product being calculated by means of an MCMC algorithms. Several strategies have been proposed to optimise the annealing schedule. The optimisation of the annealing schedule requires a measure of the complexity of generating samples with the MCMC algorithm, which typically involves using upper bounds on their mixing times. Despite progress in this direction, the estimation of chain mixing times is a difficult subject: the known bounds are, with few exceptions (of very simple models with a lot of symetries), very pessimistic. The paper presents an imporoved mean estimator for comuting the nomalizing constant of a Gibbs distribution. The proposed algorithm is claimed to improve SOTA results in this field. Compared to other techniques, the proposed algorithm is "less dependent" (the authors words) on the tightness of the bounds use to measure the mixing time. The main improvements bought by the authors is clained to be high-accuracy computaitons. 

In this paper, the authors introduce a new method to estimate the partition functions of Gibbs distributions. They try to improve the classical mean estimate computation by using an adaptive MCMC mean estimator. And they justify the priority of their methods through theoretical analysis and experiments.