figurec

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

Shahrzad Haddadanⓡ<EMAIL_ADDRESS>Yue Zhuangⓡ
<EMAIL_ADDRESS>Cyrus Cousinsⓡ<EMAIL_ADDRESS>Eli Upfal
<EMAIL_ADDRESS>
###### 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.

##  1 Introduction

The Gibbs (Boltzmann) distribution is a family of probability distributions of
exponential form. First introduced in the context of statistical mechanics
[26], Gibbs distributions are now ubiquitous in a variety of other
disciplines, such as chemistry [32, 25], economics [56, 1] 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 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 [42, 51,
53], including restricted Boltzmann machines [68, 44], Markov random fields
[41, 49], and Bayes networks [33], and are applied in the analysis of images
and graphical data [67, 45, 24, 21], topic modeling (LDA) [28, 59, 55, 64],
and more [13, 58, 2, 18, 17, 70, 27, 31, 50, 60].

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 [54, 38, 23,
69, 65, 35, 6, 43, 30, 40].

Like other MCMC methods, here various heuristics are used. The most well-known
heuristics are the _annealed importance sampling_ [54, 66, 40] or _convergence
diagnostics methods_ [20, 11, 10, 63]. 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 [54, section 4] or [7, 34]).

On the other hand, theoreticians study this problem by designing Fully
Polynomial Randomized Approximation Schemes (FPRAS) (see 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
adaptive111The 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
[35, 43] (see § 1.1), and should not confused with the work of [62].
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 [5, 65, 9,
12, 4, 29]. 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 [15, 5, 8].

###  1.1 Preliminaries and Prior Work

Assume a _sample space_ $\Omega$, _Hamiltonian function_
$H\mathrel{\mathop{\mathchar
58\relax}}\Omega\rightarrow\\{0\\}\cup[1,\infty)$, and _inverse temperature_
parameter $\beta\in\mathbb{R}$, referred to as _inverse temperature_. The
_Gibbs distribution_ on $\Omega$, $H(\cdot)$, and $\beta$ is then
characterized by probability law

$\forall x\in\Omega\mathrel{\mathop{\mathchar 58\relax}}\
\pi_{\beta}(x)\doteq{\frac{1}{Z{({\beta})}}}\exp\bigl{(}-\beta
H(x)\bigr{)}\enspace.$ (1)

Here $Z{({\beta})}$ is the _normalizing constant_ or _Gibbs partition
function_ (GPF) of the distribution, with

$Z{({\beta})}\doteq{\sum_{x\in\Omega}}\exp\bigl{(}-\beta
H(x)\bigr{)}\enspace.$ (2)

Estimating the GPF $Z{({\beta})}$, is computationally challenging, since
typically the size of $\Omega$ 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 $\Omega$, a Hamiltonian function $H$, and a parameter $\beta$,
design a Fully Polynomial Randomized Approximation Scheme (FPRAS) for
estimating the partition function
$Z{({\beta})}\doteq\sum_{x\in\Omega}\exp\bigl{(}-\beta H(x)\bigr{)}$. In other
words, for user-supplied $\varepsilon$, the task is to produce an estimate
$\hat{Z}(\beta)$, such that with probability at least $1-\delta$, we have
$(1-\varepsilon){Z}(\beta)\leq\hat{Z}(\beta)\leq(1+\varepsilon){Z}(\beta)$, in
time polynomial in $\nicefrac{{1}}{{\varepsilon}}$,
$\ln(\nicefrac{{1}}{{\delta}})$, 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 [23, 6, 69, 6], the current state of the
art is due to Huber and Schott [36], with Kolmogorov’s sharper analysis [43].
They introduce the _paired product estimator_ (PPE), see 1.1, and apply the
_tootsie-pop algorithm_ (TPA) to adaptively compute a near-optimal _cooling
schedule_ , i.e., a sequence of inverse temperatures
$\beta_{0}<\beta_{1}<\dots<\beta_{\ell-1}<\beta_{\ell}$ satisfying
$\beta_{\ell}=\beta$, and that $Z{({\beta_{0}})}$ is easy to compute, e.g.,
$\beta_{0}=0$ is often convenient, since
$Z{({0})}=\mathinner{\\!\left\lvert\Omega\right\rvert}$. We thus define
$Q\doteq\nicefrac{{Z{({\beta})}}}{{Z{({\beta_{0}})}}}$ and estimate it using
the paired product estimator.

######  Definition 1.1 (PPE [35]).

Assume a _cooling schedule_ $\beta_{0},\beta_{1},\dots,\beta_{\ell}$. For each
pair $(\beta_{i},\beta_{i+1})$ in the schedule, we define two random
variables, $X_{i}\sim\pi_{\beta_{i}}$ and $Y_{i}\sim\pi_{\beta_{i+1}}$, all
independent, and we then define
$f_{\beta_{i},\beta_{i+1}}\doteq\exp\bigl{(}-\frac{\beta_{i+1}-\beta_{i}}{2}H(X_{i})\bigr{)}$
and
$g_{\beta_{i},\beta_{i+1}}\doteq\exp\bigl{(}\frac{\beta_{i+1}-\beta_{i}}{2}H(Y_{i})\bigr{)}$.
It is easy to verify that
$\mathbb{E}[f_{\beta_{i},\beta_{i+1}}]=Z{({\frac{\beta_{i}+\beta_{i+1}}{2}})}/Z{({\beta_{i}})}$,
and
$\mathbb{E}[g_{\beta_{i},\beta_{i+1}}]=Z{({\frac{\beta_{i}+\beta_{i+1}}{2}})}/Z{({\beta_{i+1}})}$.
We then define $F\doteq\prod_{i=1}^{k}f_{\beta_{i},\beta_{i+1}}$,
$G\doteq\prod_{i=1}^{k}g_{\beta_{i},\beta_{i+1}}$. Letting $\hat{\mu}$ and
$\hat{\nu}$ denote empirical estimates of $\mathbb{E}[F]$ and $\mathbb{E}[G]$,
respectively, the _paired product estimator_ (PPE) is
$\hat{Q}\doteq\nicefrac{{\hat{\mu}}}{{\hat{\nu}}}$ .

Denote by ${\mathbb{V}}_{\rm
rel}[X]\doteq\nicefrac{{\mathbb{E}[X^{2}]}}{{\mathbb{E}[X]^{2}}}-1=\nicefrac{{{\mathbb{V}}[X]}}{{\mathbb{E}[X]^{2}}}$
the _relative variance_ of a random variable $X$. _The TPA schedule_ [37, 36]
is generated by an adaptive algorithm, which, by a proper setting of
parameters, outputs a cooling schedule guaranteeing constant $\mathbb{V}_{\rm
rel}[F]$ and $\mathbb{V}_{\rm rel}[G]$ (see alg. 3 in the supplementary
material). Kolmogorov [43] 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 ([43]).

Let $H_{\max}\doteq\max_{x\in\Omega}H(x)$, using  Tpa $(k,d)$, $k=\Theta(\log
H_{\max})$ and $d=16$ to generate cooling schedule
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})$. W.h.p., we have
$\ell=\Theta(\log(Q)\log(H_{\max}))$ and $\mathbb{V}_{\rm
rel}[F]+1=\prod_{i=1}^{\ell}({\mathbb{V}}_{\rm
rel}[f_{\beta_{i},\beta_{i+1}}]+1)=\Theta(1)$ and $\mathbb{V}_{\rm
rel}[G]+1=\prod_{i=1}^{\ell}({\mathbb{V}}_{\rm
rel}[g_{\beta_{i},\beta_{i+1}}]+1)=\Theta(1)$.

Kolmogorov [43] 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 [36] assume unit-cost for exact sampling from each
$\pi_{\beta_{i}}$, and Kolmogorov [43] 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.

Let $\cal M$ be an ergodic Markov chain with state space $S$ and stationary
distribution $\pi$. Let $\tau_{\rm mix}(\varepsilon)$ denote the
$\varepsilon$-_mixing time_ of $\cal M$, and define $\tau_{\rm
mix}\doteq\tau_{\rm mix}(\nicefrac{{1}}{{4}})$. Letting $\lambda$ denote the
_second largest absolute eigenvalue_ of $\mathcal{M}$’s transition matrix, the
_relaxation time_ of $\cal M$ is $\tau_{\rm rx}\doteq(1-\lambda)^{-1}$, and it
is related to the mixing time $\tau_{\rm mix}$, by $\left(\tau_{\rm
rx}(\mathcal{M})-1\right)\ln(2)\leq\tau_{\rm
mix}(\mathcal{M})\leq\bigl{\lceil}\tau_{\rm
rx}(\mathcal{M})\ln\bigl{(}\smash{\nicefrac{{2}}{{\sqrt{\pi_{\min}}}}}\bigr{)}\bigr{\rceil}$
[47]. Let $T$ be an upper bound on $\max\\{\tau_{\rm
rx}(\mathcal{M}),\tau_{\rm mix}(\mathcal{M})\\}$.

Consider any i.i.d. sampling concentration bound like Chebyshev’s,
Hoeffding’s, or Bernstein’s inequalities [52], with, say, sample complexity
$m_{\varepsilon}$. Using MCMC as a black-box sampling tool, we obtain the same
precision estimation guarantees, with a computational cost of
$m_{\varepsilon}\cdot\tau_{\rm mix}(\varepsilon/m_{\varepsilon})$, which is
equal to $m_{\varepsilon}\log(m_{\varepsilon}\cdot\varepsilon^{-1})\cdot T$ in
the absence of exact values for $\tau_{\rm 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_ [57, 52, 48, 14, 39], or function specific mixing time [61].
Note that since $\log(\frac{1}{2\varepsilon})(\tau_{\rm rx}-1)\leq\tau_{\rm
mix}(\varepsilon)\leq\log(\frac{1}{\varepsilon\pi_{\min}})\tau_{\rm rx}$,
using these bounds is often more efficient, saving at least
$\log(m_{\varepsilon})$ 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 1,
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 _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.

##  2 Algorithms

In this section, we develop two _doubly-adaptive_ _fully polynomial randomized
approximation schemes_ providing more efficient algorithmic solutions to 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(\tau_{\rm
mix},\tau_{\rm rx})$, and $\Lambda$ denotes a upper-bound on the second
absolute eigenvalue $\lambda$. We use ${\mathcal{G}}_{H,\beta}$ to denote any
Markov chain with Gibbs stationary distribution $\pi_{\beta}$, eq. (1). Having
the Hamiltonian $H$, we denote its maximum and minimum values as $H_{\rm max}$
and $H_{\rm min}$, i.e., $H_{\rm max}\doteq\max_{x\in\Omega}\\{H(x)\\}$ and
$H_{\rm min}\doteq\min_{x\in\Omega}\\{H(x)\\}$. Having a schedule
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})$, the paired product estimators
$f_{\beta_{i},\beta_{i+1}}$, $g_{\beta_{i},\beta_{i+1}}$,
$F=\bigotimes_{i=1}^{\ell}f_{\beta_{i},\beta_{i+1}}$ and
$G=\bigotimes_{i=1}^{\ell}g_{\beta_{i},\beta_{i+1}}$ are as in 1.1. When
writing
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})=\scalebox{0.92}[0.98]{{Tpa}}(k,d)$,
we mean the cooling schedule is obtained from running alg. 4 in the Appendix,
and we always set $k=\log H_{\rm max}$ and $d=64$, as these parameters are
shown to produce a near-optimal schedule w.h.p. [43].

We first introduce a novel MCMC-based _multiplicative_ mean estimation
procedure  RelMeanEst (see alg. 3), and analyze its computational complexity
in terms of a new quantity, which we coin _the relative trace variance_ (see
2.1).  RelMeanEst receives as input a Markov chain $\mathcal{M}$, a function
$f$, and precision parameters $\varepsilon$ and $\delta$, 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 $(\mathcal{M},f)$, leaving out the precision parameters.

