Fast Doubly-Adaptive MCMC to Estimate the Gibbs Partition Function with Weak Mixing Time Bounds
Authors: Shahrzad Haddadan
Yue Zhuang
Cyrus Cousins
Eli Upfal
Venue: NeurIPS 2021 Poster (2021)
Keywords: MCMC
Gibbs distribution
partition functions

Abstract: We present a novel method for reducing the computational complexity of rigorously estimating the partition functions of Gibbs (or Boltzmann) distributions, which arise ubiquitously in probabilistic graphical models. A major obstacle to applying the Gibbs distribution in practice is the need to estimate their partition function (normalizing constant).  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 estimate computations use MCMC as a black-box to draw approximately-independent samples. Here 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 components 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 models and Ising models. 


Fast Doubly-Adaptive MCMC to Estimate the Gibbs Partition Function with Weak Mixing Time Bounds
Shahrzad Haddadan
Brown University
The Data Science Initiative shahrzad haddadan@@maal......... com
Yue Zhuang
Cyrus Cousins r Brown University
Brown University
The Data Science Initiative yue_ zhuangl @brown. edu
Department of Computer Science
cyrus_couss coriusscousisstene www.bowoooo.oonnn edu
Eli Upfal
Brown University
Department of Computer Science
weazzrruufal@brrrre wawaalbboown.ee. edu
Abstract
We present 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 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 a black-box draw approximate samples. We develop 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: (I) 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 classic factor graphs, such as voting models and Ising models.
Introduction
The Gibbs (Boltzmann) distribution is a family of probability distributions of exponential form. First introduced in the context of statistical mechanics (25). Gibbs distributions now ubiquitous in a variety of other disciplines, such as chemistry 24 31, economics (4) and machine learning. Gibbs distributions are typically used to model the global state of a system function of collection of interdependent random variables, each representing local states in the system. The dependencies
indicates randomized ordering and equal contribution.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