Letting
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})=\scalebox{0.92}[0.98]{{Tpa}}(k,d)$,
we first present  ParallelTraceGibbs , in which we invoke both  RelMeanEst
$({\mathcal{G}}_{H,\beta_{i}},f_{\beta_{i},\beta_{i+1}})$ and  RelMeanEst
$({\mathcal{G}}_{H,\beta_{i}},g_{\beta_{i},\beta_{i+1}})$ for each
$i=1,2,\dots,\ell-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 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 [43] (thm. 2.4 and 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
$\varepsilon\geq\varepsilon_{0}$ their computation complexity is dominated by
the _true relaxation time_ $\tau_{\rm rel}$ (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 trace
variance_ depends on both the chain $\mathcal{M}$ and the function $f$, and it
generalizes the _relative variance_ , defined as ${\mathbb{V}}_{\rm
rel}[f]\doteq\nicefrac{{{\mathbb{V}}[f]}}{{\mathbb{E}[f]^{2}}}$, which depends
only on $f$, and is used in i.i.d. regimes.

######  Definition 2.1 (Relative Trace Variance).

For arbitrary $\tau$, consider a trace of length $\tau$ of a Markov chain
$\cal M$, and a real-valued function $f$. On $\cal M$, we define the relative
trace variance of $f$ as

$\mathrm{Reltrv}^{\tau}_{{\cal
M}}[f]\doteq\frac{\mathbb{E}[{\bar{f}(\vec{X}_{1\mathrel{\mathop{\mathchar
58\relax}}\tau})}^{2}]}{(\mathbb{E}[\bar{f}(\vec{X}_{1\mathrel{\mathop{\mathchar
58\relax}}\tau})])^{2}}-1\enspace,$

where $\vec{X}_{1\mathrel{\mathop{\mathchar 58\relax}}\tau}\doteq
X_{1},X_{2},\dots,X_{\tau}$ is a trace of length $\tau$ of $\cal M$, and
$\bar{f}(\vec{X}_{1\mathrel{\mathop{\mathchar
58\relax}}\tau})\doteq(\frac{1}{\tau})\sum_{i=1}^{\tau}f(X_{i})$. 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 ${\rm trv}^{(\tau)}(\mathcal{M},f)$, which they
showed it captures MCMC-mean estimation with _additive_ precision guarantees
[16]. In fact, the two terms are related as

$\mathrm{Reltrv}^{\tau}_{\cal M}[f]=\frac{{\rm
trv}^{(\tau)}(\mathcal{M},f)}{(\mathbb{E}[\bar{f}(\vec{X}_{1\mathrel{\mathop{\mathchar
58\relax}}\tau})])^{2}}\enspace.$

Note that the two terms are not easily convertible without knowing the _mean_
, $\mathbb{E}[\bar{f}(\vec{X}_{1\mathrel{\mathop{\mathchar 58\relax}}\tau})]$.

######  Lemma 2.1.

For any $\tau$ we have

$\mathrm{Reltrv}^{\tau}_{\cal M}[f]\leq\mathbb{V}_{\rm rel}[f]\enspace.$ (3)

Furthermore, for $\tau\geq\tau_{\rm rx}({\cal M})$ we have,

$\mathrm{Reltrv}^{\tau}_{\cal M}[f]=O\left(\frac{\tau_{\rm rx}({\cal
M})}{\tau}\mathrm{Reltrv}^{\tau_{\rm rx}({\cal M})}_{\cal
M}[f]\right)\enspace.$ (4)

Lemma 2.1 enables us to compare the computational complexity of our algorithms
with the state of the art [43]. In particular, using (3), we show our results
improve the state of the art (which is in terms of $\mathbb{V}_{\rm rel}$),
and using (4), we show that for high-precision estimations, the sample
complexity of our algorithms only depends on $\tau_{\rm 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 inter-trace 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., [47].

######  Definition 2.2 (Product Chain and Tensor Product Function).

Consider $k$ Markov chains $\\{{\cal M}_{i}\\}_{i=1}^{k}$ each defined on
state space $S_{i}$ and assume real valued functions
$\\{f_{i}\mathrel{\mathop{\mathchar
58\relax}}S_{i}\rightarrow{\mathbb{R}}\\}_{i=1}^{k}$. The _product chain_
${\cal M}^{\otimes}_{1\mathrel{\mathop{\mathchar 58\relax}}k}$ is defined on
the Cartesian product of $S_{i}$ as follows: at any step ${\cal
M}^{\otimes}_{1\mathrel{\mathop{\mathchar 58\relax}}k}$ chooses $i$ with
probability $\omega_{i}$ (thus $\sum_{i=1}^{k}\omega_{i}=1)$, and moves from
$(x_{1},x_{2},\dots,x_{i},\dots,x_{k})$ to
$(x_{1},x_{2},\dots,y_{i},\dots,x_{k})$, with the transition probability of
moving from $x_{i}$ to $y_{i}$ in ${\cal M}_{i}$. The _tensor product_ of
$\\{f_{i}\\}_{i=1}^{k}$, denoted by $\bigotimes_{1\mathrel{\mathop{\mathchar
58\relax}}k}f_{i}$, is defined as
$\bigl{(}\bigotimes_{1\mathrel{\mathop{\mathchar
58\relax}}k}f_{i}\bigr{)}(x_{1},x_{2},\dots,x_{k})=\prod_{i=1}^{k}f_{i}(x_{i})$.

#####  RelMeanEst

Let $T$ denote an upper bound on the relaxation time of a Markov chain
$\mathcal{M}$.  RelMeanEst receives $T$, $\mathcal{M}$, $f$ and precision
parameters $\varepsilon$ and $\delta$ as input. Before it starts collecting
samples, it runs the chain for a _warm start_ (6 of alg. 3). Starting from a
minimum sample size $m^{\downarrow}$, it runs $\mathcal{M}$ for $T\cdot
m^{\downarrow}$ steps, and collect samples $X_{1},X_{2},\dots,X_{T\cdot
m^{\downarrow}}$. It then computes for $j=1,2,\dots,m^{\downarrow}$,
$\bar{f}_{j}\doteq\sum_{i=(j-1)\cdot T+1}^{j\cdot T}f(X_{i})$; using them, it
calculates an empirical estimate of the mean, $\hat{\mu}$, and an empirical
estimation for the _trace variance_ of $\mathcal{M}$ and $f$, $\hat{v}$. Based
on these estimates, we derive an upper-bound on the current trace variance
$u_{i}$ and relative error $\hat{\bm{\varepsilon}}_{i}^{\times}$, and check
whether is smaller than the user-specified error $\varepsilon$ (lines 18-19).
If so, we return the current mean estimate, otherwise we double the sample
size and repeat.

1:procedure

RelMeanEst

2: Input: Markov chain $\mathcal{M}$, upper-bound on relaxation time $T$,
real-valued function $f$ with range $[a,b]$, letting $R=b-a$, multiplicative
precision $\varepsilon$, error probability $\delta$.

3: Output: Multiplicative approximation $\hat{\mu}$ of
$\mu=\operatorname*{\mathbb{E}}_{\pi}[f]$.

4: $T\leftarrow\left\lceil\frac{1+\Lambda}{1-\Lambda}\ln\sqrt{2}\right\rceil$;
$\Lambda^{\prime}\leftarrow\Lambda^{T}$ $\triangleright$ Choose $T$ to be an
upperbound on relaxation time

5: $I\leftarrow
1\vee\left\lfloor\log_{2}\left(\frac{bR}{2a^{2}}\cdot\frac{(1-\varepsilon)^{2}}{(1+\varepsilon)\varepsilon}\right)\right\rfloor;\
\alpha\leftarrow\frac{(1+\Lambda^{\prime})R\ln\frac{3I}{\delta}(1+\varepsilon)}{(1-\Lambda^{\prime})b\varepsilon};\
m_{0}\leftarrow 0$ $\triangleright$ Initialize _sampling schedule_

6: $T_{\rm unif}\leftarrow\left\lceil T\cdot\ln(\nicefrac{{1}}{{\pi_{\rm
min}}})\right\rceil$;
$(\vec{X}_{0,1},\vec{X}_{0,2})\leftarrow\mathcal{M}^{T_{\rm unif}}(\bot)$
$\triangleright$ Warm-start two chains for $T_{\rm unif}$ steps from arbitrary
$\bot\in\Omega$

7: for $i\in 1,2,\dots,I$ do

8: $m_{i}\leftarrow\left\lceil\alpha r^{i}\right\rceil$ $\triangleright$ Total
sample count at iteration $i$; $r$ is the geometric ratio (constant, usually
2) size

9: for $j\in(m_{i-1}+1),\dots,m_{i}$ do

10: $(\vec{X}_{j,1},\vec{X}_{j,2})\leftarrow(T$ steps of $\mathcal{M}$
starting at $\vec{X}_{j-1,1},\vec{X}_{j-1,2})$ $\triangleright$ Run two
independent copies of $\mathcal{M}$ for $T$ steps

11:
$\displaystyle\bar{f}(\vec{X}_{j,1})\leftarrow\frac{1}{T}\sum_{t=1}^{T}f\bigl{(}\vec{X}_{j,1}(t)\bigr{)}$;
$\displaystyle\bar{f}(\vec{X}_{j,2})\leftarrow\frac{1}{T}\sum_{t=1}^{T}f\bigl{(}\vec{X}_{j,2}(t)\bigr{)}$
$\triangleright$ Average $f$ over $T$-traces

12: end for

13:
$\displaystyle{\hat{\bm{\mu}}}_{i}\leftarrow\frac{1}{2m_{i}}\sum_{i=1}^{m_{i}}\bigl{(}f(\vec{X}_{j,1})+f(\vec{X}_{j,2})\bigr{)}$;
$\displaystyle\hat{v}_{i}\leftarrow\frac{{1}}{2{m_{i}}}\sum_{i=1}^{m_{i}}\bigl{(}(f(\vec{X}_{j,1})-f(\vec{X}_{j,2})\bigr{)}^{2}$
$\triangleright$ Compute empirical mean; trace variance

14: $\displaystyle
u_{i}\leftarrow\hat{v}_{i}+\frac{(11+\sqrt{21})(1+\nicefrac{{\Lambda^{\prime}}}{{\sqrt{21}}})R^{2}\ln\frac{3I}{\delta}}{(1-\Lambda^{\prime})m_{i}}+\sqrt{\frac{(1+\Lambda^{\prime})R^{2}\hat{v}_{i}\ln\frac{3I}{\delta}}{(1-\Lambda^{\prime})m_{i}}}$$\triangleright$
Variance upper bound

15:
$\displaystyle\hat{\bm{\varepsilon}}^{+}_{i}\leftarrow\frac{10R\ln\frac{3I}{\delta}}{(1-\Lambda^{\prime})m_{i}}+\sqrt{\frac{(1+\Lambda^{\prime})u_{i}\ln\frac{3I}{\delta}}{(1-\Lambda^{\prime})m_{i}}}$
$\triangleright$ Apply Bernstein bound

16:
$\displaystyle\hat{\bm{\mu}}^{\times}_{i}\leftarrow\frac{(\hat{\bm{\mu}}_{i}-\hat{\bm{\varepsilon}}^{+}_{i})\vee
a+(\hat{\bm{\mu}}_{i}+\hat{\bm{\varepsilon}}^{+}_{i})\wedge b}{2}$
$\triangleright$ Optimal mean estimate

17:
$\displaystyle\hat{\bm{\varepsilon}}^{\times}_{i}\leftarrow\frac{((\hat{\bm{\mu}}_{i}+\hat{\bm{\varepsilon}}^{+}_{i})\wedge
b-(\hat{\bm{\mu}}_{i}-\hat{\bm{\varepsilon}}^{+}_{i})\vee
a}{2\hat{\bm{\mu}}^{\times}_{i}}$ $\triangleright$ Empirical relative error
bound

18: if $(i=I)\vee(\hat{\bm{\varepsilon}}^{\times}_{i}\leq\varepsilon)$ then
$\triangleright$ Terminate if accuracy guarantee is met

19: return $\hat{\bm{\mu}}_{i}^{\times}$

20: end if

21: end for

22:end procedure

Algorithm 1

RelMeanEst

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-\delta$,  RelMeanEst will output $\hat{\mu}$
satisfying $(1-\varepsilon)\hat{\mu}\leq\mu\leq(1+\varepsilon)\hat{\mu}$.
Furthermore, with probability at least $1-\frac{\delta}{3I}$, the total Markov
chain steps of  RelMeanEst , $\hat{m}$, obeys

$\hat{m}\in\mathcal{O}\left(\ln\left(\frac{\ln\frac{b}{a\varepsilon}}{\delta}\right)\left(\frac{T\cdot
R}{\mu\varepsilon}+\frac{\tau_{\rm rx}\mathrm{Reltrv}^{\tau_{\rm
rx}}}{\varepsilon^{2}}\right)\right).$ (5)

###  2.2 Doubly adaptive algorithms:  SuperChainTraceGibbs and
ParallelTraceGibbs

Let
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})=\scalebox{0.92}[0.98]{{Tpa}}(k,d)$,
and consider a family of Gibbs chains ${\mathcal{G}}_{H,\beta_{i}}$, each
corresponding to some $\beta_{i}$, and the paired product estimators
$F=\bigotimes_{i=1}^{\ell}f_{\beta_{i},\beta_{i+1}}$
$G=\bigotimes_{i=1}^{\ell}g_{\beta_{i},\beta_{i+1}}$. The  Tpa method is
designed to ensure $\mathbb{V}_{\rm rel}$ 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\cdot R/\mu$, as the _range term_ , and to the term $\tau_{\rm
rx}\mathrm{Reltrv}^{\tau_{\rm rx}}$ as the _trace variance term_. Note that as
$\varepsilon$ 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 $\tau_{\rm rx}$
and $\mathrm{Reltrv}^{\tau_{\rm rx}}$.

In order to ensure that the ranges of estimators are small, we prove that the
length of each inverse-temperature interval in the  Tpa schedule is w.h.p.
small. Having a schedule $(\beta_{0},\beta_{1},\dots,\beta_{\ell})$ we define
and use the following notation: for $0\leq i\leq\ell-1$, _interval length_
$\Delta_{i}\doteq\beta_{i+1}-\beta_{i}$, _maximum interval length_
$\Delta_{\max}\doteq\max_{i}\Delta_{i}$, and _total length_
$\Delta\doteq\beta_{\ell}-\beta_{0}$ .

######  Lemma 2.3.

Let $z(\beta)\doteq\ln\left(Z{({\beta})}\right)$, and let $\beta_{i}$,
$\beta_{i+1}$ be two consecutive points generated by
$\scalebox{0.92}[0.98]{{Tpa}}(k,d)$. For arbitrary $\varepsilon\geq 0$, we
have:

  1. 1.

$\mathbb{P}(z(\beta_{i})-z(\beta_{i+1})\leq\varepsilon)\geq(1-\exp(-\varepsilon
k/d))^{d}$.

  2. 2.

$\mathbb{P}\left(\Delta_{i}\geq\nicefrac{{\varepsilon}}{{\mathbb{E}[H(x)]}}\right)\leq
d\exp(-\varepsilon k/d),$ where $\mathbb{E}[H(x)]$ is taken w.r.t.
$x\sim\pi_{\beta_{i+1}}$.

#####  SuperChainTraceGibbs

Let $\cal{\cal G}^{\otimes}$ the product of ${\mathcal{G}}_{H,\beta_{i}}$s
with uniform weights i.e., $\omega_{i}=\frac{1}{\ell},\forall i$ (see 2.2).
SuperChainTraceGibbs calls  RelMeanEst $({\cal{\cal G}^{\otimes}},F)$ and
RelMeanEst $({\cal{\cal G}^{\otimes}},G)$, with appropriate parameters, and
simply outputs the ratio of the two estimates (see alg. 4, left).

1:procedure

SuperChainTraceGibbs

(…)

2:
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})\leftarrow\scalebox{0.92}[0.98]{{Tpa}}(k,d)$†††
$k=\Theta(\log H_{\max})$ and $d=64$ as in [43]

3: $\varepsilon^{\prime}\leftarrow\frac{\varepsilon}{2+\varepsilon}$;
$\delta^{\prime}\leftarrow\frac{\delta}{2}$

4: for $i\in 1,2,\dots,\ell$ do

5: $f_{i}(x)\doteq\exp(-\frac{\beta_{i+1}-\beta_{i}}{2}H(x))$

6: $g_{i}(x)\doteq\exp(\frac{\beta_{i}-\beta_{i-1}}{2}H(x))$

7: end for

8: $F\doteq\bigotimes_{i=1}^{\ell}f_{i}$;
$G\doteq\bigotimes_{i=1}^{\ell}g_{i}$

9:
${\mathcal{G}}^{\otimes}\leftarrow\bigotimes_{i=1}^{\ell}{\mathcal{G}}_{H,\beta_{i}}$,
with $\omega_{i}=\frac{1}{\ell},\forall i$

10:
$R_{f}\leftarrow\exp(-\frac{\beta-\beta_{0}}{2}H_{\min})-\exp(-\frac{\beta-\beta_{0}}{2}H_{\max})$

11:
$R_{g}\leftarrow\exp(\frac{\beta-\beta_{0}}{2}H_{\max})-\exp(\frac{\beta-\beta_{0}}{2}H_{\min})$

12:
$\hat{\mu}\leftarrow\scalebox{0.92}[0.98]{{RelMeanEst}}({\mathcal{G}}^{\otimes},R_{f},T,F,\varepsilon^{\prime},\delta^{\prime})$

13:
$\hat{\nu}\leftarrow\scalebox{0.92}[0.98]{{RelMeanEst}}({\mathcal{G}}^{\otimes},R_{g},T,G,\varepsilon^{\prime},\delta^{\prime})$

14: return $\hat{Z}\leftarrow\frac{\hat{\nu}}{\hat{\mu}}$

15:end procedure

16:procedure

ParallelTraceGibbs

(…)

17:
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})=\scalebox{0.92}[0.98]{{Tpa}}(k,d)$

18:
$\varepsilon^{\prime}\leftarrow\frac{\sqrt[\ell]{1+\varepsilon}-1}{\sqrt[\ell]{1+\varepsilon}+1}$;
$\delta^{\prime}\leftarrow\frac{\delta}{2\ell}$

19: for $i\in 1,2,\dots\ell$ do

20: $f_{i}(x)\doteq\exp(-\frac{\beta_{i+1}-\beta_{i}}{2}H(x))$

21: $g_{i-1}(x)\doteq\exp(\frac{\beta_{i}-\beta_{i-1}}{2}H(x))$

22:
$R_{f}\leftarrow\exp(-\frac{\beta_{i+1}-\beta_{i}}{2}H_{\min})-\exp(-\frac{\beta_{i+1}-\beta_{i}}{2}H_{\max})$

23:
$R_{g}\leftarrow\exp(\frac{\beta_{i+1}-\beta_{i}}{2}H_{\max})-\exp(\frac{\beta_{i+1}-\beta_{i}}{2}H_{\min})$

24: $\hat{\mu}_{i}\leftarrow\scalebox{0.92}[0.98]{{RelMeanEst}}({\cal
G}_{i},R_{f},T_{i},f_{i},\varepsilon^{\prime},\delta^{\prime})$

25: $\hat{\nu}_{i}\leftarrow\scalebox{0.92}[0.98]{{RelMeanEst}}({\cal
G}_{i},R_{g},T_{i},g_{i},\varepsilon^{\prime},\delta^{\prime})$

26: end for

27: return
$\hat{Z}\leftarrow{\prod_{i=1}^{\ell}\frac{\hat{\nu}_{i}}{\hat{\mu}_{i}}}$

28:end procedure

Algorithm 2

SuperChainTraceGibbs

and

ParallelTraceGibbs

We now show the correctness and efficiency of  SuperChainTraceGibbs . Let
$\tau_{\rm prx}$ denote ${\mathcal{G}}^{\otimes}$’s true (and unknown)
relaxation time and $T$ a known upper-bound on it ($T\geq\tau_{\rm prx}$),
$\varepsilon$ and $\delta$ are user specified precision parameters. For
simlicity of presentation we use the following notation to refer to relative
ranges: ${\rm relR}={\rm Range}(F)/\mu+{\rm Range}(G)/\nu$, where
$\mu=\operatorname*{\mathbb{E}}[F]$ and $\nu=\operatorname*{\mathbb{E}}[G]$ .

######  Theorem 2.4.

With probability at least $1-\delta$, it holds that the total number $\hat{m}$
of Markov chain steps taken by  SuperChainTraceGibbs is upper-bounded by

$\mathcal{\tilde{O}}\Biggl{(}\ln\biggl{(}\frac{1}{\delta}\biggr{)}\Biggl{(}\frac{T\cdot{\rm
relR}}{\varepsilon}+\frac{{\tau_{\rm
prx}}\cdot\bigl{(}\mathrm{Reltrv}^{\tau_{\rm prx}}_{{\cal
G}^{\otimes}}(F)+\mathrm{Reltrv}^{\tau_{\rm prx}}_{{\cal
G}^{\otimes}}(G)\bigr{)}}{\varepsilon^{2}}\Biggr{)}\Biggr{)}~{}.$

######  Lemma 2.5.

Defining $\alpha_{1}=\sqrt{\frac{Z(\beta_{0})}{Z(\beta_{0}-\Delta_{\max})}}$,
we have: $\hskip 42.67912pt\frac{{\rm
Range}(F)}{\mu}\leq\alpha_{1}\sqrt{\frac{Q}{\exp(\Delta H_{\min})}}$ and
$\frac{{\rm Range}(G)}{\nu}\leq\alpha_{1}\sqrt{\frac{\exp(\Delta
H_{\max})}{Q}}$ .