in the system are modeled by Hamiltonian function, and the probability distribution is inversely proportional exponent of the Hamiltonian scaled by the temperature (see eq. $ T.I.
Gibbs distributions provide potent statistical inference tools in many machine learning applications.
They appear in prohabilistic graphical models 41 49 [11, including restricted Boltzmann machines
43 [6], Markov random fields [40 47), and Bayes networks [32]. and are applied in the analysis of
images and graphical data [2] 23
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 compute, estimate, its partition function (normalizing constant), henceforth written GPF. The partition function is defined over the Cartesian product of supports of (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
42
52 63 67
Like other MCMC methods, here various heuristics are used. The most well-known heuristics
the annealed importance sampling [39 52 64] or convergence diagnostics methods [I0 II 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 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 either unreliable computationally
intractable (see 12. section 4] or g 331
On the other hand, theoreticians study this problem by designing Fully Polynomial Randomized Approximation Schemes (FPRAS) (see problem [[]. The state of the art FPRAS for estimating the GPF is a multi-stage algorithm involving 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-bounde on various measures such as variances of estimators or mixing times of Markov chains. In static algorithms, these upper-bounds adaptive 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.
given priori, and
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 priori known upper bounds 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 [ 5
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 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 method is particularly significant in the most challenging regime of high precision estimates. We demonstrate the advantage of our method in experiments run classic factor graphs, such as voting and Ising models (5 15
1.1 Preliminaries and Prior Work
Assume sample space Q, Hamiltonian function H Ωt + (1, and inverse temperature parameter 3 R. referred to as inverse temperature. The Gibbs distribution on Ω, (-). and / is
"The 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 R4 2p (see*) [.I. and should not confused with the work of (60).

---

then characterized by probability law
= 1
Vo E E mp(x)=
exp( exp(-3H(xx)
(1)
((3)
Here 7(3) is the normalizing constant Gibbs partition function (GPF) of the distribution, with
( - Σ exp((33((x))
(2)
JE
Estimating the GPF 2(3), is computationally challenging, since typically the size of Qt is exponential in the number ofvariables, and the values of random terms the sum have large variance (due the exponential). The following problem has been extensively studied, and is the focus of this paper.
Problem 1. Given domain Q. Hamiltonian function H, and parameter 3, design Fully Polynomial Randomized Approximation Scheme (FPRAS) for estimating the partition function Z(3) = Ezea exp(-3H (x)) In other words, for user-supplied e, the task is to produce
estimate
(B). such that with probability at least 1 - 6, we have (1 -1- e ((8) < 7 S (1 + + in time polynomial in 1/c, In(V/), and all other problem parameters (e.g., the number of vertices in 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 [.]] and apply the tootsie-pop algorithm (TPA) to adaptively compute a near-optimal cooling schedule, i.c., sequence of inverse temperatures 30 < B1 <... < Be-1 33 satisfying Br and that ((8) is easy compute, e.g., 90 0 is often convenient, since Z(0) = [2). We thus define Q - (BB) and estimate it using the paired product estimator.
Definition 1.1 (PPE 344). Assume a cooling schedule Bor B2.....B. Be. For each pair (Be, B,-1) in the schedule, define two random variables, X, TG, and Y, TT all independent, and we then define f3,,944 cxp( 8 8. ((x))) and ga..8441 = cxp 1B-1-- H(:))). It is easy to verify that [[[ Z 4 1 1 / and Elg3,.041l Z( 2 1 We then define F = I[ f8,,8 G = IT 1= 98..8...1 Letting ja and i denote empirical estimates of E(F] and G respectively, the paired product estimator (PPE) is Q n/3.
Denote by [xa[XX]
v|xl/E|xT² the relative variance of a random variable X. 35 36 is generated by an adaptive algorithm. which, by proper setting of
||x"//EXX|=
The TPA schedule
parameters, outputs cooling schedule guaranteeing constant Vra[F)) and [G] (see alg. in the supplementary material). Kolmogorov 42 presents tighter analysis of Huber's TPA method, and proves that with slight modifications (see alg. 4. the Appendix) the schedule has shorter length, while preserving constant relative variance for the paired product estimators (see thm. In this paper, we use Kolmogorov's 's algorithm, and we denote it by TPA(k; k; Huber's and Kolmogorov's 's versions of TPA are presented in the Appendix.
d). For completeness, both of
We will use the following result in our analysis:
Theorem 1.1 (22)). Let Hmax max Hx),, using TPA(K, d). = (log H mmax) and d = 16 to generate cooling schedule (So, 31 , ,3).. W.h.p., have E O Q) log H max) )) and
Viel [F] (Vel
(11) and V.ddVWalGGH II (Vralgg, 1+1) O())
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 an approximate sampling oracle draws samples by running chain for T steps, there is a factor T potential improvement.
3

---

MCMC Mean- -att-mstimmtor Huber and Schott [35] assume unit-cost for exact sampling from each Tg, and Kolmogorov 42 extends their analysis include the complexity of generating approximate samples with standard MCMC processes, assuming a priori upper- houeds on their mixing times. The main contribution of our paper is 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.
Let M be an ergodic Markov chain with state space S and stationary distribution я. Let Tmix (e) denote the e-mixing time of M, and define Tmix = +mmix(/al Letting À denote the second largest absolute eigenvalue of M'ss transition matrix, the relaxation time of M is Tix = (1 A)= and it is
related to the mixing time Tmix. by - [n(2) S Tmix S [Tx In (2/
45
Let T be an upper bound on max[Trx(MM (M)).
Consider any i.idd sampling concentration bound like Chebyshev's, Hoeffding's s, or Bernstein's 's inequalities [50], with, say, sample complexity my. Using MCMC as a black-box sampling tool, we obtain the same precision estimation guarantees, with computational cost of m. (mix(E/me),
which is equal to me log(m, - T in the absence of exact values for Timix.
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( + (Trx - 1) S Tmix S log( ETTratn Tyx
using these bounds is often more efficient, saving at least log(me) steps.
Recently, Cousins al. 16 introduce 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, 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 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 [ we tailor their techniques to our setting, which requires developing new algorithms and analysis tools.
1.2 Our Main Contributions
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 priori bounds, and decays gracefully as they become looser.
The improvement of our method is partieularly 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 observa- tions 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.
2 Algorithms
In this section, we develop two doubly-adaptive fully polynomial randomized approximation schemes providing more efficient algorithmic solutions to problem T 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- hound max(Tmix: Txx), and A denotes a upper- bound on the second absolute eigenvalue X. We use GH, to denote any Markov chain with Gibbs stationary distribution Tg, eq. (I). Having the Hamiltonian H, we denote its maximum and minimum values as H Hmxx and H min- i.e., H Haxx = max ( (x)) and Hmin min {H)))
(x)). Having a
schedule (Bo, St.... Be), the paired product estimators fa..a I3: F = 8 fo. 41 and

---

G 8 -=1 9344 are as in definition 1.1 cooling schedule is obtained from running alg. 4 in the Appendix, and we always h log Hmax and d = 64, as these parameters are shown to produce near-optimal schedule w.h.p. 42
When writing (Bu. .......,e TPA(k, d). we mean the
We first introduce novel MCMC-based multiplicatiie mean estimation procedure RELMEANEST (see alg. [ and analyze its computational complexity in terms of new quantity, which we coin the relative trace variance (see definition 2.1 RELMEANEST receives as input Markov chain M, function f. and precision parameters & and 6, 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 M. f), leaving out the precision parameters.
Letting (30, ................................ , 3z) TPA(K, d), we first present PARALLELTRACEGIIE in which we invoke both RELMEANEST((GH.R A. S f and REL BELMERANET(GIIL. (a,, for each i 1,,....,, which invokes
We then present an mofeeeooe-efficint algorithm, SUPERCHAINTRACCGIIRS RELMEANEST once each on F and G on "super" product chain (see definition (.2). We prove correctness of both PARALLELTRACIGIBBS and SUPERCHAINTRACEGIBBS and bound their complexity in terms of the relative trace variance of the estimators, Furthermore, we prove SUPERCHAINTRACEIIISSS
improves the computational complexity of the state of the art 22 (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 & 2 Eu their computation complexity is dominated by the true relaxation time Trel (of each Gibbs chain or the product chain).
2.1 Relative trace variance and RELMEANEST
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 frace variance depends on both the chain M and the function J, and it generalizes the relative variance, defined as Vnalf] = VU[fE(]]", which depends only f, and is used in i.i.d. regimes. Definition 2.1 (Relative Trace Variance). For arbitrary т, consider a trace of length t of a Markov chain M. and a real-valued function f. On M, we define the relative trace variance of f as
Reltev = [EIf(X;:))1>*
where X1.r = =x,,xx..... X, is a trace of length I of M, and may drop the subscript when the chain is clear from the context.
(t) We
The above definition similar to what Cousins et al. coined as the inter-trace variance, denoted by try(r) (M,, f). which they showed it captures MCMC-mean estimation with additive precision guarantees [16]. In fact, the two terms are related as
Reltrv (rv( (M, (BIf(X1++)))
Note that the two terms are not easily convertible without knowing the mean, B[J( (8:=))
Lemma 2.1. For any T we have
Reltrv Mill Wielf]]
(3)
Furthermore, for T 2 Tex(M) we have,
Reltrvi/lf]= =
Trx(M
Reltrv (M)
To M [f]
(4)
Lemmal 2.T enables 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 V..),, and using (4) show that for high-precision estimations, the sample complexity of our algorithms only depends on Trx- which improves the state of the art (which in term of T))
The relative trace variance is a better analysis tool for estimating the GPF, because, unlike the inter- trace variance, it leads directly to relative error bounds, rather than abvolute error bounds. We now present some definitions which also be found in standard MCMC textbooks, e.g., [45].
5

---

Definition 2.2 (Product Chain and Tensor Product Function). Consider k. Markov chains {M;
each defined on state space St and assume real valued functions {fi S 3 +R The product chain Mix is defined on the Cartesian product of S; as follows: at any step M chooses with probability wi (thus E wi =1), and moves from (21,32,..- ......... "...,,i;....... :.ii:..,220 r*)
to (21, 22,..., -V..-..,, with the transition probability of moving from 2; yi in M,. The tensor product of {fi) denoted by
Dick fú defined as (8 1sk: (21, 22...., .x) = II. f((z)).
Let T denote upper bound the relaxation time of a Markov chain M. MEANESS receives T, M,, f and precision parameters & and as input. Before it starts col- lecting samples, it runs the chain for
RELMEANEST
start ($ 22 of alg. [T. Starting from minimum sample size m it runs M for T ' steps, and collect samples X1,
[lg.[])
It then computes for j = 1,2, m J, = 2i=(j-1)-T) fX;;;; using them, it calculates an empirical estimate of the mean, ft, and empirical estimation for the trace variance of M and f. 0. Based on these estimates, we derive an upper-bound on the current trace variance u, and relative error E and check whether smaller than the user-specified error (lines 18-19). If so, return the current mean estimate, otherwise we double the sample size and repeat.
Algorithm 1 RELMEANEST
1: procedure RELMEANEST
Input: Markov chain M, upper-bound on relaxution time T,, real-valued function with range Ja, 6). letting a. multiplicative precision emor probability
Output: Multiplicative upproximation of B.
4: T In 2 - A
D Choose T to be an upperbound on relaxation time
I
1V
log
a
R ((1is), 11 N
my
Initialize sampling schedude
Tuni T In((/e/emm))
] (No... Xo21 M²
Warm-start two chains for Tunit steps from arbitrary 1. 9 Total sample count iteration the geometric ratio (constant, usually 2)
1...... do aar'
m.
for (m--1 ..
10:
(X X,,2) steps starting at E
Run independent copies of for T steps b Average over T-truces
F
11:
f(X,a)
Σ 1 (Ra((t);; f(X,2) + J
Σ / (8g2(())
12:
end
13:
A
2m Σ (f( f(x;a)): 82 2m.. Σ (C f(Xxa))"
& Compute empirical mean; trace variance
(11 (1 X 2f
A
14:
34 4- à
> Variance upper bound
Nmm;
Mjm;
Et 10 o
(1 A'),,
15:
Nim.
V
(1 N'mm
Apply Bernstein bound
e!) va+a
Ab
16:
A
> Optimal mean estimate
((s: + Ett A b- (A4- +t)Va
17:
E"
Empirical relative error bound
20"""
18:
if (E" <e) then
Terminate if accuracy guarantee
19:
return P."
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 8, RELMEANEST will output ğt satisfying (1 - E < μ < (1 + E)G. Furthermore, with probabil- ity at least 31, the total Markov chain steps of RELMEA,EST, in, obeys
In
T T
Trx Reltry Tix
m E
O
In
ar
+
(5)
8

---

2.2 Doubly adaptive algorithms: SUPERCHAINTRACEEIIESS and PARALLELTRACEGIEBS
Let (Bu.S...... 3e) (P, d), and consider family of Gibbs chains GH,S,, each corresponding to some 84 and the paired product estimators F' = 8 method is designed to ensure Vel of the estimators are bounded, which can be employed by concen- tration bounds (e.g., Chebyshev's bound) to guarantee the multiplicative error is bounded with high probability for given sample size.
i- fp..0... G X am B..4 The TPA
In order to generalize the same machinery for samples generated from a Markov chain using RELMEANEST, need to bound the two terms appearing in eq. ( which dominate the com- putational complexity of MEANESS. We refer to the first term, 7 R2/1. ju, as the range term, and to the term Txx Re the trace variance term. Note that becomes smaller, the trace variance term dominates the sample complexity of RELMEANEST, thus dependence on loose hounds and R dominated by dependence on true and a priori unknown values Tix and Reltrv"
In order to ensure that the ranges of estimators are small, we prove that the length each inverse- temperature interval in the TPA schedule is w.h.p. small. Having schedule (Bo, .,...., .B)) we define and use the following notation: for - i i L interval length A; Bi+1 - Br. maximum interval length A max = max, A,, and total length A = 31 Bo
Lemma 2.3. Let z(3) In (()),, and let Br Bit1 be two consecutive points generated by (PA d). For arbitrary & 2 0, we have:
1. ((z 2(B,+1) I <) 2 (1 exp(-ek/d)
2. P (A, 2 s/E[H(x))) I d dexp(- =k/d)) where (x)( is taken w.E.t. x ~ TB:11
Let 9 the product of GH, with uniform weights ie., wi = Vi (see definition 2.2). CUPERCHAINTRACFGIIN calls RELMEANEST F) and RELMMANNETT(S G), with appropriate parameters, and simply outputs the ratio of the estimates (see alg. 2 left).
CUPPERCHAINTAACEGIIO
Algorithm 2 SUPERCHAINTRACEGIRBS. and PARALLELTRACEGIBBS.
1: procedure SUPERCHANNRRACGGIBSS.
16: procedure GRAALELTRAAFFIBBSS...
3c) + TPA(k, di] 2: (8o, ..)
17:
(8o, Uo.........8 3t) = (PA(k,
3: g'
18:
e' (TIE 1,5' 2f
V1+7-1
4: for i E 1,2...., do
19:
for 1,
8 (a)) 3 6: (xxx) - exxp 8 (()))
5: fi(z) - expl
20:
fi(z) = expl 2:11 ((()))
21:
g--( = expl 8.-8. 11 (1))
7:
end for
R2 + exp(- Bet1 3
2 Hmmn) - expl Beos 98 Hmax)
8: F - 8' G= D'_1 19
22:
Rg + exp( 3:+1- Hmm)) expl - Hmn))
9: Go 8i_1 G11.3, with
Vi
23:
10: R2 + exp(- 8-20 Hmmn) - exp( #- Hmox)
24:
á. < RELMEANEST G. R2 To, 2,,
11: R2 expl app 2- Hmax) exp( 8- Hmm)
2
25:
D. + RELMEANNET(GG, R., Tr, E'.''
26: end
12: A + RELMEANEST(CG R,, F
27:
return 2 + --I' &
13: RELMEANESTIC R., R,,,,,,,,
28:
end procedure
14: return 2 7-
15: end procedure
"k = Olog Hmax) and = 64 as in 42
We now show the correctness and efficiency of SUPERCHHAINTRACGIIBBS. Let Tprx denote S true (and unknown) relaxation time and T' a known upper- bound on (T' 2 Tpxx), and are user specified precision parameters. For simlicity of presentation we use the following notation to refer to relative ranges: relR = Range(F)/H
(G)ggeGG)/vv where μ = E[F] and v E[G..
Theorem 2.4. With probability at least 1 5, it holds that the total number in of Markov chain steps is upper-bounded by
taken by SUPERCHAINTRACEEGIII
ò 6(11 In
(
T relR
+ Tpex (Reltrv
(F) ReRlRe Go
Turx
go (G))
E2
Lemma 2.5. Defining 01
Z Z180o Amar] we have:
Z(Bo)
Rangel F)) < a1 exp( and
Rango( ()) < on
expl Amsx)

---

PARALLELTRACEUIBSS
SUPERCHAINTRACEUUBBB
Tha + PPE 42
In glu Hmax
Roltrven.o (fi)
Fpex PRoltryTpoxpe + Reltry (0])
Tj (Yea(F) P) + Vaa(G))
Reltrven. R1r,3, (9c)
(emux(T)
In Hanx
T:
Table I: Comparison of the number of Markov chain steps, when E is adequately small. In all columns, a multiplicative factor of 1/2² is omitted ease presentation, and = In Q. Note that computational complexity of both PARALLELTRACEEIIS and SUpERCHAiNTRaCECIAA only depends true relaxation times, denoted by 7s, and TPA + PPE method's complexity dependent on their upper bounds, denoted by T.
Using lemma 2.5 and thm. 2.4 identify Eo such that for S the trace variance term in will the order to make a fair comparison with the state of the art [42] employ eq. (3) of lemma 2.1 Finally use thm. [.] and conclude:
dominate computational complexity of SUPERCHAINTRACEGIINN
Corollary 2.6. Let an be in lemma 2.5
Tmax
max, Ti and
(Tpex/"))
When & I co. the number of Markov chain steps of SUPERCHAINTRACEGIBBS is dominated by ((Emaax)
expl (uil) ExplAHmaat
PARALLEITRACEGIBES For i = -- RELMEANEST (GHH,, J8,&. 11) and RELMEANEST(SH.. . - 1) independently. We show the com- (G
111, 1, PARALLELTRACEGIBES (alg, 2 right) runs
putational complexity of PARALLELTRACHGIBBB in thm.
For 1,2,..... ...,, l, assume Ti is the true (unknown) relaxation time of 9H.s, and T, a known bound on it. For simplicity of presentation we use the following notations: relR, = Rf Range (gageegga_11 .9.) /vv, where μ4 B(fa,....) and v; B(gg,,9..1).
Theorem 2.7 (Efficiency of AARALIEITRACRGGIES)) With probability least 6, holds that the total number m of Markov chain steps taken by PARALLELTRACEGiBSS is upper-bounded by
L Ti relR,
δ log
Σ 2=1
(Reltrven. (fA...)) Re.t (93)-1.8.)
22
Furthermore, for all 1 S S e, Rango(fa. $ /l/ gl/log(n) and Ran s goo(t)/ log", where ci((1) (Hmns/28(H())) for x ~ TT.
PARALLEETRRACEGIIEE and SUPER tradeoffs. PARALLELTRACEGIBEE is usually slower than SUPERCHAINTRACEGIIN because in each iteration i = 1,2,.... l, the mean estimator must acquire a higher-precision estimate so that all estimators together achieve s-ô relative-crror guarantee. Relaxation times (true values and their upper-boundsy appear in a sum in the complexity of PARALLELTTRACHEIIBS, whereas they appear in maximum in CEEEECCHAAARRCEEE of the trace variance terms in both of these algorithms occur different values of e. A comparison of the complexity of these algorithms, in the high-precision regime, with Kolmogorov's TPA + PPE (which uses MCMC black box) presented in table
make different computational complexity
(E 2i=1" Vs. max;= Tς). Furthermore, dominance
Experimental Results
In this section 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 [2
Setup. We the experiments using the single site Gibbs sampler (known also the Glauber dynamics) on two different factor graph models:
(A) The Ising model on 2D lattices. Having 2-dimension lattice of size n, the Hamiltonian defined on n² random variables having values +1 and their dependency is represented by the
Hamiltonian: H (z) -E (xx(i) x(j)). We run the algorithms on lattices of sizes 2 2.

---

3 x 3, 4 x 4, and 6 x 6. For each lattice, the parameter 3 2 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. and A.6. of supplementary material).
(B) The logical voting model. For parameter n, have 2n random variables: the query variable Q E {-1,11 1,, and the voter variables T2, Tn and F1, F P' all in {0, 1). The factors have 2n 1 weights,
:,....Ih The Hamiltonian is:
HQQ,T, = w@wmxmmmmmxmma max Ti - wQ max Fi +
Σ WT. T + Σ
WP FF where w, WF, E [-1,1]
The parameters are reported in fig. We follow De Sa et al. 17 and use hierarchy width to derive upper bounds on mixing times. To make 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. line 8). All code is available at httpt:ts////// thub. com/
eyoopiia/Doublyy_Aatttttt..t.. wwtplll_Adappteeeeeeee
zysophia/Doubly
10
algorithm
algorithm
algorithm
Kolmogorov
Kolmogorov
Kolmogorov
10"
paralleiGibbs
12
paralleiGibbs
peralleiGibbs
superGitbs
superGitbs
absolut: rezttve CITOP
superGibbs
complolly
1.0
25
16'
10'
10'
0.00
1/8
1/6
(a) .01,3x3 3x3 lattice
(b) 8 .0000, lattice
relative errors
Figure 1: Comparison of sample complexity and precision on Ising models. See also the A.6. of the supplementary material
Our experiments demonstrate the practical advantages of our doubly adaptive method, validating our theoretical analysis.
Results:
(1) We first compare the complexity of our algorithms to Kolmogorov's algorithm. Our experiments show the superiority of both versions of our methods different models and various sets of parameters. Figure [ 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(3) varying (fig. (2c,, and when Z is fixed and is varying (fig. 2ap. All of these hold while the precision of our algorithms beats 42 as
(fig. [c].
(2) To demonstrate the advantage of using the relative trace variance, in contrast to the relative variance, run both of our algorithms using simpler sampling, and we compare the results. This is done by setting T <- 1 in of RELMEANEST. In and
estimator which only uses progressive
Figure 2b show the effectiveness of trace averaging, since both SUPERCHAINTRACHCOH C.RALLEI.RAACGIIE beat their simplified versions (T +- 1) after 1/2 passes certain threshold. This is consistent for different parameters of the voting model.
(3) Comparing the performance of SUPERCHAINTRACEEIBBS and PARALLELTRACEGIBBS, we observe that in all of experiments SUPERCHAINTRACEGIII has better performance than PARALLELTRACEGIBBS In fig. 2bb we show the trace variance term PARALLELTRACEGIURT becomes dominant earlier as 1/c grows, thus it performs better in this perspective. This is consistent with our theoretical findings, because the ranges of estimators in PARALLEITRRAEGIIS are smaller than the ranges used in

---

algorithm
Kolmogonox
paraililiibos
3rper
algoniti
pana Koa
Kolmagorov
parale bes
sperGitta
superGibbs
SuperOitts
42*10" 60x13" 90"
45*"0" 6410"
12
(b) Comparison of our algorithms and the effect of trace variance.
(c) Complexity 2(8): 0.025 comparison against TPA PPE 42
Complexity comparison against TPA PPE (42).
Figure 2: Experiments voting models. In (a) and (b) the parameters / 0.1, w 0.9, wr (0.2, 0.5, 0.1) and wr (0.8, 0.2, 0.9). (c), we have n and weights and are picked randomly to generate models with various values ((3).
4 Conclusions: advantages and limitations of proposed algorithms
We develop a doubly-adaptive MCMC-based estimator for the partition function of Gibbs distributions, which resolves major impediment of prior methods that MCMC 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 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 weak upper-bounds on mixing and relaxation times.
In particular, Kolmogorov's 's method requires (+/=²) approximately independent samples, where f is the length of cooling schedule. This requires tight convergence (total variance distance of (s*//) from stationary) for each sample, which adds multiplicative In _. with f = (In Q In Hmax), to its complexity (see column 3 of table and [42), theorem 9). In contrast, our doubly adaptive method only depends relaxation times, which do not depend
Limitations. While significantly improving the state of the art, our methods suffer from several limitations. In SUPERCHAIN TRACEEGGIS the major limitation is the dependence the relative ranges of and G, which can be large, especially when the Hamiltonian range large. Another issue that the product chain's 's mixing time is dominated by E max [rax(T=) as opposed to E т;- While PARALLEEI TRW.EI.RAACEEEER TARAI.EE.TRRAFFOEO ACEGIBBS circumvents these issues by estimating each factor of the telescoping product independently, it fails to beat bound and the nigher-precision guarantees required for each subproblem. Improving performance further will likely require new estimators with smaller ranges and relative trace variances.
efficiency in general, due both to the union
Statement of Broader Impact. While probabilistic graphical models other machine learning methods that rely on MCMC estimations continue to grow in importance and popularity. But rumning the MCMC to theoretical convergence guarantees is often prohibitively expensive, while running 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-focusee maintaining theoretical guarantees of accuracy, with adaptive methods that bound efficiency in terms of (potentially unknown) convergence rate metrics and variances.
practical utility, we seek to reduce proof-burdee while
Acknowledgements. Shahrzad Haddadan is supported by NSF Award CCF-1740741 Cyrus Cousins and Eli Upfal are supported by NSF grant RI- RI-1813444 and DARPA/AFRL grant FA8750. The authors are thankful to anonymous reviewers of NeurIPS 2021 for several valuable inputs.
10

---

References
[1] Permit allocation in emissions trading using the Boltzmann distribution. Physica, 391:4883 4890, 2012.
[2] H. Afshar, S. Sanner, and Christfried Webers. Closed-form Gibbs sampling for graphical models with algebraic constraints, In AAAI, 2016.
[3] David Aldous, Geoffrey R Grimmett, Douglas Howard, Fabio Martinelli, Michael Stecle, and Laurent Aaloff-Coste. Probability on discrete structures, volume 110. Springer Science & Business Media, 2013.
[4] Y. Alimohammadi, Nima Anari, Kirankumar Shiragur, and Vuong. Fractionally log-concave and sector-stable polynomials: Counting planar matchings and more. ArXiv, abs/2102.02708 2021.
Nima Anari, Kuikui Liu, and Shayan Oveis Gharan. Spectral independence high-dimensional expanders and applications to the hardcore model. SIAM Journal on Computing, (0)FFOCS20-1, 2021.
[6] I. Bezáková, Daniel Stefankovic, V. Vazirani, and Eric Vigoda. Accelerating simulated annealing for the permanent and combinatorial counting problems. In SODA 2006, 2006.
[7] Nayantara Bhatnagar, Allan Sly. and Prasad Tetali. Reconstruction threshold for the hardcore model. In Approximation, Randomization, and Combinatorial Optimization. Algorithums and Techniques, pages 434-447. Springer, 2010.
[8] Nayantara Bhainagar, Andrej Bogdanov, and Elchanan Mossel. The computational complexity of estimating MCMC convergence time. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 424 435. Springer, 2011.
[9] A. Blanca, P. Caputo, Z. Chen, D. Parisi, Daniel Stefankovic, and Eric Vigoda. On mixing of Markov chains: Coupling, spectral independence, and entropy factorization. ArXiv, abs/2103.07459, 2021.
[10] Stephen P Brooks and Andrew Gelman. General methods for monitoring convergence of iterative simula- tions. Journal of computational and graphical statistics, 7(4):434 455, 1998.
[11] Stephen Brooks and Gareth O Roberts. Assessing convergence of Markov chain Monte Carlo algorithms.
Statistics and Computing, 8(4):319-335, 1998.
[12] Zongchen Chen, Andreas Galanis, Daniel Štefankovič, and Eric Vigoda. Rapid mixing for colorings via spectral independence. In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1548-1557. STAM, 2021.
[13] H. Cheng. L. Qu, D. Garrick, and R. Fernando. A fast and efficient Gibbs sampler for BayesB in whole-genome analyses. Genetics, Selection, Evolution : GSE, 47, 2015.
[14] Kai-Min Chung, Henry Lam, Zhemming Liu, and Michael Mitzenmacher. Chemoff-HHoeflding bounds for Markov chains: Generalized and simplified. arXiv: 22159 2012.
[15] Barry A Cipra. An introduction to the Ising model. The American Mathematical Monthly, 94(10):937-959 1987.
[16] Cyrus Cousins, Shahrzad Haddadan, and Eli Upfal. Making mcan-estimation more efficient using an MCMC trace variance approach: DynaMITE. CoRR, hts/2011..11123 2020. http:tttts://axiiv..o org/
aba/2011 11129
[17] C. De Sa, Ce Zhang, K. Olukotun, and C. Ré. Rapidly mixing Gibbs sampling for class of factor graphs using hierarchy width. Advances in neural information processing systems, 04:3079-3083 2015.
[18] Chris De Sa, Vincent Chen, and Wing Wong. Minibatch Gibbs sampling on large graphical models. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 1165-1173 PMLR, 10-15 Jul 2018. URL URtphtttp:/prroeeiiii .ttll [19] Christopher De Sa, Kunle Olukotun, and Christopher Ré. Ensuring rapid mixing and low bias for asynchronous Gibbs sampling. In Proceedingy of the 33rd International Conference International Conference on Machine Learning Volume 48, 1CML'16. 16, pages 1567-1576, 2016.
mlr. wwwppeeeee//eeeeeeeeeee
[20] Anand Dixit and Vivekananda Roy. MCMC diagnostics for higher dimensions using Kullback Leibler divergence. Journal of Statistical Computation and Simulaation. 87 (711)22222- 2017. doi: 10.1 [nttps/ts:////oo org/10 http://0.11000004400000 1080/00949655.....o 2017. 1335313
[21] H. Elliott, H. Derin, R. Cristi, and D. Geman. Application of the Gibbs distribution image segmentation.
9:678-681, 1984. doi: 10.1 1199/11///S//... 109/CASSPP.994.111 172637.
[22] G. S. Fishman. Choosing sample path length and number of sample paths when starting in steady state. Oper. Res. Lett., 16:209-211 1994.
[23] Stuart Geman and Donald Geman. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6(6):77) 6(6): 741, 1984. doi: 1ww.1009TAAMM.1844...77.77

---

[24] Charles R Gibbs. Characterization and application of ferrozine iron reagent as ferrous iron indicator.
Analytical Chemistry, 48(8):1197-1201, 1976.
[25] Josiah Willard Gibbs. Elementary Principles in Statistical Mechanics. Scribner, 1902.
[26] Joseph Gonzalez, Yucheng Low, Arthur Gretton, and Carlos Guestrin. Parallel Gibbs sampling: From colored fields to thin junction trees. In Geoffrey Gordon, David Dunson, and Miroslav Dudik, editors. Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, volume 15 of Proceedings of Machine Learning Research, pages 324-332, Fort Lauderdale, FL, USA, 11-13 Apr 2011. PMLR. URL hattt::/pprccedinng mlr. wwwpppeeeeee. htt///11/////// wwwa///nnnalela...oaaaoee [omttl
[27] Tom Griffiths. Gibbs sampling in the generative model latent Dirichlet allocation, Technical report, 2002.
[28] Shahrzad Haddadan and Peter Winkler. Mixing of permutations by biased transpositions. Theory of Computing Systems, 63(5):1068-000 2019.
[29] David G. Harris and Vladimir Kolmogorov. Parameter estimation for Gibbs distributions. CoRR, 0os/2007.10824, 2020. (Rtthhtss://arxt www.oorr/////ooooooooeeoo abs/2007 10824
[30] Bryan He, Christopher De Sa, loannis Mitliagkas, and Christopher Ré. Scan order in Gibbs sampling: Models in which it matters and bounds on how much. Advances in neural information processing sysiems, 29, 2016.
[31] Arnim Hellweg and Frank Eckert. Brick by brick computation of the Gibbs free energy of reaction in solution using quantum chemistry and COSMO-RS. AIChE Journal, 63(0):3944-3954, 2017.
[32] Tomas Hrycej. Gibbs sampling in Bayesian networks. Artificial Intelligence, 46(3):31--33, 1990. ISSN 0004-3702. doi: http:///aooo///rrtcc/////////000
URL httt:p::/www.scieeeee
[33] Daniel Hsu. Aryeh Kontorovich, David A. Levin, Yuval Peres, Csaba Szepesvári, and Geoffrey Wolfer. Mixing time estimation reversible Markov chains from single sample path. The Annals of Applied Probability, 29(4):2439 2480, 2019.
doi: 10. 10.121//11-1 AAA14457 URL lhttts:///do org/10. 1214/
18- AAP1457
[34] Mark Huber. Approximation algorithms for the normalizing constant Gibbs distributions. The Annals Applied Probability, 25(2):974-985, 2015. ISSN 10505164. URL http:://uu. jator /orge.ort//eeeeeeeeeeeeeeeoeeeee org/ stable/
24519939
[35] Mark Huber and Sarah Schott. Random construction of interpolating sets for high-dimensional integration. Journal of Applied Probability 51(1):92-105 100) 2014. doi: 10. (0.1239/ja//113577114 1416.
[36] Mark Huber, Sarah Schott, et al. Using TPA for Bayesian inference. Bayesian Statistics, 9:257-282, 2010.
[37] Mark Jerrum, Leslie G. Valiant, and Vijay V. Vazirani. Random generation of combinatorial structures from a uniform distribution. Theoretical Computer Science, 13:169-188, 1986.
[38] Bai Jiang, Qiang Sun, and Jianqing Fan. Bernstei''s incquality for general Markov chains. arXiv: 10805.007), 2018.
[39] Georgios Karagiannis and Christophe Andrieu. Annealed importance sampling reversible jump MCMC algorithms. Journal of Computational and Graphical Statistics, 22(3):623-648, 2013.
[40] John G Kemeny, Laurie Snell, and Anthony Knapp. Denumerable Markov chains: with a chapter of Markov random fields by David Griffeath, volume 40. Springer Science & Business Media, 2012.
[41] Daphne Koller and Nir Friedman. Probabilistic Graphicul Models: Principles and Techniques Adaptive Computation and Machine Learning. The MIT Press, 2009. ISBN 0262013193.
[42] Vladimir Kolmogorov. A faster approximation algorithm for the Gibbs partition function. In Conference On Learning Theory, pages 228 249. PMLR, 2018
[43] Oswin Krause, Asja Fischer, and Christian Igel. Algorithms for estimating the partition function of restricted Boltzmann machines. Arrificial Intelligence, 278: 103-195, 10 2019. doi: 10. 10.1016/j.arimn artint. 103195.
[44] Patricio S La Rosa, Terrence Brooks, Elena Deych, Berkley Shands, Fred Prior, Linda Larson Prior,, and William D Shannon. Gibbs distribution for statistical analysis of graphical data with a sample application to femri brain images. Statistics in medicine, 35(4):566-580 February 2016. ISSN 0277-67 15. doi: 10. 00.1002/sim.675 6..57.. htRLshttts::/// org/ 10. 1002/s0m..m.mmm 00022/sim6675.. 6757
[45] David A Levin and Yuval Peres. Markov chains and mixing times, volume 107. American Mathematical Soc., 2017.
[46] Pascal Lezaud. Chernoff-type bound for finite Markov chains. The Annals of Applied Probability, 8, 08 1998. doi: 10. 121 12144aap/102889045333
[47] Xianghang Liu and Justin Domke. Projecting Markov random field parameters for fast mixing. NIPS' 14, pages 1377-1385, Cambridge. MA, USA, 2014. MIT Press.
12

---

[48] D. Lunn, D. Spiegelhalter, A. Thomas. and N. Best. The BUGS project: Evolution, critique and future directions. Stat Med., 28(25):3049-67, 2009 Nov 10. doi: 10.1002/sim.3680.
[49] Andrew McCallum, Karl Schultz, and Sameer Singh. FACTORIE: Probabilistic programming via imperatively defined factor graphs. In Advances in Neural Information Processing Systems, vol- ume 22. Curran Associates, Inc., 2009. URLLtttp://proceeddii heweoeeednnngsnuurssss........ http:/pp..cpppppe/2/////e// 2009//11e/////////////
847 https//00ooeeeeeeeeoeaeeee-----e
[50] Michael Mitzenmacher and Eli Upfal. Probability and computing: Randomization and probabilistic techniques in algorithmss and data analysis. Cambridge university press, 2017.
[51] Joris M. Mooij and Soon Ong. libDAI: A free/open source C++ library for discrete approximate inference methods, 2008
[52]
Radford M Neal. Annealed importance sampling. Statistics and computing, 11(2):12--139, 2001.
[53]
David Newman, Padhraic Smyth, Max Welling, and Arthur Asuncion. Distributed inference for latent Dirichlet allocation. In Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2008 URL attppa:/pproceedinge wtw.cceddnnn....uriss..... httprips.cccppppee2/00/
[54]
Marco Patriarca, Anirban Chakraborti, and Kimmo Kaski. Gibbs versus non-gibbs distributions in money
dynamics. Physica A: Statistical Mechanics and its Applications. 340(1):334 339, 2004, ISSN 0378-4371 doi: dott:ttt:://000o0000000011...ppp
04.0024. URL (ttL::/ts:://// / www.sccddreect............. com/stt.oomcccc scie/ce/// att.cle//i1/0077771100005333 News and Expectations in Thermostatistics.
[55] Daniel Paulin. Concentration incqualities for Markov chains by Marton couplings and spectral methods. Electron. J. Probab., 20, 2015.
[56] M. Plummer. JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. 2003.
[57] Ian Porteous, David Newman, Alexander Ihler. Arthur Asuncion, Padhraic Smyth, and Max Welling. Fast collapsed Gibbs sampling for latent Dirichlet allocation. Proceedings of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD '08, pages 569-577, New York, NY, USA, 2008. Association for Computing Machinery. ISBN 978 188160058199 1934. doi: 10.1 01144/1//18899 1401960. URL /RLhttps::/ ftopss://doi. org///......//... 1145/1401890 1401960
[58] Adarsh Prasad, Vishwak Srinivasan, Sivaraman Balakrishnan, and Pradeep Ravikumar. On learning Ising models under Huber's is contamination model. Advances in neural information processing systems, 33, 2020.
[59] Maxim Rabinovich, Aaditya Ramdas. Michael Jordan, and Martin Wainwright. sunction-speciite mixing times and concentration away from equilibrium. Bayesian Analysis, 15, 05 2016. doi: 10. 2214/1-BB1155
[60] Gareth O Roberts and Jeffrey Rosenthal. Examples of adaptive MCMC. Journal of computational and graphical statistics, 18(2):399-377, 2009.
[61] Jeffrey Rosenthal et al. Quantitative convergence rates of Markov chains: A simple account. Electronic Communications in Probability,
2002.
[62] Alexander Smola and Shravan Narayanamurthy. An architecture for parallel topic models. Proceedings of the VLDB Endowment, 3(1-2) September 2010. ISSN 2150-8097. doi: 1920931. URt https:://ast org. 10. 144778/19208// 1920931
[63] D. Stefankovic, S. Vempala, and E. Vigoda. Adaptive simulated annealing: A ncar-optimal connection between sampling and counting. In 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS'07), pages 183-193, 2007. doi: 10.1109/OOSS.200
[64] Andreas Stordal and Ahmed Elsheikh. Iterative ensemble smoothers in the annealed importance sampling framework. Advances in Water Resources, 86:231-239, 2015.
[65] Lucas Theis, Jascha Sohl-Dickstein, and Mathias Bethge. Training sparse natural image models with fast Gibbs sampler of an extended state space. In Proceedings of the 25th International Conference on Neural Information Processing Systems- Volume 1, NIPS' 12, pages 1124-1132, Red Hook, NY, USA,
2012. Curran Associates Inc.
[66] Christopher Tosh. Mixing rates for the alternating Gibbs sampler over restricted Boltzmann machines and friends. In JCML, 2016.
[67] Daniel Štefankovič, Santosh Vempala, and Eric Vigoda. Adaptive simulated annealing: A near-optimal connection between sampling and counting. ACM, 56(3), May 2009. ISSN 0004-5411, doi: 10.1 145/ 1516512.151655 URL htttl:://ddoooi org/10. 1111111111111 /1145//51111..... 1516520
[68] Ruqi Zhang and C. De Sa. Poisson miiistnniinnbathhiing for Gibbs sampling with convergence rate guarantees. In NeurIPS, 2019.
13


Reviews:
[{"review_id": "zc85nLVmZMC", "reviewer": "Reviewer_YZ7S", "summary": "The authors propose a method for estimating the partition function of Gibbs distributions. The method primarily relies on an adaptive mean estimator, and is shown to have theoretical and practical advantages over previous methods.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The problem of approximating the partition function of discrete probabilistic models is relevant and has a lot of ML applications. This paper mainly focuses on improving previous methods based on paired product estimators [1] by introducing an improved mean estimator that does not require black box (approximately) independent samples. Based on previous work on trace averaging [2], they introduce an adaptive mean estimator, and theoretically show that the resulting method is superior to the state of the art. The final algorithms (Alg. 2 in the paper) and their upper bounds are an interesting contribution, albeit not completely novel, since both adaptive parts (cooling and mean estimation) are largely based on previous work.\n\nI see two major downsides with this paper:\n- The exposition of the main theoretical results in section 2 suffers from lack of clarity. Both sections 2.1 and 2.2 are heavy on notation, and contain very little intuitive explanation.\n- The experimental section only discusses very small toy examples, which makes me wonder whether the proposed algorithms would actually be applicable to any practically relevant ML problem.\n\n&nbsp;\n### References\n[1] Mark Huber and Sarah Schott. Random construction of interpolating sets for high-dimensional integration. Journal of Applied Probability, 2014.\\\n[2] Cyrus Cousins, Shahrzad Haddadan, and Eli Upfal. Making mean-estimation more efficient using an MCMC trace variance approach: DynaMITE. ArXiv, 2020.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "eW1g6jZ-14", "reviewer": "Reviewer_DZFN", "summary": "The authors propose a doubly \"adaptive\" MCMC algorithm to approximate the partition function in gibbs distributions.\nIn related work, the partition function of Gibbs distributions are computed using paired product estimators [1].\nSuch 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.\nThe ratio of these two products is then used to compute the partition function.\n\nIn the current work, the first \"adaptive\" refers to the adaptive TPA cooling schedule to decide $[\\beta_{(\\cdot)}]$.\nThe 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$.\n\n[1] - Huber, Mark. \"Approximation algorithms for the normalizing constant of Gibbs distributions.\" The Annals of Applied Probability 25.2 (2015): 974-985.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### Strengths\n\n- The method is well motivated and the theoretical analysis is very thorough.\n\n### Weaknesses\n\n- The usage of the word adaptive (as an auto-decision on the number of steps) is highly misleading. An entire line of work [2] exists where the word adaptive refers to changing the transition kernel to improve the efficiency of the MCMC.\n- In line 77 describing Huber and Schott's estimator, [1] makes statements on estimating $z=\\log(Z)$. Please refer to the discussion 2 paragraphs before algorithm 2.2 in [1] ([29] in reviewed paper). (Maybe I am missing something but in order to apply linearity of expectations to a product and to compute an unbiased estimator of a quotient, the log domain makes more sense. However please correct me if I am missing something.)\n- Please consider comparison with Annealed Importance Sampling (AIS) [3] which samples a single transition at each temperature such that stationary samples at each temperature aren't needed.\n  - I may have missed a nuance which prohibits the use of AIS, and if that is the case I am willing to change my rating.\n  - If AIS can indeed be used in this scenario, the second \"adaptation\" can only be justified thorough empirical evaluation.\n- Please add error bars to all plots not just Fig 1.e.\n- Please use larger font in plots.\n\n### Other criteria\n\n- Originality : Existing problem, novel solution which combines [1,4] using original tools.\n- Quality : Theory is very rigorous, empirical evaluation lacking.\n- Clarity : Motivations are well explained, plots need work and experimental evaluation is lacking.\n- Significance : Unclear, refer to comment regarding AIS in Weaknesses.\n\n### Note\n\nI am willing to update my review, especially my overall score, based on a strong rebuttal arguing that either AIS is not applicable in this task or is demonstrably bad when empirically evaluated against the method.\n\n[1] - Huber, Mark. \"Approximation algorithms for the normalizing constant of Gibbs distributions.\" The Annals of Applied Probability 25.2 (2015): 974-985.\n\n[2] - Roberts, Gareth O., and Jeffrey S. Rosenthal. \"Examples of adaptive MCMC.\" Journal of computational and graphical statistics 18.2 (2009): 349-367.\n\n[3] - Neal, Radford M. \"Annealed importance sampling.\" Statistics and computing 11.2 (2001): 125-139.\n\n[4] - Cousins, Cyrus, Shahrzad Haddadan, and Eli Upfal. \"Making mean-estimation more efficient using an MCMC trace variance approach: DynaMITE.\" arXiv preprint arXiv:2011.11129 (2020).", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Adequately addressed", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "MbBi90BGAgo", "reviewer": "Reviewer_AZX6", "summary": "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.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the reviewer thinks this is a good paper. As claimed by the authors, they introduce a new method that is proven to outperform the state of art algorithms. And the experiments also justify this. \n\nHowever, the reviewer thinks this paper is not very clear. The reviewer is not familiar with this problem and feels it's a little hard to follow this paper. The reviewer will point out some confusing parts in \"Limitations\".", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Some limitations of this paper are listed as follows: (The reviewer is not an expert in this area and might misunderstand something.)\n1. In the contribution section (Sec 1.2), the reviewer thinks the authors can delete lines 121-129 since all four points explain how the paper achieves an improvement in computational complexity. \nIt's good that the authors explain how they obtain the improvement. But they are actually not the main contributions, the authors should move them somewhere else.\n2. In the conclusion section (Sec 4), Table 1 compares the computation complexity of different methods. But the authors omit the $1/\\epsilon$ term. It's true that this term can be omitted when $\\epsilon$ is small enough, but the authors should mention the conditions that this holds. The authors also mention $RelR_i$ is a large number, then it seems we should choose $\\epsilon$ super small to see the priority of new methods.\n3. According to Table 1, the authors don't give any comparison between $\\text{Reltrv}^{\\tau_i}$ and $V_{rel}$. If they are in the same order, the reviewers think the new method is better because $\\ln(l/\\epsilon)>l^2$ when $\\epsilon$ is small. Then it seems that the main improvement comes from the new bound eliminates the $\\ln(1/\\epsilon)$ factor. This seems a small improvement. Can the authors give some explanation why this is important?\n\nAs a potential reader, the reviewer thinks that there are some confusing places in the paper:\n1. In line 59, \"the values of random terms in the sum have large variance\". According to (1), the reviewer didn't see any randomness. Does the author mean the variance of $\\pi_{\\beta}$ is large so $Z(\\beta)$ is hard to estimate?\n2. In line 83, $\"O(log(Q)log(H_{\\max})\"$, this is the first time $H_{\\max}$ appears. But the authors define it in Lemma 2.6 line 231. The authors can move the definition forward or at least give an explanation in line 83.\n3. In line 84, \"TPA(k,d)\", this is the first time $d$ appears. Again, there is no definition and explanation. \n4. In lines 100 and 108, what is $\\pi_\\min$? The reviewer didn't find its definition in the paper. \n5. In the results of Section 2.2, what is $\\epsilon$? \n6. In line 303, according to Table 1, the reviewers think \"$\\ln(l)$\" should be \"$\\ln(l/\\epsilon)$\".\n7. In line 305 Table 1, what is $T_i$? Is it the same as $\\tau_i$?\n8. As said before, the comparison in Table 1 is very unclear. The authors should spend more space to do the comparison and make it clear.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2 hours"}, "extra_text": {}}, {"review_id": "AyD4Veae_G1", "reviewer": "Reviewer_TAYx", "summary": "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.\nThe 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. ", "questions": "", "limitations": "", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\nMajor comments:\nExcept for the presentation of the basic algorithm (which is fairly classical), the paper is difficult to read. It is difficult to understand the contributions, the results, I was  quite confuse even after reading the full paper. The discussion and the proofs use a lot of external references, which require each time to understand other conventions, notations, concepts. It would be much better to recap the required results in a Section in the supplement paper, to ease the reading. \n\nI have tried to check the proof of Theorem 2.4 in the supplementary version for example, and I found it quite difficut to follow. I would definitely prefer in the supplement paper to get  complete and detailed arguments, not a sketch of proofs with external references (when you cite an external result in a paper is definitely better to get a precise reference to the theorem or proposition which is used, it is very complicated otherwise). \n\nI am surprised not to see in the complexity bounds the initial distribution of the Markov chain. The authors use the results  of Jiang, Sun, Fan (2018) [Theorem 1 and 2] which both assumed that the MC is stationary. It would be interesting to know how the authors plan to work with the non-stationary case.  \n\n\nThe assumption of knowing an upper bound of the spectral gap seems difficult to me. The estimators are, except for examples, generally extremely pessimistic, especially for general state chains. \nMinor comments\n- I find the discussion at the beginning of Section 2 difficult to follow\n- Lemma 1 is dificult to understand because some quantities which are presented are not defined (what is  $\\Delta_i$ ? )\n- line 177: \"estimating means.. union bound).\" difficult to understand as it is stated. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "none", "review_text": "", "extra_scores": {"ethical_concerns": "none", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The goal of this paper is a better algorithm to estimate the partition function / normalization constant in a Gibbs distribution. These are often of interest in many different forms in the NeurIPS community. The best current methods are based on annealing, essentially starting at a normalized distribution and moving through a series of distributions of different temperatures and estimating the ratio of each using MCMC. This paper proposes a new algorithm in this line which is intended to be less sensitive to tightness in the mixing times of MCMC.

This is a highly mathematical paper, where the goal is improved theoretical guarantees (rather than, say, practical performance improvement.) Unfortunately, all reviewers signaled that they were not confidently able to verify the correctness of the core results. In the end, the SAC stepped in and did a careful reading of the paper and recommends acceptance.

The SAC's reasoning is as follows:

The presentation is indeed terse and dense. But that is par for the course for theoretical papers at NeurIPS given the page limit. Everything was well defined (as far as I saw), and all the arguments seemed credible. The authors outlined the strategy of their arguments, hiding the considerable complexity of all the details in the supplementary material. I actually found the writing pretty good compared to the average published Neurips paper.

Regarding the concern by one reviewer about the assumption re spectral gap: I do not see this as a downside. It seems inevitable that one will need some assumptions on the dynamics of these chains in order to say anything theoretical. The goal (in my mind) is not to attempt the impossible task of theoretically explaining what "really" happens in practice, because one never knows. Assuming something about a spectral gap in the transition matrix of a markov chain seems to me to be less fraught epistemologically that (say) assuming one's data is drawn iid from some distribution ... an assumption that I doubt is true of almost all practical applications of machine learning, but upon which almost all ML theory is built... (Thus my view is aligned to that of the authors in their "Note to chair")

So it looks like this is a real advance on a problem of interest. The analysis looks credible (all one can ever really tell without spending a week working through every single step). It is easily at least as clear, correct, and significant as the majority of published Neurips papers so I do not have any qualms recommending acceptance....
Decision: Accept (Poster)