Using lemma 2.5 and thm. 2.4, we identify $\varepsilon_{0}$ such that for
$\varepsilon\leq\varepsilon_{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 [43] we employ eq. 3 of lemma 2.1.
Finally we use thm. 1.1 and conclude:

######  Corollary 2.6.

Let $\alpha_{1}$ be as in lemma 2.5, $\tau_{\max}\doteq\max_{i}\tau_{i}$ and
$\varepsilon_{0}\doteq({\tau_{\rm prx}}/T)\cdot\left(\sqrt{\frac{\exp(\Delta
H_{\min})}{Q}}+\sqrt{\frac{Q}{\exp(\Delta H_{\max})}}\right)\cdot\alpha_{1}$.
When $\varepsilon\leq\varepsilon_{0}$, the number of Markov chain steps of
SuperChainTraceGibbs is dominated by $\tilde{O}(\ell\tau_{\max})$ .

#####  ParallelTraceGibbs

For $i=1,2,\dots,\ell-1$,  ParallelTraceGibbs (alg. 4, right) runs
$\scalebox{0.92}[0.98]{{RelMeanEst}}({\mathcal{G}}_{H,\beta_{i}},f_{\beta_{i},\beta_{i+1}})$
and
$\scalebox{0.92}[0.98]{{RelMeanEst}}({\mathcal{G}}_{H,\beta_{i}},g_{\beta_{i},\beta_{i+1}})$
independently. We show the computational complexity of  ParallelTraceGibbs in
thm. 2.7.

For $i=1,2,\dots,\ell$, assume $\tau_{i}$ is the true (unknown) relaxation
time of ${\mathcal{G}}_{H,\beta_{i}}$ and $T_{i}$ is a known bound on it. For
simplicity of presentation we use the following notations: ${\rm
relR}_{i}\doteq{\rm Range}(f_{\beta_{i},\beta_{i+1}})/\mu_{i}+{\rm
Range}(g_{\beta_{i-1},\beta_{i}})/\nu_{i}$, where
$\mu_{i}=\mathbb{E}(f_{\beta_{i},\beta_{i+1}})$ and
$\nu_{i}=\mathbb{E}(g_{\beta_{i},\beta_{i+1}})$.

######  Theorem 2.7 (Efficiency of  ParallelTraceGibbs ).

With probability at least $1-\delta$, it holds that the total number $\hat{m}$
of Markov chain steps taken by  ParallelTraceGibbs is upper-bounded by

$\displaystyle\mathcal{\tilde{O}}\Biggl{(}\log\left(\frac{\ell}{\delta}\right)\sum_{i=1}^{\ell}\left(\frac{\ell\cdot{T_{i}\cdot\rm
relR}_{i}}{\varepsilon}+\frac{\ell^{2}}{\varepsilon^{2}}\tau_{i}\cdot\left(\mathrm{Reltrv}^{\tau_{i}}_{{\mathcal{G}}_{H,\beta_{i}}}(f_{\beta_{i},\beta_{i+1}})+\mathrm{Reltrv}^{\tau_{i}}_{{\mathcal{G}}_{H,\beta_{i}}}(g_{\beta_{i-1},\beta_{i}})\right)\right)\Biggr{)}~{}.$

Furthermore, for all $1\leq i\leq\ell$, ${\rm
Range}(f_{\beta_{i},\beta_{i+1}})/\mu_{i}\leq\ell^{1/\log(n)}$ and ${\rm
Range}(g_{\beta_{i-1},\beta_{i}})/\nu_{i}\leq\ell^{\alpha_{0}(i)/\log n}$,
where $\alpha_{0}(i)=(\nicefrac{{H_{\max}}}{{2\mathbb{E}[H(x)]}})-1,$ for
$x\sim{\pi_{\beta_{i}}}~{}$.

ParallelTraceGibbs

and  SuperChainTraceGibbs make different computational complexity tradeoffs.
ParallelTraceGibbs is usually slower than  SuperChainTraceGibbs , because in
each iteration $i=1,2,\dots,\ell$, the mean estimator must acquire a higher-
precision estimate so that _all estimators together_ achieve an
$\varepsilon$-$\delta$ 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 ($\sum_{i=1}^{\ell}\tau_{i}$ vs.
$\max_{i=1,\dots,\ell}\tau_{i}$). Furthermore, dominance of the trace variance
terms in both of these algorithms occur at different values of $\varepsilon$.
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.

ParallelTraceGibbs |  SuperChainTraceGibbs |  Tpa \+ PPE [43]
---|---|---
$\displaystyle\ell^{2}\sum_{i=1}^{\ell}\tau_{i}\left({\rm Reltrv}^{\tau_{i}}_{{\mathcal{G}}_{H,\beta_{i}}}(f_{i})+\hbox{}\right.$  |  $\displaystyle\tau_{\rm prx}\left({\rm Reltrv}^{\tau_{\rm prx}}[F]+{\rm Reltrv}^{\tau_{\rm prx}}[G]\right)$ |  $\displaystyle\ln\frac{q\ln H_{\max}}{\varepsilon}\sum_{i=1}^{\ell}{T_{i}}\cdot\left({\mathbb{V}}_{\rm rel}(F)+{\mathbb{V}}_{\rm rel}(G)\right)$
$\left.{\rm Reltrv}^{\tau_{i}}_{{\mathcal{G}}_{H,\beta_{i}}}(g_{i})\right)$ | $\displaystyle=O\left(\ell\max\\{\tau_{i}\\}_{i=1\mathrel{\mathop{\mathchar 58\relax}}\ell}\right)$ | $\displaystyle=O\left(\ln\frac{q\ln H_{\max}}{\varepsilon}\sum_{i=1}^{\ell}{T_{i}}\right)$
Table 1: Comparison of the number of Markov chain steps, when $\varepsilon$ is
adequately small. In all columns, a multiplicative factor of
$\nicefrac{{1}}{{\varepsilon^{2}}}$ is omitted to ease presentation, and
$q=\ln Q$. Note that computational complexity of both  ParallelTraceGibbs and
SuperChainTraceGibbs only depends on true relaxation times, denoted by
$\tau_{i}$, and the  Tpa \+ PPE method’s complexity is dependent on their
upper bounds, denoted by $T_{i}$.

##  3 Experimental Results

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 [43].

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\times n$, the Hamiltonian is defined on $n^{2}$ random variables having
values $\pm 1$ and their dependency is represented by the Hamiltonian:
$H(x)=-\sum_{(i,j)\in E}{\mathbbm{1}(x(i)=x(j))}$. We run the algorithms on
lattices of sizes $2\times 2$, $3\times 3$, $4\times 4$, and $6\times 6$. For
each lattice, the parameter $\beta\geq 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\in\\{-1,1\\}$, and the voter variables
$T_{1},T_{2},\dots,T_{n}$ and $F_{1},F_{2},\dots,F_{n}$ all in $\\{0,1\\}$.
The factors have $2n+1$ weights,
$\omega,\omega_{T_{i}},\omega_{F_{i}},i=1,\dots,n$. The Hamiltonian is:

$\displaystyle H(Q,T,F)=\omega Q\max_{i}{T_{i}}-\omega
Q\max_{i}{F_{i}}+\sum_{i=1}^{n}\omega_{T_{i}}T_{i}+\sum_{i=1}^{n}\omega_{F_{i}}F_{i}\textrm{
, where }\omega,\omega_{T_{i}},\omega_{F_{i}}\in[-1,1]$

The parameters are reported in fig. 2. We follow De Sa et al. [19] 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 [43]. At each iteration of  RelMeanEst , the sample size is extended
with geometric ratio $1.1$ (see alg. 3 line 8). All code is available at
https://github.com/zysophia/Doubly_Adaptive_MCMC.

(a) $\beta=.01$, 3$\times$3 lattice

(b) $\beta=.002$, 6$\times$6 lattice

(c) relative errors

Figure 1: Comparison of sample complexity and precision
$\frac{1}{\varepsilon}$ on Ising models. See also the A.6. of the
supplementary material

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. 2(c) and 2(a) for the voting model, when
$\varepsilon$ is fixed and $Z{({\beta})}$ is varying (fig. 2(c)), and when $Z$
is fixed and $\varepsilon$ is varying (fig. 2(a)). All of these hold while the
precision of our algorithms beats [43] as $\varepsilon\rightarrow 0$ (fig.
1(c)).

(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\leftarrow 1$ in line 4 of  RelMeanEst
. In Figure 2(b), we show the effectiveness of _trace averaging_ , since both
SuperChainTraceGibbs and  ParallelTraceGibbs beat their simplified versions
($T\leftarrow 1$) after $\nicefrac{{1}}{{\varepsilon}}$ 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. 2(b), we show the trace
variance term  ParallelTraceGibbs becomes dominant earlier as
$\nicefrac{{1}}{{\varepsilon}}$ 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 .

(a) Complexity vs. $\varepsilon$ comparison against  Tpa \+ PPE [43].

(b) Comparison of our algorithms and the effect of trace variance.

(c) Complexity vs. $Z(\beta)$; $\varepsilon=0.025$ comparison against  Tpa \+
PPE [43].

Figure 2: Experiments on voting models. In (a) and (b) the parameters are
$\beta=0.1,n=3,\omega=0.9,\omega_{T}=\langle 0.2,0.5,0.1\rangle$ and
$\omega_{F}=-\langle 0.8,0.2,0.9\rangle$. In (c), we have $n=5$, and the
weights and $\beta$ are picked randomly to generate models with various values
of $Z(\beta)$.

##  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 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
$\Theta(\nicefrac{{\ell}}{{\varepsilon^{2}}})$ approximately independent
samples, where $\ell$ is the length of cooling schedule. This requires tight
convergence (total variance distance of
$O(\nicefrac{{\varepsilon^{2}}}{{\ell}})$ from stationary) for each sample,
which adds a multiplicative $\ln\frac{\ell}{\varepsilon^{2}}$, with
$\ell=\Theta(\ln Q\ln H_{\max})$, to its complexity (see column 3 of table 1
and [43], theorem 9). In contrast, our doubly adaptive method _only_ depends
on relaxation times, which do not depend on $\varepsilon$.

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
$\ell\max\\{\tau_{i}\\}_{i=1}^{\ell}$, as opposed to
$\sum_{i=1}^{\ell}\tau_{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.

1:procedure

RelMeanEst

2: Input: Markov chain $\mathcal{M}$, upper-bound on relaxation time $T$,
real-valued function $f$ with range $[a,b]$, letting $R=b-a$, multiplicative
precision $\varepsilon$, error probability $\delta$.

3: Output: Multiplicative approximation $\hat{\mu}$ of
$\mu=\operatorname*{\mathbb{E}}_{\pi}[f]$.

4: $T\leftarrow\left\lceil\frac{1+\Lambda}{1-\Lambda}\ln\sqrt{2}\right\rceil$;
$\Lambda^{\prime}\leftarrow\Lambda^{T}$ $\triangleright$ Choose $T$ to be an
upperbound on relaxation time

5: $I\leftarrow
1\vee\left\lfloor\log_{2}\left(\frac{bR}{2a^{2}}\cdot\frac{(1-\varepsilon)^{2}}{(1+\varepsilon)\varepsilon}\right)\right\rfloor;\
\alpha\leftarrow\frac{(1+\Lambda^{\prime})R\ln\frac{3I}{\delta}(1+\varepsilon)}{(1-\Lambda^{\prime})b\varepsilon};\
m_{0}\leftarrow 0$ $\triangleright$ Initialize _sampling schedule_

6: $T_{\rm unif}\leftarrow\left\lceil T\cdot\ln(\nicefrac{{1}}{{\pi_{\rm
min}}})\right\rceil$;
$(\vec{X}_{0,1},\vec{X}_{0,2})\leftarrow\mathcal{M}^{T_{\rm unif}}(\bot)$
$\triangleright$ Warm-start two chains for $T_{\rm unif}$ steps from arbitrary
$\bot\in\Omega$

7: for $i\in 1,2,\dots,I$ do

8: $m_{i}\leftarrow\left\lceil\alpha r^{i}\right\rceil$ $\triangleright$ Total
sample count at iteration $i$; $r$ is the geometric ratio (constant, usually
2) size

9: for $j\in(m_{i-1}+1),\dots,m_{i}$ do

10: $(\vec{X}_{j,1},\vec{X}_{j,2})\leftarrow(T$ steps of $\mathcal{M}$
starting at $\vec{X}_{j-1,1},\vec{X}_{j-1,2})$ $\triangleright$ Run two
independent copies of $\mathcal{M}$ for $T$ steps

11:
$\displaystyle\bar{f}(\vec{X}_{j,1})\leftarrow\frac{1}{T}\sum_{t=1}^{T}f\bigl{(}\vec{X}_{j,1}(t)\bigr{)}$;
$\displaystyle\bar{f}(\vec{X}_{j,2})\leftarrow\frac{1}{T}\sum_{t=1}^{T}f\bigl{(}\vec{X}_{j,2}(t)\bigr{)}$
$\triangleright$ Average $f$ over $T$-traces

12: end for

13:
$\displaystyle{\hat{\bm{\mu}}}_{i}\leftarrow\frac{1}{2m_{i}}\sum_{i=1}^{m_{i}}\bigl{(}f(\vec{X}_{j,1})+f(\vec{X}_{j,2})\bigr{)}$;
$\displaystyle\hat{v}_{i}\leftarrow\frac{{1}}{2{m_{i}}}\sum_{i=1}^{m_{i}}\bigl{(}(f(\vec{X}_{j,1})-f(\vec{X}_{j,2})\bigr{)}^{2}$
$\triangleright$ Compute empirical mean; trace variance

14: $\displaystyle
u_{i}\leftarrow\hat{v}_{i}+\frac{(11+\sqrt{21})(1+\nicefrac{{\Lambda^{\prime}}}{{\sqrt{21}}})R^{2}\ln\frac{3I}{\delta}}{(1-\Lambda^{\prime})m_{i}}+\sqrt{\frac{(1+\Lambda^{\prime})R^{2}\hat{v}_{i}\ln\frac{3I}{\delta}}{(1-\Lambda^{\prime})m_{i}}}$$\triangleright$
Variance upper bound

15:
$\displaystyle\hat{\bm{\varepsilon}}^{+}_{i}\leftarrow\frac{10R\ln\frac{3I}{\delta}}{(1-\Lambda^{\prime})m_{i}}+\sqrt{\frac{(1+\Lambda^{\prime})u_{i}\ln\frac{3I}{\delta}}{(1-\Lambda^{\prime})m_{i}}}$
$\triangleright$ Apply Bernstein bound

16:
$\displaystyle\hat{\bm{\mu}}^{\times}_{i}\leftarrow\frac{(\hat{\bm{\mu}}_{i}-\hat{\bm{\varepsilon}}^{+}_{i})\vee
a+(\hat{\bm{\mu}}_{i}+\hat{\bm{\varepsilon}}^{+}_{i})\wedge b}{2}$
$\triangleright$ Optimal mean estimate

17:
$\displaystyle\hat{\bm{\varepsilon}}^{\times}_{i}\leftarrow\frac{((\hat{\bm{\mu}}_{i}+\hat{\bm{\varepsilon}}^{+}_{i})\wedge
b-(\hat{\bm{\mu}}_{i}-\hat{\bm{\varepsilon}}^{+}_{i})\vee
a}{2\hat{\bm{\mu}}^{\times}_{i}}$ $\triangleright$ Empirical relative error
bound

18: if $(i=I)\vee(\hat{\bm{\varepsilon}}^{\times}_{i}\leq\varepsilon)$ then
$\triangleright$ Terminate if accuracy guarantee is met

19: return $\hat{\bm{\mu}}_{i}^{\times}$

20: end if

21: end for

22:end procedure

Algorithm 3

RelMeanEst

1:procedure

SuperChainTraceGibbs

(…)

2:
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})\leftarrow\scalebox{0.92}[0.98]{{Tpa}}(k,d)$‡‡‡
$k=\Theta(\log H_{\max})$ and $d=64$ as in [43]

3: $\varepsilon^{\prime}\leftarrow\frac{\varepsilon}{2+\varepsilon}$;
$\delta^{\prime}\leftarrow\frac{\delta}{2}$

4: for $i\in 1,2,\dots,\ell$ do

5: $f_{i}(x)\doteq\exp(-\frac{\beta_{i+1}-\beta_{i}}{2}H(x))$

6: $g_{i}(x)\doteq\exp(\frac{\beta_{i}-\beta_{i-1}}{2}H(x))$

7: end for

8: $F\doteq\bigotimes_{i=1}^{\ell}f_{i}$;
$G\doteq\bigotimes_{i=1}^{\ell}g_{i}$

9:
${\mathcal{G}}^{\otimes}\leftarrow\bigotimes_{i=1}^{\ell}{\mathcal{G}}_{H,\beta_{i}}$,
with $\omega_{i}=\frac{1}{\ell},\forall i$

10:
$R_{f}\leftarrow\exp(-\frac{\beta-\beta_{0}}{2}H_{\min})-\exp(-\frac{\beta-\beta_{0}}{2}H_{\max})$

11:
$R_{g}\leftarrow\exp(\frac{\beta-\beta_{0}}{2}H_{\max})-\exp(\frac{\beta-\beta_{0}}{2}H_{\min})$

12:
$\hat{\mu}\leftarrow\scalebox{0.92}[0.98]{{RelMeanEst}}({\mathcal{G}}^{\otimes},R_{f},T,F,\varepsilon^{\prime},\delta^{\prime})$

13:
$\hat{\nu}\leftarrow\scalebox{0.92}[0.98]{{RelMeanEst}}({\mathcal{G}}^{\otimes},R_{g},T,G,\varepsilon^{\prime},\delta^{\prime})$

14: return $\hat{Z}\leftarrow\frac{\hat{\nu}}{\hat{\mu}}$

15:end procedure

16:procedure

ParallelTraceGibbs

(…)

17:
$(\beta_{0},\beta_{1},\dots,\beta_{\ell})=\scalebox{0.92}[0.98]{{Tpa}}(k,d)$

18:
$\varepsilon^{\prime}\leftarrow\frac{\sqrt[\ell]{1+\varepsilon}-1}{\sqrt[\ell]{1+\varepsilon}+1}$;
$\delta^{\prime}\leftarrow\frac{\delta}{2\ell}$

19: for $i\in 1,2,\dots\ell$ do

20: $f_{i}(x)\doteq\exp(-\frac{\beta_{i+1}-\beta_{i}}{2}H(x))$

21: $g_{i-1}(x)\doteq\exp(\frac{\beta_{i}-\beta_{i-1}}{2}H(x))$

22:
$R_{f}\leftarrow\exp(-\frac{\beta_{i+1}-\beta_{i}}{2}H_{\min})-\exp(-\frac{\beta_{i+1}-\beta_{i}}{2}H_{\max})$

23:
$R_{g}\leftarrow\exp(\frac{\beta_{i+1}-\beta_{i}}{2}H_{\max})-\exp(\frac{\beta_{i+1}-\beta_{i}}{2}H_{\min})$

24: $\hat{\mu}_{i}\leftarrow\scalebox{0.92}[0.98]{{RelMeanEst}}({\cal
G}_{i},R_{f},T_{i},f_{i},\varepsilon^{\prime},\delta^{\prime})$

25: $\hat{\nu}_{i}\leftarrow\scalebox{0.92}[0.98]{{RelMeanEst}}({\cal
G}_{i},R_{g},T_{i},g_{i},\varepsilon^{\prime},\delta^{\prime})$

26: end for

27: return
$\hat{Z}\leftarrow{\prod_{i=1}^{\ell}\frac{\hat{\nu}_{i}}{\hat{\mu}_{i}}}$

28:end procedure

Algorithm 4

SuperChainTraceGibbs

and

ParallelTraceGibbs

##  Appendix A Appendix

###  A.1  Algorithms used in the literature

####  A.1.1 The  Tpa method [37, 43]

We refer to Huber and Schott’s algorithm as the original  Tpa , and
Kolmogorov’s, which is used in our algorithms and referred to as  Tpa $(k,d)$
in the main manuscript, as the  Tpa method.

1:output a schedule $(\beta_{1},\dots,\beta_{l})$ of values in the interval
$[\beta_{\min},\beta_{\max}]$.

2:$\beta_{0}\leftarrow\beta_{\min}$

3:for $i=0\mathrel{\mathop{\mathchar 58\relax}}\infty$ do

4: sample $X\sim\pi_{\beta_{i}}$ draw $U\in[0,1]$ uniformly,
$\beta_{i+1}=\beta_{i}-\log U/H(X)$ (or $+\infty$ if $H(X)=0$.)

5: if $\beta_{i+1}\notin[\beta_{\min},\beta_{\max}]$ then Terminate

6: end if

7:end for

Algorithm 5 The Original Tpa-method [37]

1:input integers $k$ and $d$

2:output a schedule $(\beta_{0},\beta_{1},\dots,\beta_{l})$ of values in the
interval $[\beta_{\min},\beta_{\max}]$.

3:for $i=1\mathrel{\mathop{\mathchar 58\relax}}k$ do

4: ${\cal B}_{i}\leftarrow\textsc{The Original TPA-method}()$.

5: let ${\cal B}\leftarrow{\cal B}\cup{\cal B}_{i}$

6:end for

7:sort $\cal B$, keep one sample uniformly from the initial $d$ elements, and
keep every $d$th successive value in the remaining sequence.

8:add $\beta_{\min}$ and $\beta_{\max}$ to $\cal B$ return ${\cal B}$

Algorithm 6 Tpa-Method [43]

####  A.1.2 Single site Gibbs sampler (Glauber dynamics chain)

Consider $\beta$ and $H$ defined as above. Let $X=(X_{1},X_{2},\dots,X_{n})$
be the set of all variables in the Gibbs distribution with inverse temperature
$\beta$ and Hamiltonian $H$, thus, the domain of $H$ is
$\Omega=\Omega_{1}\times\Omega_{2}\times\dots\Omega_{n}$, and each
$\Omega_{i}$ is the range of random variable $X_{i}$. At each time step $t$,
assume the current state is $x^{(t)}=(x_{1},x_{2},\dots,x_{n})$. Take $i\sim
1,\dots,n$ uniformly at random. Sample $y$ from the following distribution:

$\pi_{\beta}(y|x^{(t)}_{-i})=\frac{\exp(-\beta H(x^{(t)};x_{i}\leftarrow
y))}{\sum_{\omega\in\Omega_{i}}\exp(-\beta
H(x^{(t)};x_{i}\leftarrow\omega))}\enspace,$ (6)

where for an arbitrary $\omega\in\Omega_{i}$ we define
$(x^{(t)};x_{i}\leftarrow\omega)$ be the vector in which all the elements
except the $i$th element are equal to $x_{i}$ and the $i$th element is
replaced with $\omega$.

In other words, for any arbitrary vectors $x^{(t)}$ and $x^{(t+1)}$, the
transition probability is:

${\mathcal{G}}_{H,\beta}(x^{(t)},x^{(t+1)})=\begin{cases}(\nicefrac{{1}}{{n}})\pi_{\beta}(y|x^{(t)}_{-i}),&\exists
y,i\text{ such that }x_{i}\neq y\text{ and }x^{(t+1)}=(x^{(t)};x_{i}\leftarrow
y),\\\
\sum_{i=1}^{n}(\nicefrac{{1}}{{n}})\pi_{\beta}(x_{i}|x_{-i}^{(t)})&\text{if
}x^{(t)}=x^{(t+1)},\\\ 0&\text{otherwise~{}.}\par\end{cases}$

###  A.2 Missing proofs: TPA and relative trace variance properties

######  Lemma A.1.

Let $z(\beta)\doteq\log\left(Z{({\beta})}\right)$, $d$ and $k$ the parameters
of the  Tpa method, and $\beta_{i}$ and $\beta_{i+1}$ two consecutive points
generated by $\scalebox{0.92}[0.98]{{Tpa}}(k,d)$, we have:

  1. 1.

For any $\varepsilon\geq 0$, we have
$\mathbb{P}(z(\beta_{j})-z(\beta_{j+1})\leq\varepsilon)\geq(1-\exp(-\varepsilon
k/d))^{d}\simeq 1-d\exp(-\varepsilon k/d)$ ,

  2. 2.

For any $\varepsilon\geq 0$,
$\mathbb{P}\left(\Delta_{i}\geq\nicefrac{{\varepsilon}}{{\mathbb{E}[H(x)]}}\right)\leq
d\exp(-\varepsilon k/d),$ where the expectation of $H(x)$ is taken with
respect to distribution $x\sim\pi_{\beta_{i+1}}$.

###### Proof of lemma A.1.

Note that $\scalebox{0.92}[0.98]{{Tpa}}(k,d)$ of [43] consists of $k$ parallel
runs of the original  Tpa of [37] and outputting a sub-sequence of elements
which are $d$ apart.

Let $(b_{i})$ be the sequence generated by $k$ parallel copies of the original
Tpa , thus $\Delta_{j}=\beta_{j+1}-\beta_{j}=b_{j+d}-b_{j}$.

We first show item 1 by bounding
$\mathbb{P}\left(z(b_{j})-z(b_{j+d})\geq\varepsilon\right)$, and using

$\mathbb{P}(b_{j+d}-b_{j}<\varepsilon)\geq\prod_{i=1}^{d}\mathbb{P}(b_{j+i}-b_{j+i-1}<\varepsilon/d)~{}.$

With the definition of the PPP, and using [30] we have
$z{(b_{i})}-z{(b_{i+1})}$ follows the exponential distribution with mean
$1/k$, thus
$\mathbb{P}(z{(b_{i})}-z{(b_{i+1})}\geq\varepsilon/d)=\exp(-\varepsilon
k/d)~{}.$ Therefore,

$\mathbb{P}(z(b_{j+d})-z(b_{j})<\varepsilon)\geq\prod_{i=1}^{d}\mathbb{P}\left(z(b_{j+i})-z(b_{j+i-1}\right)<\varepsilon/d)=(1-\exp(-\varepsilon
k/d))^{d}~{}.$

To see item 2 of the Lemma let $z^{\prime}(\beta)$ be the derivative of
$z(\cdot)$ with respect to $\beta$, which is
$z^{\prime}(\beta)=\sum_{x\in\Omega}-H(x)\exp(-\beta H(x))/Z(\beta)$, thus
$z^{\prime}(\beta)\leq 0$. Using the Cauchy–Schwarz inequality we have
$z^{\prime\prime}(\beta)=(\sum_{x\in\Omega}H^{2}(x)\exp({-\beta
H(x))}\sum_{x\in\Omega}\exp({-\beta H(x))}-(\sum_{x\in\Omega}-H(x)\exp({-\beta
H(x)})^{2})/Z^{2}(\beta)\geq 0$. Therefore,

$z^{\prime}(\beta_{i})<\frac{z(\beta_{i+1})-z(\beta_{i})}{\beta_{i+1}-\beta_{i}}<z^{\prime}(\beta_{i+1}),$

Thus,
$\beta_{i+1}-\beta_{i}<\frac{z(\beta_{i})-z(\beta_{i+1})}{-z^{\prime}(\beta_{i})}~{}$.
Note that $-z^{\prime}(\beta_{i})=\mathbb{E}[H(x)],x\sim\pi_{\beta_{i}}$.
Therefore, we have:

$\displaystyle\mathbb{P}\left(\Delta_{i}\leq\frac{\epsilon}{\mathbb{E}[H]}\right)$
$\displaystyle\geq\mathbb{P}\left(\frac{z(\beta_{i})-z(\beta_{i+1})}{-z^{\prime}(\beta_{i+1})}\leq\frac{\epsilon}{\mathbb{E}[H]}\right)$
$\displaystyle=\mathbb{P}\left({z(\beta_{i})-z(\beta_{i+1})}\leq{\epsilon}\right)$
$\displaystyle\geq(1-\exp(-\epsilon k/d))^{d}$

Thus $\mathbb{P}\left(\Delta_{i}\geq\frac{\epsilon}{\mathbb{E}[H]}\right)\geq
1-(1-\exp(-\epsilon k/d))^{d}\approx d\exp(-\epsilon k/d))\enspace.$ ∎

###### Proof of Lemma 2.1..

Note that by Thm 3.1. of [57] we have,
$\mathbb{E}[(\bar{f}(\vec{X}_{1\mathrel{\mathop{\mathchar
58\relax}}\tau})-\mathbb{E}(f))^{2}]\leq\frac{2\tau_{\rm
rx}}{\tau}\mathbb{V}[f]$. Dividing both sides by
$\left(\mathbb{E}(f)\right)^{2}$ we get the second part of the premise. The
first part concludes from setting $\tau=\tau_{\rm rx}$.

∎

###  A.3  RelMeanEst

RelMeanEst

in summary To employ progressive sampling, we start by a small sample size and
calculate the _empirical estimation_ of the variance at each iteration. We
estimate an upper bound on the trace variance based on its empirical
estimation, and using that we check a termination condition.

Our variance estimator is what Cousins et al. introduced, and is based on
running two independent chains. Each sample is obtained by taking a trace of
length $T$ (given upper-bound on relaxation time) and taking the average over
all observed values on that trace. Thus, _half the square difference_ of the
averages on the two chains is an _unbiased estimate_ of the _trace variance_.

Before showing the result, we state two key theorems from the literature,
which describe how our tail bounds work.

######  Theorem A.2 (Hoeffding-Type Bounds for Mixing Processes, (see Thm. 2.1
of [22])).

For any $\delta\in(0,1)$, we have

$\operatorname*{\mathbb{P}}\\!\left(\lvert\hat{\mu}-\mu\rvert\geq\sqrt{\frac{2(1+\lambda)({\frac{R^{2}}{4}})\ln(\frac{2}{\delta})}{(1-\lambda)m}}\right)\leq\delta\enspace.$
(7)

This implies sample complexity

$m_{H}(\lambda,R,\varepsilon,\delta)=\frac{1+\lambda}{1-\lambda}\ln({\frac{2}{\delta}})\frac{R^{2}}{2\varepsilon^{2}}\in\Theta\Bigl{(}\tau_{\rm
rx}\ln({\frac{1}{\delta}})\frac{R^{2}}{\varepsilon^{2}}\Bigr{)}\enspace.$

######  Theorem A.3 (Bernstein-Type Bound for Mixing Process [39, Thm. 1.2]).

For any $\delta\in(0,1)$, we have

$\operatorname*{\mathbb{P}}\left(|\hat{\mu}-\mu|\geq\frac{10R\ln(\frac{2}{\delta})}{(1-\lambda)m}+\sqrt{\frac{2(1+\lambda)v_{\pi}\ln(\frac{2}{\delta})}{(1-\lambda)m}}\right)\leq\delta\enspace.$
(8)

This implies sample complexity

$m_{B}(\lambda,R,v,\varepsilon,\delta)=\frac{2}{1-\lambda}\ln({\frac{2}{\delta}})\Bigl{(}\frac{5R}{\varepsilon}+\frac{(1+\lambda)v_{\pi}}{\varepsilon^{2}}\Bigr{)}\in\Theta\Bigl{(}\tau_{\rm
rx}\ln({\frac{1}{\delta}})\Bigl{(}\frac{R}{\varepsilon}+\frac{v_{\pi}}{\varepsilon^{2}}\Bigr{)}\Bigr{)}\enspace.$

We now show the main result.

###### Proof of Theorem 2.2.

Suppose confidence interval $[a,b]$. The interval endpoints, multiplicative
error $\varepsilon_{\times}$, and additive error $\varepsilon_{+}$ are related
as
$2\varepsilon_{+}=a\frac{1+\varepsilon_{\times}}{1-\varepsilon_{\times}}-a=a\frac{2\varepsilon_{\times}}{1-\varepsilon_{\times}}$,
depicted graphically below.

$0$$a$$\mu$$b$ Worst Case:
$\varepsilon_{+}=\frac{a\varepsilon_{\times}}{1-\varepsilon_{\times}}$
---
Arbitrary Case: $\varepsilon_{+}=\mu\varepsilon_{\times}$
---
Best Case:
$\varepsilon_{+}=\frac{b\varepsilon_{\times}}{1+\varepsilon_{\times}}$
---

We derive a geometric progressive sampling schedule such that the algorithm
draws sample sizes, ranging between optimistic and pessimistic (over unknown
variance and mean) upper and lower bounds on the sufficient sample size.

Using the Markov chain Bennett inequality [39], the best-case complexity,
assuming maximal expectation, and minimal variance, is

$\displaystyle m^{\downarrow}$ $\displaystyle\geq
m_{B}(\Lambda,R,0,\varepsilon_{+},\frac{2\delta}{3I})$
$\displaystyle\geq\frac{(1+\Lambda)R\ln\frac{3I}{\delta}}{(1-\Lambda)\varepsilon_{+}}=\frac{(1+\Lambda)R\ln\frac{3I}{\delta}(1+\varepsilon_{\times})}{b(1-\Lambda)\varepsilon_{\times}}\enspace.$

The worst-case complexity, then assuming minimal expectation, and maximal
variance, is

$\displaystyle m^{\uparrow}$ $\displaystyle\geq
m_{H}(\Lambda,R,\varepsilon_{+},\frac{2\delta}{3I})$
$\displaystyle\geq\frac{(1+\Lambda)R^{2}\ln\frac{3I}{\delta}}{2(1-\Lambda)\varepsilon_{+}^{2}}=\frac{(1+\Lambda)R^{2}\ln\frac{3I}{\delta}(1-\varepsilon_{\times})^{2}}{2(1-\Lambda)a^{2}\varepsilon_{\times}^{2}}\enspace,$

via the Markov chain Hoeffding’s inequality [46].

Consequently, a doubling schedule requires
$I=\left\lfloor\log_{2}\left(\frac{m^{\uparrow}}{m^{\downarrow}}\right)\right\rfloor=\left\lfloor\log_{2}\left(\frac{bR}{2a^{2}}\cdot\frac{(1-\varepsilon_{\times})^{2}}{(1+\varepsilon_{\times})\varepsilon_{\times}}\right)\right\rfloor$
steps.

All tail bounds on variances and means are hold simultanously with probability
at least $1-\delta$ (by union bound), and the doubling schedule never
overshoots the sufficient sample size by more than a constant factor, which
yields the stated guarantees.

The proof consists of two parts, in both we make derive our new bounds by
writing an $\varepsilon_{\times}$-multiplicative approximation in terms of an
$\varepsilon_{+}$-additive approximation.

In the worst-case, we _underestimate_ the true mean $\mu$ by a factor
$(1-\varepsilon_{\times})$, and thus require a radius
$\varepsilon_{+}=\varepsilon_{\times}(1-\varepsilon_{\times})\mu$ additive
confidence interval.

We first show the _correctness guarantee_.

Observe that the sampling schedule is selected such that the final iteration
$I$ of the algorithm will draw a sufficiently large sample (size
$m^{\uparrow}$) such that the Hoeffding inequality will yield such a
confidence interval, even for worst-case (minimal) $\mu$. Now observe that
over the course of the algorithm, in each iteration, 3 tail bounds are
applied; one to upper-bound the variance, and then two to upper and lower
bound the mean in terms of the variance bound) as in [16]. By union bound, all
$3I$ tail-bounds hold simultaneously with probability at least $1-\delta$,
thus when the algorithm terminates, it produces a correct answer with at least
said probability.

We now show the _efficiency guarantee_. Suppose we get $\hat{\mu}$ from
RelMeanEst , by guarantee of correctness of the algorithm, we have a lower
bound on $\hat{\mu}$, $\hat{\mu}\geq\mu(1-\varepsilon_{\times})$ with
probability at least $1-\delta$.

Furthermore, we have $\varepsilon_{+}=\mu\varepsilon_{\times}$ and
$\mathrm{trv}^{(\tau\mathrm{rx})}=(\mathrm{Reltrv}^{\tau_{\rm
rx}}-1)\times\hat{\mu}^{2}\geq(\mathrm{Reltrv}^{\tau_{\rm
rx}}-1)\mu^{2}(1-\varepsilon_{\times})^{2}$. For this ${\varepsilon_{+}}$, we
have via the Bernstein inequality that

$m^{*}\in\mathcal{O}\left(\log\left(\frac{\log(R/(\mu\varepsilon_{\times}))}{\delta}\right)\left(\frac{R/\mu}{(1-\Lambda)\varepsilon_{\times}}+\frac{\tau_{\rm
rx}(\mathrm{Reltrv}^{\tau_{\rm
rx}}-1)}{\varepsilon_{\times}^{2}}\right)\right)$

would be a sufficient sample size if (1) the algorithm were to draw a sample
of this size, and (2) we were to use the _true trace variance_ instead of the
_estimated upper-bound on trace variance_.

Fortunately, correcting for (1) adds a constant factor to the sample
complexity, as the first sample size $\alpha$ is selected to be twice the
minimal sufficient sample size $m^{\downarrow}$ (i.e., the sample size such
that no smaller sample size would be sufficient), and at each iteration the
sample size selected is double the previous (line 8). In other words, this
geometric grid will never overshoot any sample size by more than a factor 2.

Resolving (2) is a bit more subtle, but we now show that there is no
asymptotic change in replacing the variance with the estimated variance upper
bound (w.h.p.). First, note that the Bernstein bound is _bidirectional_ , so
it can just as well be used to upper-bound empirical variance with true
variance as to upper-bound true variance with empirical variance. We bound
true variance in terms of empirical variance on line 14, and note that here we
have

$v\leq
u\in\hat{v}+\mathcal{O}\left(\frac{R^{2}\ln\frac{I}{\delta}}{m}+\sqrt{\frac{R^{2}\hat{v}\ln\frac{I}{\delta}}{m}}\right)\enspace.$

Fortunately, the latter terms are negligible, as in line 15, we bound

$\displaystyle\varepsilon_{+}$
$\displaystyle\in\mathcal{O}\left(\frac{R\ln\frac{I}{\delta}}{m}+\sqrt{\frac{u\ln\frac{I}{\delta}}{m}}\right)$
$\displaystyle=\mathcal{O}\left(\frac{R\ln\frac{I}{\delta}}{m}+\sqrt{\frac{\biggl{(}\hat{v}+\mathcal{O}\biggl{(}{\frac{R\ln\frac{I}{\delta}}{m}}+{\sqrt{\frac{\hat{v}\ln\frac{I}{\delta}}{m}}}\biggr{)}\biggr{)}\ln\frac{I}{\delta}}{m}}\right)$
$\displaystyle=\mathcal{O}\left(\frac{R\ln\frac{I}{\delta}}{m}+\sqrt{\frac{\biggl{(}v+\mathcal{O}\biggl{(}{\frac{R\ln\frac{I}{\delta}}{m}}+{\sqrt{\frac{v\ln\frac{I}{\delta}}{m}}}\biggr{)}+\mathcal{O}\biggl{(}{\frac{R\ln\frac{I}{\delta}}{m}}+{\sqrt{\frac{\hat{v}\ln\frac{I}{\delta}}{m}}}\biggr{)}\biggr{)}\ln\frac{I}{\delta}}{m}}\right)$
(w.h.p.)
$\displaystyle=\mathcal{O}\left(\frac{R\ln\frac{I}{\delta}}{m}+\sqrt{\frac{{v}\ln\frac{I}{\delta}}{m}}\right)\enspace.$
(w.h.p.)

Putting these together, we thus have that, w.h.p., sample consumption is
bounded as

$\hat{m}\in
2\mathcal{O}(m^{*})=\mathcal{O}\left(\log\left(\frac{\log(R/(\mu\varepsilon_{\times}))}{\delta}\right)\left(\frac{R/\mu}{(1-\Lambda)\varepsilon_{\times}}+\frac{\tau_{\rm
rx}(\mathrm{Reltrv}^{\tau_{\rm
rx}}-1)}{\varepsilon_{\times}^{2}}\right)\right)\enspace.$

To conclude, we need only relate $T(\mathrm{Reltrv}^{T}-1)$ and $\tau_{\rm
rx}(\mathrm{Reltrv}^{\tau_{\rm rx}}-1)$. Letting $T$ as in line
LABEL:alg:tbound, note that since $T\geq\tau_{\rm rx}$, it holds that
$T(\mathrm{Reltrv}^{T}-1)\geq\tau_{\rm rx}(\mathrm{Reltrv}^{\tau_{\rm
rx}}-1)$, by the trace variance inequalities, which yields the result. ∎

###  A.4 Missing proofs from analysis of  SuperChainTraceGibbs

###### Proof of thm 2.4.

Follows immediately from thm. 2.2 and plugging in the values for paired
product estimators and the product chain. ∎

###### Full Proof of Lemma 2.8.

Let $\bar{\beta}_{i,i+1}\doteq\frac{\beta_{i}+\beta_{i+1}}{2}$, we have
${\mu_{i}}=\frac{Z(\bar{\beta}_{i,i+1})}{Z(\beta_{i})}$ and
${\nu_{i}}=\frac{Z(\bar{\beta}_{i,i+1})}{Z(\beta_{i+1})}$. Thus we have
$\nu=\frac{\prod_{i=1}^{\ell-1}Z(\bar{\beta}_{i,i+1})}{\prod_{i=1}^{\ell-1}Z(\beta_{i+1})}>1$,
$\mu=\frac{\prod_{i=1}^{\ell-1}Z(\bar{\beta}_{i,i+1})}{\prod_{i=1}^{\ell-1}Z(\beta_{i})}<1$.

Note that $\nu=\mu\frac{Z(\beta_{0})}{Z(\beta_{\max})}$, thus we proceed by
bounding $\mu$.

$\displaystyle\log\prod_{i=1}^{\ell-1}Z(\bar{\beta}_{i,i+1})$
$\displaystyle=\sum_{i=1}^{\ell-1}z(\bar{\beta}_{i,i+1})$ Taking $\log$
$\displaystyle\geq\sum_{i=1}^{\ell-1}z(\beta_{i})-\frac{\Delta_{i}}{2}\operatorname*{\mathbb{E}}_{x\sim\pi_{\beta_{i}}}[H(x)]$
Taylor expansion & that $\frac{\partial^{2}}{\partial\beta^{2}}z(\beta)>0$

Thus, by taking exponents we get:

$\displaystyle\prod_{i=1}^{\ell-1}Z(\bar{\beta}_{i,i+1})$
$\displaystyle\geq\exp\left(\sum_{i=1}^{\ell-1}z(\beta_{i})-\frac{\Delta_{i}}{2}\operatorname*{\mathbb{E}}_{x\sim\pi_{\bar{\beta}_{i,i+1}}}[H(x)]\right)$
$\displaystyle\geq\left(\prod_{i=1}^{\ell-1}Z(\beta_{i})\right)\exp\left(-\sum_{i=1}^{\ell-1}\frac{\Delta_{i}}{2}\operatorname*{\mathbb{E}}_{x\sim\pi_{\beta_{i}}}[H(x)]\right)$

Therefore,
$\mu=\frac{\prod_{i=1}^{\ell-1}Z(\bar{\beta}_{i,i+1})}{\prod_{i=1}^{\ell-1}Z(\beta_{i})}\geq\exp\left(-\sum_{i=1}^{\ell-1}\frac{\Delta_{i}}{2}\operatorname*{\mathbb{E}}_{x\sim\pi_{\beta_{i}}}[H(x)]\right)$.
Using this form, we now employ the fundamental theorem of calculus to prove
the premise:

Let $\Delta_{\max}\dot{=}\max_{i}\Delta_{i}$.

$\displaystyle\mu$
$\displaystyle\geq\exp\left(-\sum_{i=1}^{\ell-1}\frac{\Delta_{i}}{2}\operatorname*{\mathbb{E}}_{x\sim\pi_{\beta_{i}}}[H(x)]\right)$
$\displaystyle=\exp\left(-\sum_{i=1}^{\ell-1}\frac{\Delta_{i}}{2}\operatorname*{\mathbb{E}}_{x\sim\pi_{\beta_{i}}}[H(x)]\right)$
$\displaystyle\geq\exp\left(\frac{1}{2}\int_{\beta_{\min}-\Delta_{\max}}^{\beta_{\max}-\Delta_{\max}}-\operatorname*{\mathbb{E}}_{x\sim\pi_{\beta}}[H(x)]\,\mathrm{d}\beta\right)$
Increasing Integrand
$\displaystyle=\exp\left(\frac{1}{2}\bigl{(}z(\beta_{\max}-\Delta_{\max})-z(\beta_{\min}-\Delta_{\max})\bigr{)}\right)$
FTOC and that $z^{\prime}(\beta)=\mathbb{E}_{x\sim\pi_{\beta}H}$
$\displaystyle\geq\exp\left(\frac{1}{2}\bigl{(}z(\beta_{\max})-z(\beta_{\min}-\Delta_{\max})\bigr{)}\right)$
$z$ is Decreasing
$\displaystyle=\exp\left(\frac{1}{2}\bigl{(}z(\beta_{\max})-z(\beta_{\min})+z(\beta_{\min})-z(\beta_{\min}{-}\Delta_{\max})\bigr{)}\right)$
$\displaystyle\geq
Q^{-\frac{1}{2}}\sqrt{\frac{Z(\beta_{\min})}{Z(\beta_{\min}-\Delta_{\max})}}~{}.$

From the above we also conclude that $\nu\geq
Q^{1/2}\sqrt{\frac{Z(\beta_{\min})}{Z(\beta_{\min}-\Delta_{\max})}}$. Note
that ${\rm
Range}(f)=\exp(-\frac{\Delta}{2}H_{\min})-\exp(-\frac{\Delta}{2}H_{\max})\leq\sqrt{\exp(-\Delta
H_{\min})}$ and ${\rm
Range}(g)=\exp(\frac{\Delta}{2}H_{\max})-\exp(\frac{\Delta}{2}H_{\min})\leq\sqrt{\exp(\Delta
H_{\max})}$ . Thus the lemma is concluded.

∎

###### Proof of Corollary 2.6.

The corollary follows from thm 2.2 plugging in $R$ from lemma 2.5 and setting
$\tau_{\rm prx}=\ell\max_{i=1}^{\ell}\tau_{i}$ (see, e.g., [47]). ∎

###  A.5 Analysis of  ParallelTraceGibbs

Let $(\beta_{0},\beta_{1},\dots\beta_{l})$ be a cooling schedule generated by
Tpa $(k,d)$, where $k$ and $d$ are chosen as in [43]. For each $i$ let
$f_{\beta_{i},\beta_{i+1}}$, $g_{\beta_{i-1},\beta_{i}}$ be the paired
estimators corresponding to this schedule, and
$\mu_{i}=\mathbb{E}[f_{\beta_{i},\beta_{i+1}}]$,
$\nu_{i}=\mathbb{E}[g_{\beta_{i-1},\beta_{i}}]$ .  ParallelTraceGibbs
estimates $Q$ by running  RelMeanEst on each ${\mathcal{G}}_{H,\beta_{i}}$, to
estimate $\mu_{i}$ and $\nu_{i}$s each with precision
$\varepsilon^{\prime}=\nicefrac{{(\sqrt[l]{1+\varepsilon}-1)}}{{(\sqrt[l]{1+\varepsilon}+1)}}$.
Note that by this setting, $Q$ will be approximated within multiplicative
factor of
$\left(\nicefrac{{1+\varepsilon^{\prime}}}{{1-\varepsilon^{\prime}}}\right)^{\ell}$.
Assume $\tau_{i}$ is the true relaxation time of ${\mathcal{G}}_{H,\beta_{i}}$
and suppose $\Lambda_{i}$ is a known upper bound on the second eigenvalue of
${\mathcal{G}}_{H,\beta_{i}}$, thus $(\Lambda_{i}-1)^{-1}\log(2)\geq\tau_{i}$.
The following hold and thm 2.7 is immediately concluded from it:

######  Lemma A.4.

Let $H_{\max}\doteq\max_{x\in\Omega}H(x)$. we have:

  1. 1.

for all $1\leq i\leq\ell$, ${\rm
Range}(f_{\beta_{i},\beta_{i+1}})/\mu_{i}\leq\ell^{1/\log(n)}$,

  2. 2.

for all $1\leq i\leq\ell$, ${\rm
Range}(g_{\beta_{i-1},\beta_{i}})/\nu_{i}\leq\ell^{\alpha_{0}(i)/\log n}$,
where
$\alpha_{0}(i)=(\nicefrac{{H_{\max}}}{{2\mathbb{E}[H(x)]}})-1,~{}x\sim{\pi_{\beta_{i}}}$.

###### Proof.

Let $\Delta_{i}=\beta_{i+1}-\beta_{i}$. Thus,
$f_{i}(x)=\exp\left(\frac{-\Delta_{i}}{2}H(x)\right)$ and
$g_{i}(x)=\exp\left(\frac{\Delta_{i}}{2}H(x)\right)$. So we have:

${\rm
Range}(f_{i})=\exp\left(\frac{-\Delta_{i}}{2}\min_{x}H(x)\right)-\exp\left(\frac{-\Delta_{i}}{2}\max_{x}H(x)\right)\leq\exp\left(\frac{-\Delta_{i}}{2}\min_{x}H(x)\right)$

and

${\rm
Range}(g_{i})=\exp\left(\frac{\Delta_{i}}{2}\max_{x}H(x)\right)-\exp\left(\frac{\Delta_{i}}{2}\min_{x}H(x)\right)\leq\exp\left(\frac{\Delta_{i}}{2}\max_{x}H(x)\right)$

$\mu_{i}=Z(\beta_{i}+\Delta_{i}/2)/Z(\beta_{i})\quad\&\quad\nu_{i}=Z(\beta_{i+1}-\Delta_{i}/2)/Z(\beta_{i+1})$

$\displaystyle\frac{{\rm Range}(f_{i})}{\mu_{i}}$
$\displaystyle\leq\frac{\exp\left(\frac{-\Delta_{i}}{2}\min_{x}H(x)\right)}{\exp\left(z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})\right)}~{},\frac{{\rm
Range}(g_{i})}{\nu_{i}}$
$\displaystyle\leq\frac{\exp\left(\frac{\Delta_{i}}{2}\max_{x}H(x)\right)}{\exp\left(z(\beta_{i+1}-\Delta_{i}/2)-z(\beta_{i+1})\right)}$
(9)

Writing
$\Delta_{i}/2=\frac{\Delta_{i}/2}{z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})}(z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i}))$,
we get:

$\displaystyle\frac{{\rm Range}(f_{i})}{\mu_{i}}$
$\displaystyle\leq\exp\left(-\frac{\Delta_{i}}{2}\min_{x}H(x)-\left(z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})\right)\right)$
$\displaystyle\leq\exp\left(\left(z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})\right)\left(\frac{-\Delta_{i}\cdot\min_{x}H(x)}{2\left(z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})\right)}-1\right)\right)$

and

$\frac{{\rm
Range}(g_{i})}{\nu_{i}}\leq\exp\left(\left(z(\beta_{i+1}-\Delta_{i}/2)-z(\beta_{i+1})\right)\left(\frac{\Delta_{i}\cdot\max_{x}H(x)}{2\left(z(\beta_{i+1}-\Delta_{i}/2)-z(\beta_{i+1})\right)}-1\right)\right)$

Let $z^{\prime}$ and $z^{\prime\prime}$ be the first and second derivative of
$z$ with respect to $\beta$. Note that
$z^{\prime}(\beta)=\mathbb{E}_{x\sim\pi_{\beta}}[-H(x)]$. Since
$z^{\prime\prime}\geq 0$ we have:

$z^{\prime}(\beta_{i})<\frac{z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})}{\Delta_{i}/2}<z^{\prime}(\beta_{i}+\Delta_{i}/2)$

and

$z^{\prime}(\beta_{i+1}-\Delta_{i}/2)<\frac{z(\beta_{i+1})-z(\beta_{i+1}-\Delta_{i}/2)}{\Delta_{i}/2}<z^{\prime}(\beta_{i+1}).$

Which are equivalent to
$\frac{1}{z^{\prime}(\beta_{i}+\Delta_{i}/2)}\leq\frac{\Delta_{i}/2}{z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})}\leq\frac{1}{z^{\prime}(\beta_{i})}$
and
$\frac{1}{z^{\prime}(\beta_{i+1})}\leq\frac{\Delta_{i}/2}{z(\beta_{i+1})-z(\beta_{i+1}-\Delta_{i}/2)}\leq\frac{1}{z^{\prime}(\beta_{i+1}-\Delta_{i}/2)}$.

Therefore,

$\displaystyle\frac{{\rm Range}(f_{i})}{\mu_{i}}$
$\displaystyle\leq\exp\left(\left(z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})\right)\left(\frac{-\min_{x}H(x)}{2}\frac{1}{z^{\prime}(\beta_{i})}-1\right)\right)$
(10)
$\displaystyle=\exp\left(\left(z(\beta_{i}+\Delta_{i}/2)-z(\beta_{i})\right)\left(\frac{\min_{x}H(x)}{2}\frac{1}{\mathbb{E}[H]}-1\right)\right)$
(11)
$\displaystyle\leq\exp\left(z(\beta_{i})-z(\beta_{i}+\Delta_{i}/2)\right)$
(12)

Similarly for range of $g_{i}$s we have:

$\displaystyle\frac{{\rm Range}(g_{i})}{\nu_{i}}$
$\displaystyle\leq\exp\left(\left(z(\beta_{i+1}-\Delta_{i}/2)-z(\beta_{i+1})\right)\left(\frac{-\max_{x}H(x)}{2}\frac{1}{z^{\prime}(\beta_{i})}-1\right)\right)$
(13)
$\displaystyle\leq\exp\left(\left(z(\beta_{i+1}-\Delta_{i}/2)-z(\beta_{i+1})\right)\left(\frac{\max_{x}H(x)}{2\mathbb{E}_{\pi_{\beta_{i}}}[H(X)]}-1\right)\right)$
(14)

We now use (10) together with lemma A.1. Setting $d=1$ we have,

$\displaystyle\mathbb{P}\left(z{(\beta_{i})}-z{(\beta_{i+1})}>\frac{\log(3l/4)}{\log
n}\right)=\exp(-\frac{\log(3l/4)}{\log n}\cdot k)$
$\displaystyle=(3/4)\exp(-\log l/\log n(\log n))$ $\displaystyle=(3/4)(1/l).$

Using union bound over all $1\leq i\leq\ell$ and that
$z(\beta_{i})-z(\beta_{i+1})\geq z(\beta_{i})-z(\beta_{i}+\Delta_{i}/2)$, we
conclude that with probability at least $3/4$ we have that for all $f_{i}$,
${\rm Range}(f_{i})/\mu_{i}\leq\ell^{1/\log(n)}$.

Similarly using (13), the union bound, lemma A.1 and that
$z(\beta_{i})-z(\beta_{i+1})\geq z(\beta_{i}-\Delta_{i}/2)-z(\beta_{i+1})$, we
can show that with constant probability all $g_{i}$s generated by the  Tpa
schedule obey: $\forall g_{i};1\leq i\leq\ell,~{}{\rm
Range}(g_{i})/\nu_{i}\leq\exp\left((\log l/\log
n)\cdot(\alpha)\right)=\ell^{\alpha_{0}/\log n}$, where
$\alpha_{0}=\frac{\max_{x}H(x)}{2\mathbb{E}_{\pi_{\beta_{i}}}[H(X)]}-1$. ∎

The following corollary is concluded from lemma A.4 and relative trace
variance bounds:

######  Corollary A.5.

When
$\varepsilon\leq\ell^{1/\log(n)}(1+\ell^{\alpha_{0}(i)})\cdot\frac{\ell\tau_{\beta_{i}}}{(1-\Lambda_{i})^{-1}}$,
RelMeanEst invoked on the $i$th iteration will stop using sample consumption
of $\tilde{O}\left({\ell^{2}}\tau_{i}{\rm Reltrv}_{i}\right)$ note that this
is improvement over classic bounds which are
$\tilde{O}\left((1-\Lambda_{i})^{-1}\mathbb{V}{\rm rel}_{i}\right)$. In total
the sample complexity of  ParallelTraceGibbs for
$\varepsilon\leq\ell^{1/\log(n)}\min_{i}(1+\ell^{\alpha_{0}(i)})\cdot\frac{\ell\tau_{\beta_{i}}}{(1-\Lambda_{i})^{-1}}$
is dominated by $\tilde{O}\left({\ell^{2}}\sum_{i=1}^{\ell}\tau_{i}{\rm
Reltrv}_{i}\right)~{}.$

###  A.6 Further experimental results

(a) $\beta=.05$, 2$\times$2 lattice

(b) $\beta=.01$, 3$\times$3 lattice

(c) $\beta=.02$, 4$\times$4 lattice

(d) $\beta=.002$, 6$\times$6 lattice

Figure 3: Comparison of sample complexity on Ising models.

## References

  * [1] Permit allocation in emissions trading using the Boltzmann distribution.  Physica, A 391:4883–4890, 2012.
  * [2] H. Afshar, S. Sanner, and C. Webers.  Closed-form Gibbs sampling for graphical models with algebraic constraints.  In AAAI, 2016.
  * [3] D. Aldous, G. R. Grimmett, C. D. Howard, F. Martinelli, J. M. Steele, and L. Saloff-Coste.  Probability on discrete structures, volume 110.  Springer Science & Business Media, 2013.
  * [4] Y. Alimohammadi, N. Anari, K. Shiragur, and T. Vuong.  Fractionally log-concave and sector-stable polynomials: Counting planar matchings and more.  ArXiv, abs/2102.02708, 2021.
  * [5] N. Anari, K. Liu, and S. O. Gharan.  Spectral independence in high-dimensional expanders and applications to the hardcore model.  SIAM Journal on Computing, (0):FOCS20–1, 2021.
  * [6] I. Bezáková, D. Stefankovic, V. Vazirani, and E. Vigoda.  Accelerating simulated annealing for the permanent and combinatorial counting problems.  In SODA 2006, 2006.
  * [7] N. Bhatnagar, A. Bogdanov, and E. Mossel.  The computational complexity of estimating MCMC convergence time.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 424–435. Springer, 2011.
  * [8] N. Bhatnagar, A. Sly, and P. Tetali.  Reconstruction threshold for the hardcore model.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 434–447. Springer, 2010.
  * [9] A. Blanca, P. Caputo, Z. Chen, D. Parisi, D. Stefankovic, and E. Vigoda.  On mixing of Markov chains: Coupling, spectral independence, and entropy factorization.  ArXiv, abs/2103.07459, 2021.
  * [10] S. P. Brooks and A. Gelman.  General methods for monitoring convergence of iterative simulations.  Journal of computational and graphical statistics, 7(4):434–455, 1998.
  * [11] S. P. Brooks and G. O. Roberts.  Assessing convergence of Markov chain Monte Carlo algorithms.  Statistics and Computing, 8(4):319–335, 1998.
  * [12] Z. Chen, A. Galanis, D. Štefankovič, and E. Vigoda.  Rapid mixing for colorings via spectral independence.  In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1548–1557. SIAM, 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] K.-M. Chung, H. Lam, Z. Liu, and M. Mitzenmacher.  Chernoff-Hoeffding bounds for Markov chains: Generalized and simplified.  arXiv:1201.0559, 2012.
  * [15] B. A. Cipra.  An introduction to the Ising model.  The American Mathematical Monthly, 94(10):937–959, 1987.
  * [16] C. Cousins, S. Haddadan, and E. Upfal.  Making mean-estimation more efficient using an MCMC trace variance approach: DynaMITE.  CoRR, abs/2011.11129, 2020.
  * [17] C. De Sa, V. Chen, and W. 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.
  * [18] C. De Sa, K. Olukotun, and C. Ré.  Ensuring rapid mixing and low bias for asynchronous Gibbs sampling.  In Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48, ICML’16, pages 1567–1576, 2016.
  * [19] C. De Sa, C. Zhang, K. Olukotun, and C. Ré.  Rapidly mixing Gibbs sampling for a class of factor graphs using hierarchy width.  Advances in neural information processing systems, 28:3079–3087, 2015.
  * [20] A. Dixit and V. Roy.  MCMC diagnostics for higher dimensions using Kullback Leibler divergence.  Journal of Statistical Computation and Simulation, 87(13):2622–2638, 2017.
  * [21] H. Elliott, H. Derin, R. Cristi, and D. Geman.  Application of the Gibbs distribution to image segmentation.  9:678–681, 1984.
  * [22] J. Fan, B. Jiang, and Q. Sun.  Hoeffding’s lemma for Markov chains and its applications to statistical learning.  arXiv:1802.00211, 2018.
  * [23] G. S. Fishman.  Choosing sample path length and number of sample paths when starting in steady state.  Oper. Res. Lett., 16:209–219, 1994.
  * [24] S. Geman and D. Geman.  Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6(6):721–741, 1984.
  * [25] C. R. Gibbs.  Characterization and application of ferrozine iron reagent as a ferrous iron indicator.  Analytical Chemistry, 48(8):1197–1201, 1976.
  * [26] J. W. Gibbs.  Elementary Principles in Statistical Mechanics.  Scribner, 1902.
  * [27] J. Gonzalez, Y. Low, A. Gretton, and C. Guestrin.  Parallel Gibbs sampling: From colored fields to thin junction trees.  In G. Gordon, D. Dunson, and M. Dudík, 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.
  * [28] T. Griffiths.  Gibbs sampling in the generative model of latent Dirichlet allocation.  Technical report, 2002.
  * [29] S. Haddadan and P. Winkler.  Mixing of permutations by biased transpositions.  Theory of Computing Systems, 63(5):1068–1088, 2019.
  * [30] D. G. Harris and V. Kolmogorov.  Parameter estimation for Gibbs distributions.  CoRR, abs/2007.10824, 2020.
  * [31] B. He, C. De Sa, I. Mitliagkas, and C. Ré.  Scan order in Gibbs sampling: Models in which it matters and bounds on how much.  Advances in neural information processing systems, 29, 2016.
  * [32] A. Hellweg and F. Eckert.  Brick by brick computation of the Gibbs free energy of reaction in solution using quantum chemistry and COSMO-RS.  AIChE Journal, 63(9):3944–3954, 2017.
  * [33] T. Hrycej.  Gibbs sampling in Bayesian networks.  Artificial Intelligence, 46(3):351–363, 1990.
  * [34] D. Hsu, A. Kontorovich, D. A. Levin, Y. Peres, C. Szepesvári, and G. Wolfer.  Mixing time estimation in reversible Markov chains from a single sample path.  The Annals of Applied Probability, 29(4):2439–2480, 2019.
  * [35] M. Huber.  Approximation algorithms for the normalizing constant of Gibbs distributions.  The Annals of Applied Probability, 25(2):974–985, 2015.
  * [36] M. Huber and S. Schott.  Random construction of interpolating sets for high-dimensional integration.  Journal of Applied Probability, 51(1):92–105, 2014.
  * [37] M. Huber, S. Schott, et al.  Using TPA for Bayesian inference.  Bayesian Statistics, 9:257–282, 2010.
  * [38] M. Jerrum, L. G. Valiant, and V. V. Vazirani.  Random generation of combinatorial structures from a uniform distribution.  Theoretical Computer Science, 43:169–188, 1986.
  * [39] B. Jiang, Q. Sun, and J. Fan.  Bernstein’s inequality for general Markov chains.  arXiv:1805.10721, 2018.
  * [40] G. Karagiannis and C. Andrieu.  Annealed importance sampling reversible jump MCMC algorithms.  Journal of Computational and Graphical Statistics, 22(3):623–648, 2013.
  * [41] J. G. Kemeny, J. L. Snell, and A. W. Knapp.  Denumerable Markov chains: with a chapter of Markov random fields by David Griffeath, volume 40.  Springer Science & Business Media, 2012.
  * [42] D. Koller and N. Friedman.  Probabilistic Graphical Models: Principles and Techniques - Adaptive Computation and Machine Learning.  The MIT Press, 2009.
  * [43] V. Kolmogorov.  A faster approximation algorithm for the Gibbs partition function.  In Conference On Learning Theory, pages 228–249. PMLR, 2018.
  * [44] O. Krause, A. Fischer, and C. Igel.  Algorithms for estimating the partition function of restricted Boltzmann machines.  Artificial Intelligence, 278:103–195, 10 2019.
  * [45] P. S. La Rosa, T. L. Brooks, E. Deych, B. Shands, F. Prior, L. J. Larson-Prior, and W. D. Shannon.  Gibbs distribution for statistical analysis of graphical data with a sample application to fcmri brain images.  Statistics in medicine, 35(4):566—580, February 2016.
  * [46] C. Leon and F. Perron.  Optimal Hoeffding bounds for discrete reversible Markov chains.  The Annals of Applied Probability, 14, 05 2004.
  * [47] D. A. Levin and Y. Peres.  Markov chains and mixing times, volume 107.  American Mathematical Soc., 2017.
  * [48] P. Lezaud.  Chernoff-type bound for finite Markov chains.  The Annals of Applied Probability, 8, 08 1998.
  * [49] X. Liu and J. Domke.  Projecting Markov random field parameters for fast mixing.  NIPS’14, pages 1377–1385, Cambridge, MA, USA, 2014. MIT Press.
  * [50] 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.
  * [51] A. McCallum, K. Schultz, and S. Singh.  FACTORIE: Probabilistic programming via imperatively defined factor graphs.  In Advances in Neural Information Processing Systems, volume 22. Curran Associates, Inc., 2009.
  * [52] M. Mitzenmacher and E. Upfal.  Probability and computing: Randomization and probabilistic techniques in algorithms and data analysis.  Cambridge university press, 2017.
  * [53] J. M. Mooij and S. Ong.  libDAI: A free/open source C++ library for discrete approximate inference methods, 2008.
  * [54] R. M. Neal.  Annealed importance sampling.  Statistics and computing, 11(2):125–139, 2001.
  * [55] D. Newman, P. Smyth, M. Welling, and A. Asuncion.  Distributed inference for latent Dirichlet allocation.  In Advances in Neural Information Processing Systems, volume 20. Curran Associates, Inc., 2008.
  * [56] M. Patriarca, A. Chakraborti, and K. Kaski.  Gibbs versus non-gibbs distributions in money dynamics.  Physica A: Statistical Mechanics and its Applications, 340(1):334–339, 2004.  News and Expectations in Thermostatistics.
  * [57] D. Paulin.  Concentration inequalities for Markov chains by Marton couplings and spectral methods.  Electron. J. Probab., 20, 2015.
  * [58] M. Plummer.  JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling.  2003\.
  * [59] I. Porteous, D. Newman, A. Ihler, A. Asuncion, P. Smyth, and M. Welling.  Fast collapsed Gibbs sampling for latent Dirichlet allocation.  In 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.
  * [60] A. Prasad, V. Srinivasan, S. Balakrishnan, and P. Ravikumar.  On learning Ising models under Huber’s contamination model.  Advances in neural information processing systems, 33, 2020.
  * [61] M. Rabinovich, A. Ramdas, M. Jordan, and M. Wainwright.  Function-specific mixing times and concentration away from equilibrium.  Bayesian Analysis, 15, 05 2016.
  * [62] G. O. Roberts and J. S. Rosenthal.  Examples of adaptive MCMC.  Journal of computational and graphical statistics, 18(2):349–367, 2009.
  * [63] J. Rosenthal et al.  Quantitative convergence rates of Markov chains: A simple account.  Electronic Communications in Probability, 7:123–128, 2002.
  * [64] A. Smola and S. Narayanamurthy.  An architecture for parallel topic models.  Proceedings of the VLDB Endowment, 3(1–2):703–710, Sept. 2010\.
  * [65] D. Stefankovic, S. Vempala, and E. Vigoda.  Adaptive simulated annealing: A near-optimal connection between sampling and counting.  In 48th Annual IEEE Symposium on Foundations of Computer Science (FOCS’07), pages 183–193, 2007.
  * [66] A. S. Stordal and A. H. Elsheikh.  Iterative ensemble smoothers in the annealed importance sampling framework.  Advances in Water Resources, 86:231–239, 2015.
  * [67] L. Theis, J. Sohl-Dickstein, and M. Bethge.  Training sparse natural image models with a 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.
  * [68] C. Tosh.  Mixing rates for the alternating Gibbs sampler over restricted Boltzmann machines and friends.  In ICML, 2016.
  * [69] D. Štefankovič, S. Vempala, and E. Vigoda.  Adaptive simulated annealing: A near-optimal connection between sampling and counting.  J. ACM, 56(3), May 2009.
  * [70] R. Zhang and C. De Sa.  Poisson-minibatching for Gibbs sampling with convergence rate guarantees.  In NeurIPS, 2019.