A Faster Approximation Algorithm for the Gibbs Partition Function
Vladimir Kolmogorov
Institute of Science and Technology Austria
vnk@ist.ac.at

arXiv:1608.04223v4 [cs.DS] 27 Dec 2017

Abstract
P
We consider the problem of estimating the partition function Z(β) =
x exp(−βH(x)) of
a Gibbs distribution with a Hamilton H(·), or more precisely the logarithm of the ratio q =
ln Z(0)/Z(β). It has been recently shown how to approximate q with high probability assuming
the existence of an oracle that produces samples from the Gibbs distribution for a given parameter
value in [0, β]. The current best known approach due to Huber [9] uses O(q ln n·[ln q +ln ln n+ε−2 ])
oracle calls on average where ε is the desired accuracy of approximation and H(·) is assumed to lie
in {0} ∪ [1, n]. We improve the complexity to O(q ln n · ε−2 ) oracle calls. We also show that the
same complexity can be achieved if exact oracles are replaced with approximate sampling oracles
2
that are within O( q εln n ) variation distance from exact oracles. Finally, we prove a lower bound of
Ω(q · ε−2 ) oracle calls under a natural model of computation.

1

Introduction

It is known that for large classes of problems, e.g. self-reducible problems [14], there is an intimate
connection between approximate counting and sampling: the ability to solve one problem allows solving
the other one. This paper explores this connection for Gibbs distributions.
Let Ω be some finite set and H(·) be some real-valued function on Ω called a Hamiltonian. The
Gibbs distribution for such a system is a family of distributions {µβ } on Ω parameterized by β, where
µβ (x) =

1
exp(−βH(x))
Z(β)

The normalizing constant Z(β) is called the partition function:
X
Z(β) =
exp(−βH(x))

∀x ∈ Ω

(1)

(2)

x∈Ω

Estimating this function for a given value of β is a widely studied computational problem with applications in many areas. In particular, it is a key computational task in statistical physics. Evaluations of
Z(·) yield estimates of important thermodynamical quantities, such as the free energy. Note, parameter β corresponds to the inverse temperature. A classical example of a Gibbs distribution in physics
is the Ising model.
P
Example 1. Given an undirected graph (V, E), let Ω = {−1, +1}V and H(x) = {i,j}∈E [xi 6= xj ]
where [·] is 1 if its argument is true, and 0 otherwise. Distribution (1) for such a Hamiltonian is called
the Ising model. It is ferromagnetic if β > 0, and antiferromagnetic if β < 0 (although in the latter
case the function H ′ (x) = −H(x) is usually treated as the Hamiltonian). Computing Z(β) exactly is
a #P-complete problem, and is even hard to approximate in the antiferromagnetic case [13].
The problem of counting various combinatorial objects such as proper k-coloring and matchings
in graphs can also be naturally phrased as estimating the partition function.
|V | be the set of all colorings in an undirected graph G = (V, E).
Example 2. Let
P Ω = {1, . . . , k}
Define H(x) = {i,j}∈E [xi = xj ], then Z(+∞) gives the number of proper k-colorings.

1

Example 3. Let Ω be the set of matchings M ⊆ E in an undirected graph G = (V, E). Define
H(M ) = |M |, then Z(0) = |Ω|.
A related problem is that of sampling from the distribution µβ for a given value of β. There is
a vast literature on designing sampling algorithms from Gibbs distributions, see e.g. [17, 19, 8, 6]
or [4] for an overview. For the ferromagnetic Ising model there exists a polynomial-time approximate
sampling algorithm [13] and also an exact sampling algorithm that appears to run efficiently at or above
the critical temperature [18]. Approximate sampling of k-colorings in low-degree graphs is addressed
in [12, 21] (for β = +∞, though techniques are potentially extendable to other values of β), and for
matchings polynomial-time approximate sampling is described in [16, Section 2.3.5].
It is known that the ability to sample can be used for designing a randomized approximation scheme
for estimating the partition function. By definition, it ishan algorithm ithat for a given ε > 0 produces
Q
, Q(1 + ε) with probability at least 3/4.
an estimate Q̂ of the desired quantity Q such that Q̂ ∈ 1+ε
(The value 3/4 is arbitrary: by repeating the algorithm multiple times and taking the median of
the outputs the probability can be boosted to any other constant in (0, 1)). This paper studies the
following question: how many samples are needed to approximate Z(β) with a given accuracy ε?
Formal description To state the complexity of different approaches, we need to introduce several
quantities. First, we assume that H(x) ∈ {0} ∪ [1, n] for any x ∈ Ω where n is a known number.
Non-negativity of the Hamiltonian implies that Z(·) is a decreasing function. Our goal will be to
estimate the ratio Q = Z(βmin )/Z(βmax ) for given values βmin < βmax . Note that computing Z(β) for
some specific value of β is usually an easy task, so this will allow estimating Z(β) for any other β. In
particular, in Examples 1, 2 and 3 we have Z(0) = 2|V | , Z(0) = k|V | and Z(+∞) = 1 respectively.
Let us denote q = log Q, and assume that there exists an oracle that can produce a sample X ∼ µβ
for a given value β ∈ [βmin , βmax ]. When stating asymptotic complexities, we will always assume
that q = Ω(1), n = 1 + Ω(1) and ε = O(1) to simplify the expressions. Bezáková et al. [2] showed
that Q can be estimated using O(q 2 (ln n)2 ) oracle calls in the worst case (for a fixed ε). This was
improved to O(q(ln q + ln n)5 ε−2 ) expected number of calls by Štefankovič et al. [22] and then to
O(q ln n · [ln q + ln ln n + ε−2 ]) by Huber [9].
The first contribution of this paper is to improve the complexity further to O(q ln n · ε−2 ) oracle
calls (on average). This is achieved by a better analysis of the algorithm in [9]. The formal statement
of our result is given in Section 3 as Theorems 5 and 7.
In many applications we only have an access to approximate sampling oracles. Using a standard
coupling argument, in Section 3.1 we show that the same complexity can be achieved with approximate
2
oracles assuming that they are within O( q εln n ) variation distance from exact oracles.
As our final contribution, we prove a lower bound of Ω(q · ε−2 ) oracle calls under a natural model
of computation. The precise statement of the result is given as Theorem 12 in Section 5.
Remark 1. The assumption that H(·) lies in {0} ∪ [1, n] can be relaxed using a standard trick.
Suppose, for example, that H(x) ∈ {hmin , hmin + 1, . . . , hmax } where hmin and hmax are known integers.
Let n = hmax − hmin . We claim that the problem can be solved using O(q ′ ln n · ε−2 ) oracle calls (on
average), where either (i) q ′ = q − (βmax − βmin ) · hmin , or (ii) q ′ = −q + (βmax − βmin ) · hmax .
Indeed, to achieve the first complexity, define new Hamiltonian H ′ (x) = H(x)−hmin . The partition
function for the new problem is Z ′ (β) = eβhmin · Z(β), and so q ′ is as defined in (i). (We use “primes”
to denote all quantities related to the new problem). We have H ′ (x) ∈ {0, 1, . . . , n}, so the algorithm
claimed above can be applied to give an estimate of q ′ and thus of q. Note that distributions µ′β and
µβ coincide, and so sampling from µβ allows to sample from µ′β .
To achieve the second complexity, define H ′ (x) = −H(x) + hmax and also change the bounds:
′
′
= −βmin . There holds Z ′ (β) = e−βhmax · Z(−β), and q ′ is as defined in (ii).
βmin
= −βmax and βmax
′
We again have H (x) ∈ {0, 1, . . . , n}, and distributions µ′β and µ−β coincide. We can now use the
same argument as before.

2

2

Background and preliminaries

We will assume for simplicity that H(·) 6= const. Let us denote z(β) = ln Z(β). It can be easily
checked that
z ′ (β) = EX∼µβ [−H(X)]
Since H(·) is non-negative and non-constant, we have z ′ (β) < 0 for any β and thus z(·) and Z(·) are
strictly decreasing functions. It is also known [23, Proposition 3.1] that function z(·) is convex for any
H(·), and in fact strictly convex if H(·) 6= const.
Next, we discuss previous approaches for estimating Z(βmin )/Z(βmax ), closely following [9].
It is well-known that for given values β1 , β2 an unbiased estimator W of Z(β2 )/Z(β1 ) can be
obtained as follows: first sample X ∼ µβ1 and then set W = exp((β1 − β2 )H(X)). Indeed,
E[W ] =

X exp(−β2 H(x))
X exp(−β1 H(x))
Z(β2 )
· exp((β1 − β2 )H(x)) =
=
Z(β1 )
Z(β1 )
Z(β1 )
x∈Ω

x∈Ω

Applying this estimator directly to (β1 , β2 ) = (βmin , βmax ) or to (β1 , β2 ) = (βmax , βmin ) is problematic
since it usually has a huge relative variance. A standard approach to reduce the relative variance is
via the multistage sampling method of Valleau and Card [20]. First, a sequence βmin = β0 ≤ β1 ≤
. . . ≤ βℓ = βmax is selected; it is called a cooling schedule. We then have
Z(βmin )
Z(β0 ) Z(β1 )
Z(βℓ−1 )
=
·
· ... ·
Z(βmax )
Z(β1 ) Z(β2 )
Z(βℓ )
Throughout the paper we refer to [βi , βi+1 ] as “interval i”, where i ∈ {0, 1, . . . , ℓ − 1}. The ratio
Z(βi )/Z(βi+1 ) for each such interval can be estimated independently as described above, and then
multiplied to give the final estimate. Fishman calls an estimate of this form a product estimator [7].
Its mean and variance are given by the lemma below. In this lemma we use the following notation: if
def E[X 2 ]
Var(X)
X is a random variable then S[X] = (E[X])
2 = (E[X])2 + 1 (the relative variance of X plus 1).
Q
Lemma 1 ([5, page 136]). For P = i Pi where the Pi are independent,
Y
Y
E[P ] =
E[Pi ],
S[P ] =
S[Pi ]
i

i

Using a fixed cooling schedule, Bezáková et al. [2] obtained an approximation algorithm that needs
O(q 2 (ln n)2 ) samples in the worst case (for a fixed ε). Štefankovič et al. [22] asymptotically improved
this to 108 q(ln q + ln n)5 ε−2 samples on average. They used an adaptive cooling schedule where the
values βi depend on the outputs of sampling oracles. A further improvement to O(q ln n·[ln q +ln ln n+
ε−2 ]) was given by Huber [9]. One of the key ideas in [9] was to replace the product estimator with
the paired product estimator, which is described next.

2.1

Paired product estimator

One run of this estimator can be described as follows:
• sample Xi ∼ µβi for each i ∈ [0, ℓ]
• for each interval i ∈ [0, ℓ − 1] compute
Wi = exp(− βi+12−βi H(Xi )),
• compute W =

Q

i Wi

and V =

Q

i

Vi .
3

Vi = exp( βi+12−βi H(Xi+1 ))

An easy calculation (see [9]) shows that
E[Wi ] =

Z(β̄i,i+1 )
,
Z(βi )

E[Vi ] =

where we denoted β̄i,i+1 =

Z(β̄i,i+1 )
,
Z(βi+1 )

βi +βi+1
.
2

S[Wi ] = S[Vi ] =

E[Vi ]/E[Wi ] =

Z(βi )
,
Z(βi+1 )

E[V ]/E[W ] =

Z(βmin )
Z(βmax )

Also,

Z(βi )Z(βi+1 )
,
Z(β̄i,i+1 )2

S[W ] = S[V ] =

Y Z(βi )Z(βi+1 )
i

Z(β̄i,i+1 )2

(3)

Z(βmin )
Although E[V ]/E[W ] = Z(β
= Q, using V /W as the estimator of Q would be a poor choice since
max )
it is biased in general. Instead, [9] uses the following procedure.

Algorithm 1: Paired product estimator. Input: schedule (β0 , . . . , βℓ ), integer r ≥ 1.
1
2

compute r independent samples of (W, V ) as described above
take their sample averages W̄ and V̄ and output Q̂ = V̄ /W̄ as the estimator of Q
The argument from [9] gives the following result.

Lemma 2. Suppose that
S[W ] = S[V ] ≤ 1 + 12 γrε̃ 2

(4)

1
, 1 + ε)) ≥ 1 − γ.
where ε̃ = 1 − (1 + ε)−1/2 = 21 ε + O(ε2 ) and γ > 0. Then P(Q̂/Q ∈ ( 1+ε

Proof. We have E[W̄ ] = E[W ] and Var(W̄ ) = 1r Var(W ), and so S[W̄ ] = 1r (S[W ] − 1) + 1. By
Chebyshev’s inequality, P(|W̄ /E[W̄ ] − 1| ≥ ε̃) ≤ (S[W̄ ] − 1)/ε̃ 2 = 1r (S[W ] − 1)/ε̃ 2 ≤ γ/2. Similarly,
P(|V̄ /E[V̄ ] − 1| ≥ ε̃) ≤ γ/2.
Denote S = W̄ /E[W̄ ] and T = V̄ /E[V̄ ]. The union bound gives P(max{|S − 1|, |T − 1|} ≥ ε̃) ≤ γ.
1
1/2 )
Observe that condition max{|S − 1|, |T − 1|} < ε̃ implies {S, T } ⊂ (1 − ε̃, 1 + ε̃) ⊆ ( (1+ε)
1/2 , (1 + ε)
and thus

Q̂
Q

=

T
S

1
∈ ( 1+ε
, 1 + ε). The claim follows.

Recall that S[W ] = S[V ] is a deterministic function of the schedule (β0 , . . . , βℓ ) (see eq. (3)).
We say that the schedule is good (with respect to fixed constants r and γ) if the resulting quantity
S[W ] = S[V ] satisfies (4). Huber presented in [9] a randomized algorithm that produces a good
schedule with probability at least 0.95 (with respect to r = Θ(ε−2 ) and γ = 0.2). By Lemma 2, the
Q
output Q̂ of the resulting algorithm lies in ( 1+ε
, Q(1 + ε)) with probability at least 0.95 · (1 − γ) > 0.75,
as desired.
Huber’s algorithm for producing schedule (β0 , . . . , βℓ ) is reviewed in the next section. It makes
O(q ln n · [ln q + ln ln n]) calls to the sampling oracle (on average). Then in Section 3 we will describe
how to reduce the number of oracle calls to O(q ln n) while maintaining the desired guarantees.

2.2

TPA method

The algorithm of [9] for producing a schedule is based on the TPA method of Huber and Schott [10, 11].
(The abbreviation stands for the “Tootsie Pop Algorithm”). Let us review the application of the
method to the Gibbs distribution with a non-negative Hamiltonian H(·).
Its key subroutine is procedure TPAstep(β) that for a given constant β produces a random variable
in [β, +∞] as follows:
• sample X ∼ µβ , draw U ∈ [0, 1] uniformly at random, return β−ln U/H(X) (or +∞ if H(X) = 0).
The motivation for this sampling rule comes from the following fact (which we prove here for
completeness).
4

Lemma 3. Consider random variable U = Z(TPAstep(β)). If H(·) is strictly positive (implying that
Z(+∞) = 0) then U has the uniform distribution on [0, Z(β)]. If H(x) = 0 for some x ∈ Ω (implying
that Z(+∞) > 0)) then U has the same distribution as the following random variable U ′ : sample
U ′ ∈ [0, Z(β)] uniformly at random and set U ′ ← max{U ′ , Z(+∞)}.
Proof. It suffices to prove P(TPAstep(β) ≥ α) = Z(α)/Z(β) for any α ∈ [β, +∞). We have
[TPAstep(β) ≥ α] = [ln U/H(X) ≤ β − α] = [ln U ≤ (β − α)H(X)] = [U ≤ exp((β − α)H(X))]
Therefore,
P(TPAstep(β) ≥ α) =
=

X

P(TPAstep(β) ≥ α|X = x)P(X = x)

X

P(U ≤ exp((β − α)H(x))) ·

X

exp((β − α)H(x)) ·

x∈Ω

x∈Ω

=

x∈Ω

=

exp(−βH(x))
Z(β)

exp(−βH(x))
Z(β)

X exp(−αH(x))
Z(α)
=
Z(β)
Z(β)

x∈Ω

Roughly speaking, the TPA method counts how many steps are needed to get from βmin to βmax .
Algorithm 2: One run of TPA. Output: a multiset B of values in the interval [βmin , βmax ].
1
2
3
4

set β0 = βmin , let B be the empty multiset
for i = 1 to +∞ do
sample βi = TPAstep(βi−1 )
if βi ∈ [βmin , βmax ] then add βi to B, otherwise output B and terminate

The output of Algorithm 2 will be denoted as TPA(1), and the union of k independent runs of
def
TPA(1) as TPA(k). For a multiset B we define multiset z(B) = {z(β) | β ∈ B} in a natural way. (Recall
that z(·) is a continuous strictly decreasing function). It is known [10, 11] that z(TPA(k)) is a Poisson
Point Process (PPP) on [z(βmax ), z(βmin )] of rate k, starting from z(βmin ) and going downwards. In
other words, the random variable z = z(TPA(k)) is generated by the following process.
Algorithm 3: Equivalent process for generating z(TPA(k)).
1
2
3
4

set z0 = z(βmin ), let z be the empty multiset
for i = 1 to +∞ do
draw η from the exponential distribution of rate k (and with the mean k1 ), set zi = zi−1 − η
if zi ∈ [z(βmax ), z(βmin )] then add zi to z, otherwise output z and terminate

One way to use the TPA method is to simply count the number of points in TPA(k). Indeed, |TPA(k)|
is distributed according to the Poisson distribution with rate k ·(z(βmin )−z(βmax )) = k ·q, so k1 |TPA(k)|
is an unbiased estimator of q. Unfortunately, obtaining a good estimate of q with this approach
requires a fairly large number of samples, namely O(q 2 ) for a given accuracy and the probability of
failure [10, 11]. A better application of TPA was proposed in [9], where the method was used for
generating a schedule (β0 , . . . , βℓ ) as follows.
Algorithm 4: Generating a schedule (β0 , . . . , βℓ ). Input: integers k, d ≥ 1.
1
2
3

sample B ∼ TPA(k)
sort the values in B and then keep every dth successive value
add values βmin and βmax and output the resulting sequence (β0 , . . . , βℓ ) = (βmin , . . . , βmax )
5

Note that the resulting sequence (z1 , . . . , zℓ−1 )=(z(β1 ), . . . , z(βℓ−1 )) can be described by a process
in Algorithm 3 where η is drawn as the sum of d exponential distributions each of rate k; this is the
gamma (Erlang) distribution with shape parameter d and rate parameter k.
Huber showed in [9] that if d = Θ(ln q + ln ln n) and k = Θ(d ln n) (with appropriate constants)
then Algorithm 4 produces a good schedule with high probability. Since q is unknown in practice,
+ 1 is computed, which is shown to be
[9] uses a two-stage procedure: first an estimate q̂ = 2 · |TPA(5)|
5
an upper bound on q with probability at least 0.99. This estimate is then used for setting d and k.
In the next section we prove that the algorithm has desired guarantees for smaller parameter
values, namely d = Θ(1) and k = Θ(ln n). This allows to reduce the complexity of Algorithm 4 by a
factor of Θ(ln q + ln ln n), and also eliminates the need for a two-stage procedure.

3

Our results

For technical reasons we will need to make the following assumption for line 2 of Algorithm 4: if
β1 , β2 , . . . is the sorted sequence of points in B then the index of the first point to be taken is sampled
uniformly from {1, . . . , d} (and after that the index is always incremented by d).
Denote m = kd and zi = z(βi ) for i ∈ [0, ℓ]. We treat m and d as being fixed, and k = md as their
function. Also let δ = ln S[W ] = ln S[V ]. From (3) we get


X
i+1
+ z(βi+1 ).
(5)
δ=
δi ,
δi = z(βi ) − 2z βi +β
2
i

Since z(·) is convex, we have δi ≥ 0 for all i.

Case I: H(x) ∈ [1, n] for all x ∈ Ω First, let us assume that H(·) does not take value 0. In this
case the proofs become somewhat simpler, and we will get slightly smaller constants.
Huber showed that for d = Θ(ln(q ln n)) the schedule is well-balanced with probability Θ(1), mean1
1
for a constant τ = 34 . (Note that E[zi − zi+1 ] ≈ m
,
ing that all intervals i satisfy zi − zi+1 ≤ τ · m
ignoring boundary effects). It was then proved1 that a well-balanced schedule satisfies δ ≤ τ2 · lnmn ,
leading to condition (4). We improve on this result as follows.
1
,
Choose a constant τ > 0 (to be specified later), and say that interval i is large if zi − zi+1 > τ · m
+
−
and small otherwise. Let δ be the sum of δi over large intervals and δ be the sum of δi over small
intervals (so that δ = δ+ + δ− ). In Section 4.1 we prove the following fact. (Recall that δ+ , δ− are
deterministic functions of the schedule (β0 , . . . , βℓ )).

d) ln n
Lemma 4. There holds δ− ≤ τ2 · lnmn and E[δ+ ] ≤ Γ(d+2,τ
2d · d! · m for the schedule (β0 , . . . , βℓ ) produced
by Algorithm 4 with parameters k = md and d, where Γ(·, ·) is the upper incomplete gamma function:
Z +∞
(with Γ(a, 0) = Γ(a) = (a − 1)!)
ta−1 e−t dt
Γ(a, b) =
b

Using Markov’s inequality, we can now conclude that for any τ + > 0 we have
P(δ+ ≥

τ+
2

·

ln n
m )

≤

1
τ+
2

·

ln n
m

· E[δ+ ] ≤

Γ(d + 2, τ d)
τ + · d · d!
+

d)
ln n
τ +τ
Thus, with probability at least 1 − Γ(d+2,τ
2 · m .
τ + · d · d! Algorithm 4 produces a schedule satisfying δ ≤
0.75
Recall that we want Algorithm 4 to succeed with probability at least ρ = 1−γ
to make the overall
probability of success at least 0.75. (Here γ is the constant from Lemma 2). Let us define function
τρ (d) as follows:
i
h
o
n
Γ(d+2,τ d)
d)
≤
1
−
ρ
=
min
τ
+
τρ (d) = min
τ + τ + | Γ(d+2,τ
+
τ · d · d!
(1−ρ) · d · d!
τ ≥0,τ + >0

1

τ ≥0

More precisely, this is what the argument of [9] would give assuming that H(·) does not take value 0.

6

The table below shows some values of this function for γ = 0.24 and ρ =
with the code of [3]).
d
upper bound on τρ (d)
achieved with τ =

1
9.903
8.645

2
6.052
5.384

4
4.000
3.634

8
2.860
2.653

16
2.197
2.075

32
1.794
1.720

64
1.539
1.492

0.75
1−γ

=

128
1.372
1.342

75
76

(computed

256
1.260
1.241

512
1.184
1.170

We can now formulate our main result for case I.
Theorem 5. Let Q̂ be the estimate given by Algorithm 1 (with parameter r) applied to the schedule
produced by Algorithm 4 (with parameters k = md and d). Suppose that
m≥

τρ (d) · ln n

2 ln 1 + 12 γrε̃ 2

for some γ ∈ (0, 0.25) and ρ =

0.75
1−γ

(6)

Q
, Q(1 + ε)) with probability at least 0.75. The
where ε̃ = 1 − (1 + ε)−1/2 = 21 ε + O(ε2 ). Then Q̂ ∈ ( 1+ε
expected number of oracle calls that this algorithm makes is mq(r + d) + 2r + 1.
In particular, (6) will be satisfied for d = 64, m ≥ 3.6 · ln n and r = 2ε̃ −2 = 8(1 + o(1))ε−2 .

Proof. As we just showed,


P δ≤

τρ (d)
2

·

ln n
m



≥ ρ

(7)


τ (d)
Condition δ ≤ ρ2 · lnmn implies condition δ ≤ ln 1 + 12 γrε̃ 2 (by (6)), which is in turn equivalent to
S[W ] ≤ 1 + 12 γrε̃ 2 . Thus, from Lemma 2 we get


τ (d)
Q
P Q̂ ∈ ( 1+ε
, Q(1 + ε)) | δ ≤ ρ2 · lnmn ≥ 1 − γ
(8)

Multiplying (7) and (8) gives the first claim.
A PPP of rate k on an interval [z(βmax ), z(βmin )] produces k[z(βmin ) − z(βmax )] = mdq points on
average. Thus, Algorithm 4 makes mdq +1 oracle calls on average and produces a sequence (β0 , . . . , βℓ )
with E[ℓ] = mq + 1. Algorithm 1 then makes (ℓ + 1)r oracle calls, i.e. (mq + 2)r calls on average. This
gives the second claim.
Case II: H(x) ∈ {0} ∪ [1, n] for all x ∈ Ω
prove the following fact.

We now consider the general case. In Section 4.2 we

Lemma 6. For any constant λ ∈ (0, 1) there exists a decomposition δ = δ− + δ+ with δ− , δ+ ≥ 0 such
that
τ 2 + ln nλ
Γ(d + 2, τ d) 2 + ln nλ
1
+ ·
E[δ+ ] ≤
·
δ− ≤ ln
1−λ 2
m
2d · d!
m
As in the first case, we conclude from the Markov’s inequality that with probability at least
n
d)
τ +τ + 2+ln λ
1
+
·
1 − Γ(d+2,τ
Algorithm
4
produces
a
schedule
satisfying
δ
≤
ln
+
1−λ
2
m . This leads to
τ · d · d!

Theorem 7. The conclusion of Theorem 5 holds if
m≥

τρ (d) · (2 + ln nλ )



2 ln 1 + 12 γrε̃ 2 (1 − λ)

for some γ ∈ (0, 0.25), ρ =

0.75
1−γ

and λ ∈ (0, 1)

(9)



For example, (9) will be satisfied for d = 64, m ≥ 3.6 · (9 + ln n) and r = 2ε̃ −2 = 8(1 + o(1))ε−2
(where we used γ = 0.24 and λ = e−7 ).

7

3.1

Approximate sampling oracles

So far we assumed that exact sampling oracles µβ are used. For many applications, however, we
only have approximate sampling oracles µ̃β that are sufficiently close to µβ in terms of the variation
distance || · ||T V defined via
||µ̃β − µβ ||T V = max |µ̃β (A) − µβ (A)| =
A⊆Ω

1X
|µ̃β (x) − µβ (x)|.
2
x∈Ω

The analysis can be extended to approximate oracles using a standard trick (see e.g. [22, Remark 5.9]).
Theorem 8. Let Q̂ be the output of the algorithm with parameters d, m, r satisfying the conditions
of Theorem 5 or 7 (depending on whether H(·) ∈ [1, n] or H(·) ∈ {0} ∪ [n]), where exact sampling
κ
oracles µβ are replaced with approximate sampling oracles µ̃β satisfying ||µβ − µ̃β ||T V ≤ mq(r+d)+3r+1
.
Q
Then Q̂ ∈ ( 1+ε
, Q(1 + ε)) with probability at least 0.75 − κ.

As mentioned in the introduction, probability 0.75 − κ can be boosted to any other probability in
(0.5, 1) by repeating the algorithm a constant number of times and taking the median (assuming that
κ is a constant in (0, 0.25)). Alternatively, one can tweak parameters in Theorems 5 and 7 to get the
desired probability directly.
Proof. It is known that there exists a coupling between µβ and µ̃β such that they produce identical
κ
. Let
samples with probability at least 1 − ||µ̃β − µβ ||T V ≥ 1 − δ, where we denoted δ = mq(r+d)+3r+1

A and Ã be the algorithms that use respectively exact and approximate samples, where the k-th call
to µβ in A is coupled with the k-th call to µ̃β̃ in Ã when β = β̃. We say that the k-th call is good if the
produced samples are identical. Note, P[k-th call is good | all previous calls were good] ≥ 1 − δ, since
the conditioning event implies β = β̃. Also, if all calls are good then A and Ã give identical results.
Let N be the number of points inside [z(βmax ), z(βmin )] produced by the call TPA(md) in Algon e−λ
rithm 4. Then N follows the Poisson distribution of rate λ = mdq, i.e. P(N = n) = λ n!
. Algorithm 4
N
makes N + 1 oracle calls, and produces a sequence (β0 , . . . , βℓ ) with ℓ ≤ d + 2. Thus, the total number
of oracle calls is N + 1 + (ℓ + 1)r ≤ N c + 3r + 1 where c = 1 + dr . Denoting µ = λ(1 − δ)c , we can write
P[all calls are good] ≥
=

∞
X

n=0
∞
X

nc+3r+1

P(N = n) · (1−δ)
µn e−µ
n!

n=0

n!

· (1−δ)nc+3r+1
c

· eµ−λ (1−δ)3r+1 = eµ−λ (1−δ)3r+1 = e−λ(1−(1−δ) ) (1−δ)3r+1

n=0
−λ(1−(1−cδ))

≥ e

=

∞
X
λn e−λ

(1−δ)3r+1 ≥ (1−λcδ)(1−δ)3r+1 ≥ 1−λcδ−(3r+1)δ ≥ 1−κ

where we used the facts that (1 − x)c ≥ 1 − cx and e−x ≥ 1 − x for x ≥ 0 and c ≥ 1. Using the union
bound, we obtain the claim of the theorem.

4
4.1

Proofs
Proof of Lemma 4

We will assume that the sequence (β0 , . . . , βℓ ) is strictly increasing (this holds with probability 1).
Accordingly, the sequence (z0 , . . . , zℓ ) is strictly decreasing. The following has been shown in [22, 9].
Lemma 9. For any i ∈ [0, ℓ − 1] there holds δi ≤ zi − zi+1 and also
−z ′ (βi )
≥ exp(2δi /(zi − zi+1 ))
−z ′ (βi+1 )
8

z(β)

z(β)
z0

z0

z1

z1

+
+

z2

z2

z3

z3

β

βmin

−
s(β)
s3

βmax

ηs+

(a)

s2

s1

s0

s

(b)

(c)

Figure 1: (a) z(β) = ln Z(β) is a strictly convex decreasing function. Four dots show a possible
output of Algorithm 4. Here ℓ = 3 and [β0 , βℓ ] = [βmin , βmax ]. (b) Definitions of the sets ∆i (in gray).
Intervals 0 and 1 are assumed to be large, while interval 2 is small. (c) Definition of the variable ηs+
(in the case of a large interval).
Proof. Denote β̄ = (βi + βi+1 )/2 and z̄ = z(β̄), then δi = zi+1 − 2z̄ + zi . Since z(·) is a convex strictly
decreasing function, we have
−z ′ (βi ) ≥

zi − z̄
β̄ − βi

− z ′ (βi+1 ) ≤

z̄ − zi+1
βi+1 − β̄

Since β̄ − βi = βi+1 − β̄, taking the ratio gives the second claim of the lemma:
1
1+λ
zi − z̄
−z ′ (βi )
2 (zi − zi+1 + δi )
=
≥
≥ e2λ
= 1
−z ′ (βi+1 )
z̄ − zi+1
1
−
λ
2 (zi − zi+1 − δi )
δi
i+1
≥ 0 and observed that λ < 1 since 1 − λ = 2 zz̄−z
> 0. The fact that
where we denoted λ = zi −z
i+1
i −zi+1
λ < 1 also gives the first claim of the lemma.

Let us define s(β) = ln[−z ′ (β)] and si = ln[−z ′ (βi )] for i ∈ [0, ℓ], then function s(·) and the
sequence (s0 , . . . , sℓ ) are strictly decreasing. Since z(β) and s(β) are continuous strictly decreasing
functions of β, we can uniquely express z via s and define a continuous strictly increasing function
def
z(s) on the interval S = [sℓ , s0 ] (Fig. 1(b)). Note, with some abuse of notation we use z(·) for two
different functions: one of argument β, and another one of argument s. The exact meaning should
always be clear from the context.
The inequality in the last lemma for an interval i ∈ [0, ℓ − 1] can be rewritten as follows:
2δi ≤ (zi − zi+1 ) · (si − si+1 )

(10)

Equivalently, we have 2δi ≤ Area(∆i ) where ∆i ⊆ [sℓ , s0 ] × [zℓ , z0 ] is the rectangle with the top
right corner at (si , zi ) and the bottom left corner at (si+1 , zi+1 ) (Fig. 1(b)). Let ∆+ be the union of
1
rectangles ∆i corresponding to large intervals i (with |zi − zi+1 | > τ · m
), and ∆− be the union of ∆i
corresponding to small intervals i. Then 2δ+ ≤ Area(∆+ ) and 2δ− ≤ Area(∆− ).
By geometric considerations it should be clear that
Area(∆− ) ≤ max {|zi − zi+1 | : i is small} · |S| ≤ τ ·

1
m

· |S|

Observe that −z ′ (β) = EX∼µβ [H(X)] ∈ [1, n] for any β, and therefore S = [sℓ , s0 ] ⊆ [0, ln n] and so
|S| ≤ ln n. This establishes the first claim of Lemma 4. Next, we focus on proving the second claim.
For a point s ∈ S let ηs be the length of the interval (zi+1 , zi ) into which z(s) falls (or 0, if
1
, and
z(s) ∈ {zℓ , . . . , z0 }). Also let ηs+ = ψ[ηs ] where ψ[·] is the following function: ψ[a] = a if a > τ · m
9

ψ[a] = 0 otherwise. Thus, if z(s) ∈ (zi+1 , zi ) for some large interval i then ηs+ = zi − zi+1 (Fig. 1(c)),
otherwise ηs+ = 0. We have
Z
ηs+ ds
Area(∆+ ) =
S

The linearity of expectation gives
+

+

2E[δ ] ≤ E(Area(∆ )) =

Z

S

E[ηs+ ]ds ≤ max E[ηs+ ] · |S|
s∈S

(11)

Now let X0 , X1 , X2 , . . . be a Poisson process on [0, +∞) and X−1 , X−2 , . . . be a Poisson process on
(−∞, 0] (both with rate k). Thus, Xi = ξ0 + . . . + ξi for i ≥ 0 and Xi = −ξ−1 − ξ−2 − . . . − ξi for
i ≤ −1, where ξj are i.i.d. variables from the exponential distribution of rate k. By the superposition
theorem for Poisson processes [15, page 16], bidirectional sequence X = . . . , X−2 , X−1 , X0 , X1 , X2 , . . .
is a Poisson process on (−∞, +∞) (again with rate k), and in particular it is translation-invariant.
Let Y = . . . , Y−2 , Y−1 , Y0 , Y1 , Y2 , . . . be the following process: draw an integer c ∈ {0, . . . , d − 1}
uniformly at random and then set Yi = Xdi+c for each i. It can be seen that Y models the output
(β0 , . . . , βℓ ) of Algorithm 4 as follows: take the sequence z(βmin )− Y0 , z(βmin )− Y1 , z(βmin )− Y2 , . . ., restrict to [z(βmax ), z(βmin )] and append z(βmin ) and z(βmax ). Then the resulting sequence has the same
distribution as (z(β0 ), . . . , z(βℓ )). We assume below that (β0 , . . . , βℓ ) is generated by this procedure.
For a point a ∈ R let θa be the length of the interval (Yi , Yi+1 ) into which a falls (or 0, if no such
interval exists). Note, the distribution of random variable θa does not depend on a (since process Y is
translation-invariant). We also denote θa+ = ψ[θa ], and let θ and θ + = ψ[θ] be random variables with
the same distributions as θa and θa+ , respectively (for any fixed a). Clearly, for each s ∈ [sℓ , s0 ] we have
ηs ≤ θa and ηs+ ≤ θa+ for a suitably chosen a, namely, a = z(βmin ) − z(s). (Note, if z(s) ∈ (zℓ−1 , z1 )
then ηs = θa and ηs+ = θa+ , but at the boundaries the inequalities may be strict). We thus have
E[ηs+ ] ≤ E[θ + ]

(12)

Lemma 10. Variable θ has the gamma (Erlang) distribution with shape parameter d + 1 and rate k,
whose probability density is f (t) = kd+1 td e−kt /d! for t ≥ 0.
Proof. We prove this fact for variable θa with a = 0. We know that Y−1 = Xc−d ≤ 0 and Y0 = Xc ≥ 0,
so θ0 = Y0 − Y−1 (with probability 1). By construction, Xc − Xc−d = ξc−d + ξc−d+1 + . . . + ξc , i.e. θ0
is a sum of d + 1 i.i.d. exponential random variables each of rate k. This implies the claim.
Remark 2. It may seem counterintuitive that all intervals ζi = Yi − Yi−1 are distributed as a sum of
d exponential random variables with the exception of i = 0, in which case it is a sum of d + 1 variables
(even though Y is translation-invariant). This can be viewed as an instance of the “inspection paradox”, discussed e.g. at [1]. Below we describe an alternative approach, which may help to understand
this phenomenon.
Let ζ be a sum of d exponential random variables each of rate k and g(·) be the probability density
of ζ. Then the following (non-rigorous) argument shows that the probability density of θ is tg(t)/E[ζ]
(after which a simple calculation would prove the claim).
Let L be some large number. Since the distribution of θa does not depend on a, we can define θ
as the output of the following process: sample Y, sample a ∈ [0, L] uniformly at random, and then
set θ = θa (i.e. the length of the interval in Y into which a falls). Let us compute the probability
that θ ∈ [t, t + dt]. Process Y will have L/E[ζ] intervals in [0, L] on average, and out of those
(L/E[ζ]) · (g(t)dt) intervals will have length in the range [t, t + dt]. The combined length of such
intervals is (L/E[ζ]) · (g(t)dt) · t. Thus, point a will fall into one of those intervals with probability
(L/E[ζ]) · (g(t)dt) · t/L = (tg(t)/E[ζ])dt. Therefore, the density of θ is tg(t)dt/E[ζ].

10

Recall that θ + = θ if θ > τ /m, and θ + = 0 otherwise. Lemma 10 now gives
+

E[θ ] =

+∞
Z

tf (t)dt =

1
k · d!

+∞
Z

kd+1 td+1 e−kt
dt =
d!

τ /m

τ /m

=

+∞
Z

+∞
Z

τ d/k

ud+1 e−u du =

τd

(kt)d+1 e−(kt)
d(kt)
k · d!

Γ(d + 2, τ d)
Γ(d + 2, τ d)
=
k · d!
md · d!

Combining this with (11) and (12) and observing again that |S| ≤ ln n finally gives the second
claim of Lemma 4.

4.2

Proof of Lemma 6

We will use the same notation as in the previous section. Since H(·) can now take value 0, we have
−z ′ (β) = EX∼µβ [H(X)] ∈ [0, n] and so [sℓ , s0 ] ⊆ [−∞, ln n] (instead of [sℓ , s0 ] ⊆ [0, ln n], as in the
previous section). We will deal with small values of s(β) exactly as in [9].
Recall that z ′ (β) is a strictly increasing function of β. Let β̂ be the unique value with z ′ (β̂) = −λ.
(If it does not exist, then we take β̂ ∈ {−∞, +∞} using the natural rule). Denote ẑ = z(β̂) and
ŝ = ln[−z ′ (β̂)]. Now introduce the following terminology for an interval i ∈ [0, ℓ − 1]:
• interval i is steep if βi+1 ≤ β̂, or equivalently si+1 ≥ ŝ;
• interval i is flat if βi ≥ β̂, or equivalently si ≤ ŝ;
• interval i is crossing if β̂ ∈ (βi , βi+1 ), or equivalently ŝ ∈ (si+1 , si ).
If steep intervals exist then β̂ ≥ βmin and z ′ (β̂) ≤ −λ. (The inequality may be strict if β̂ = +∞).
We thus have [si+1 , si ] ⊆ [ŝ, s0 ] ⊆ [ln λ, ln n] for all steep intervals i. The argument from the previous
section gives that


X
X
Γ(d + 2, τ d) ln nλ
τ ln nλ
δi  ≤
E
·
δi ≤ ·
2 m
2d · d!
m
i: i is steep and large

i: i is steep and small

(We just need to assume that βmax was replaced with min{βmax , β̂}, then we would have S = [sℓ , s0 ] ⊆
[ln λ, ln n] and |S| ≤ ln nλ instead of |S| ≤ ln n, the rest is the same as in the previous section).
Let us now consider flat intervals. The argument from [9] gives the following fact.
1
Lemma 11. The sum of δi over flat intervals i is at most ln 1−λ
.

Proof. Assume that flat intervals exist, then β̂ ≤ βmax and z ′ (β̂) ≥ −λ. (The inequality may be strict
if β̂ = −∞). Denote Ω0 = {x ∈ Ω | H(x) = 0} and Ω+ = {x ∈ Ω | H(x) ≥ 1}, then Ω = Ω0 ∪ Ω+ and
EX∼µβ̂ [H(X)] =

P

x∈Ω+

H(x)e−β̂H(x)
Z(β̂)

≥

P

x∈Ω+

e−β̂H(x)

Z(β̂)

On the other hand, EX∼µβ̂ [H(X)] = −z ′ (β̂) ≤ λ and so

=1−

Z(βmax )
Z(β̂)

P

x∈Ω0

e−β̂H(x)

Z(β̂)

≥1−

Z(βmax )
Z(β̂)

1
≥ 1 − λ and z(β̂) − z(βmax ) ≤ ln 1−λ
.

For all flat intervals i we have [zi+1 , zi ] ⊆ [z(βmax ), z(β̂)] and also δi ≤ zi − zi+1 . This gives the claim
of the lemma.

11

It remains to consider crossing intervals. Let us define values δc− and δc+ as follows. If there are no
crossing intervals then δc− = δc+ = 0. Otherwise let i be the unique crossing interval; if i is small then
τ
(since
set (δc− , δc+ ) = (δi , 0), and if i is large then set (δc− , δc+ ) = (0, δi ). In all cases we have δc− ≤ m
Γ(d+2,τ d)
+
+
δi ≤ zi − zi+1 ). Also, E[δc ] ≤ E[ψ(zi − zi+1 )] ≤ E[θ ] ≤ md · d! where function ψ(·) and random
variable θ + were defined in the previous section.
We can finally prove Lemma 6. Define δ− as δc− plus the sum of δi over small steep intervals i
and flat intervals i. Define δ+ as δc+ plus the sum of δi over large steep intervals i. By collecting
inequalities above we obtain the desired claim.

5

Lower bound

In this section we establish a lower bound on the number of calls to the sampling oracles for estimating
Z(βmin )
q = ln Z(β
. First, we describe our model of computation and the set of instances that we allow.
max )
We assume that the estimation algorithm only receives values H(x) from the sampling oracle,
and not individual states x ∈ Ω. This means that an instance can be defined by counts ch = |{x ∈
Ω | H(x)P= h}| for values h in the range of H; these counts uniquely specify the partition function
Z(β) = h ch e−βh and the distribution of sampling oracle outputs for a given β. We will thus view
an instance as a triplet Γ = (c[Γ], βmin [Γ], βmax [Γ]) where c[Γ] : R → Z≥0 is a function with a finite
non-empty support. When the instance is clear from the context, we will omit the square brackets
and write simply Γ = (c, βmin , βmax ). For a value h ∈ supp(c) let ψ(β, h | Γ) be the probability that
the sampling oracle returns value h when queried at β in instance Γ:
ψ(β, h | Γ) = ch e−βh /Z(β)
For a finite subset H ⊆ R let I(H) be the set of instances Γ = (c, 0, βmax ) satisfying supp(c) ⊆ H.
Z(0)
.
Also for a subset Q ⊆ R let I(H, Q) = {Γ ∈ I(H) | q ∗ (Γ) ∈ Q}, where we denoted q ∗ (Γ) = ln Z(β
max )
An estimation algorithm A applied to instance Γ = (c, 0, βmax ) ∈ I(H) is assumed to have the
following form. At step i (for i = 1, 2, . . .) it does one of the following two actions:
1. Call the samping oracle for some value βi ∈ R. The oracle then returns a random variable hi ∈ H
with P(hi = h) = ψ(βi , h | Γ) for each h ∈ H.
2. Output some estimate q̂ and terminate.
The i-th action is a random variable that can depend only on the set supp(c), values βmin , βmax , and
on the previously observed sequence (β1 , h1 ), . . . , (βi−1 , hi−1 ). The output q̂ of the algorithm will be
denoted as q A (Γ), and the expected number of calls to the sampling oracle as T A (Γ).
We say that algorithm A is an (ε, δ)-estimator for instance Γ if P[|q A (Γ) − q ∗ (Γ)| > ε] < δ. We
can now formulate our main theorem.
Theorem 12. There exist positive numbers qmin , nmin , c1 , c2 , c3 such that the following holds for all
1
q ≥ qmin , n ≥ nmin
l withm n ∈ Z, ε ∈ (0, c1 q), δ ∈ (0, 4 ).
√
c2 q
n

and Hnm = {h ∈ [1, n] : mh ∈ Z}. Suppose that A is an (ε, δ)-estimator for
h
i
i
4q
m , 2q , 4q ) such that T A (Γ) ≥
,
).
Then
there
exists
instance
Γ
∈
I(H
all instances in I(Hnm , 2q
n
3 3
3 3
Denote m =

c3 qε−2 ln δ−1 .

5.1

h

Proof of Theorem 12

The proof will be based on the following result. For brevity, we use notation a ± b to denote the closed
interval [a − b, a + b].
12

Lemma 13. Suppose that A is an (ε, δ)-estimator for instances Γ ∈ I(H, {q}) and Γ1 , . . . , Γd ∈
I(H, R \ (q ± 2ε)), where βmax [Γi ] = βmax [Γ] and supp(c[Γi ]) = supp(c[Γ]) for i ∈ [d]. Suppose that
Y ψ(β, h | Γi )
≥γ
ψ(β, h | Γ)

∀β ∈ R, h ∈ supp(c[Γ])

i∈[d]

for some constant γ ∈ (0, 1). Then T A (Γ) ≥

(1−δ′ −δ)d ln(δ′ /δ)
ln(1/γ)

(13)

for any constant δ′ ∈ [δ, 1 − δ].

Proof. The run of the algorithm can be described by a random variable X = ((β1 , h1 ), . . . , (βt , ht ), q̂),
where t is the number of oracle calls (possibly infinite, in which case q̂ is undefined). Let X be the set
all possible runs, and PA (· | Γ̃) be the probability measure of this random variable conditioned on Γ̃
being the input instance. The structure of the algorithm implies that this measure can be decomposed
as follows:
dPA (x | Γ̃) = ψ(x | Γ̃) dµA (x)

∀x ∈ {((β1 , h1 ), . . . , (βt , ht ), q̂) ∈ X | t is finite}

(14)

where µA (·) is some measure on X that depends only on the algorithm A, and function ψ(·) is defined
via
Y
ψ((β1 , h1 ), . . . , (βt , ht ), q̂ | Γ̃) =
ψ(βi , hi | Γ̃)
i∈[t]

Denote τ =

d ln(δ′ /δ)
ln(1/γ) ,

and define the following subsets of X :

X ∗ = {((β1 , h1 ), . . . , (βt , ht ), q̂) ∈ X | t ≤ τ and q̂ ∈ q ± ε }
X ′ = {((β1 , h1 ), . . . , (βt , ht ), q̂) ∈ X | t > τ }

/ q±ε }
X ′′ = {((β1 , h1 ), . . . , (βt , ht ), q̂) ∈ X | t ≤ τ and q̂ ∈
Suppose the claim of Lemma 13 is false, i.e. T A (Γ) ≤ (1 − δ′ − δ) · τ . We have T A (Γ) ≥ PA (X ′ | Γ) · τ ,
and therefore
PA (X ′ | Γ) ≤ 1 − δ − δ′

(15)

Since A is a (ε, δ)-estimator for instances Γ, Γ1 , . . . , Γd , we have
PA (X ′′ | Γ)
A

< δ

∗

P (X | Γi ) < δ

(16)
∀i ∈ [d]

(17)

Set X is a disjoint union of X ∗ , X ′ , X ′′ , therefore PA (X ∗ | Γ) = 1 − PA (X ′ | Γ) − PA (X ′′ | Γ) >
1 − (1 − δ − δ′ ) − δ = δ′ . Combining this with (17) gives
1X A ∗
P (X | Γi ) <
d
i∈[d]

δ A ∗
P (X | Γ)
δ′

(18)

Assumption (13) of the lemma gives that
Y ψ(x | Γi )
≥ γt ≥ γτ
ψ(x | Γ)

i∈[d]

∀x = ((β1 , h1 ), . . . , (βt , ht ), q̂) ∈ X ∗

(19)

We can now write
1/d

Y
X
δ
1
ψ(x | Γi )
≥ γ τ /d ψ(x | Γ) = ′ ψ(x | Γ)
ψ(x | Γi ) ≥ 
d
δ
i∈[d]

i∈[d]

13

∀x ∈ X ∗

(20)

where the first inequality is a relation between arithmetic and geometric means of non-negative numbers, and the second inequality follows from (19). We can write
Z
XZ
(b) δ
1X A ∗
(a) 1
(c) δ
ψ(x | Γi ) dµA (x) ≥ ′
P (X | Γi ) =
ψ(x | Γ) dµA (x) = ′ PA (X ∗ | Γ)
d
d
δ X∗
δ
X∗
i∈[d]

i∈[d]

where (a,c) follow from (14) and (b) follows from (20). We obtained a contradiction to (18).
Recall that by definition coefficients of instances should be non-negative integers. When using
Lemma 13, we can relax this requirement to non-negative rationals (since multiplying coefficients by
a constant does not affect quantities in Lemma 13) and further to non-negative reals (since they can
be approximated by rationals with an arbitrary precision).
We will use Lemma 13 with d = 2 and three instances Γ, Γ+ , Γ− . First, we will describe the
construction of Γ+ and Γ− given an instance Γ. After stating some properties of this construction, we
will define the instance Γ.
Instances Γ+ and Γ−
Suppose that Γ = (c, 0, βmax ) ∈ I(H). We set Γ+ = (c+ , 0, βmax ) and
−
Γ− = (c , 0, βmax ) where functions c+ , c− are given by
hν
c+
h = ch · e ,

−hν
c−
h = ch · e

∀h ∈ R

where ν > 0 is some constant. Let Z(·), Z+ (·), Z− (·) be the partition functions corresponding to Γ,
Γ+ , Γ− , respectively. One can check that
X
Z(β) =
ch e−βh
Z+ (β) = Z(β − ν)
Z− (β) = Z(β + ν)
h∈supp(c)

Denote z(β) = ln Z(β) and zdiff (β) = z(β) − z(βmax + β). Then
q = q ∗ [Γ] = zdiff (0)

q ∗ [Γ+ ] = zdiff (−ν)

q ∗ [Γ− ] = zdiff (ν)

Condition Γ+ , Γ− ∈ I(H, R \ (q ± 2ε)) can thus be written as follows:
|zdiff (±ν) − zdiff (0)| > 2ε

(21)

Condition (13) after cancellations becomes
Z 2 (β)
≥γ
Z(β − ν)Z(β + ν)

∀β ∈ R

or equivalently
z(β − ν) − 2z(β) + z(β + ν) ≤ ln

1
γ

∀β ∈ R

(22)

Let us define the following quantities; note that they depend only on instance Γ:
′
ρ = |zdiff
(0)|

κ = sup z ′′ (β)

(23)

β∈R
2

ρ
).
Lemma 14. Let Γ be an instance with values q = q ∗ (Γ), ρ, κ as described above. Fix ε ∈ (0, 10κ
Suppose that algorithm A is an (ε, δ)-estimator for all instances in I(H, q ± 4ε). Then for any constant
δ′ ∈ [δ, 1 − δ] we have
2(1 − δ′ − δ)ρ2 ln(δ′ /δ)
T A (Γ) ≥
9κε2

14

Proof. Non-negativity of function c implies that function z(·) is convex, and so z ′′ (β) ∈ [0, κ] for all
β ∈ R. Define ν = 3ε/ρ. For β = ±ν we can write
(a)

′
′′
|zdiff (β) − zdiff (0)| = zdiff
(0)β + zdiff
(β̃)

β2
2

(b)

∈ [|ρβ| − κβ 2 , |ρβ| + κβ 2 ]

where in (a) we used Taylor’s theorem with the Lagrange form of the remainder (here β̃ ∈ R), and
′′
(β̃)| = |z ′′ (β̃) − z ′′ (βmax + β̃)| ≤ 2κ. Observing that |ρβ| = 3ε
in (b) we used the fact that |zdiff
9κε
and κβ 2 = ε · ρ2 < ε, we get |zdiff (β) − zdiff (0)| ∈ (2ε, 4ε). Thus, condition (21) holds, and
Γ+ , Γ− ∈ I(H, q ± 4ε).
Denote f (β) = z(β) − z(β − ν). Using twice the mean value theorem, we get
˜ 2 ≤ κν 2
z(β − ν) − 2z(β) + z(β + ν) = f (β + ν) − f (β) = f ′ (β̃)ν = [z ′ (β̃) − z ′ (β̃ − ν)]ν = z ′′ (β̃)ν
˜
where β̃, β̃ ∈ R. Thus, condition (22) will be satisfied if we set γ ∈ (0, 1) so that ln γ1 = κν 2 =
Lemma 14 now follows from Lemma 13.

9κε2
ρ2 .

Instance Γ We now need to construct instance Γ such that q = q ∗ [Γ] is close to a given value q̄, and
2
the ratio ρκ is large. We will use an instance with the following partition function:
−β

Z(β) = e

N
Y

(ak + e−β/m )

(24)

k=1

where N is some
[m(n − 1)] and a1 , . . . , aN are non-negative numbers. Expanding terms
P integer in
−βh
ch e
for some coefficients ch ≥ 0, so this is indeed a valid definition of an
yields Z(β) = h∈Hm
n
instance Γ ∈ I(Hnm ). In Section 5.2 we prove the following fact.
ρ2
κ

> ( N4 −1)2 .
q
This will imply Theorem 12. Indeed, let q̄ be the value chosen in Theorem 12. Set N̂ = ln22 q̄
l m
and N = N̂ . Note that

Lemma 15. There exist values a1 , . . . , aN , βmax > 0 such that q =

N̂
≤
m(n − 1)
Thus, setting c2 >

q

q

√

2
ln 2 q̄

c2 q̄
n (n

− 1)

= const ·

n
n−1

with

ln 2 2
2 N ±O(mN )

const =

q

2
ln 2

and

/ c2

2
ln 2

will ensure that N ∈ [m(n − 1)] for sufficiently large n.

 
√
We have q = ln22 N 2 ± O(mN ) = ln22 N̂ 2 ± O(mN̂ ) = q̄ ± O(m q̄) = q̄ 1 ± O √mq̄ . Recalling
l √ m
h
i
c q
5q̄
,
that m = 2n , we conclude that q ∈ 3q̄
if q̄, n are sufficiently large. Furthermore, we have
4 4
ρ2
κ

1
> ( N4 − 1)2 > 61 q̄ if q̄ is sufficiently large (note that 61 < 8 ln
2 ).
1
1
that
the
preconditions of Theorem 12 hold. It
We set c1 = 60 , so that ε ∈ (0, 60 q̄). Now hsuppose
i
2

ρ
can be checked that ε ∈ (0, 10κ
) and q ± 4ε ⊆

Setting

δ′

=

1
2

and recalling that δ ∈ (0,

T A (Γ) ≥

2(1 −

1
2

1
4 ),

2q̄ 4q̄
3 , 3

, so the preconditions of Lemma 14 hold as well.

we obtain the desired result:

2
−1
− δ) ln( 12 /δ) ρ2
1 − ln
ln δ−1 − ln 2 q̄
ln 4 q̄ ln δ
≥
≥
·
·
·
9ε2
κ
18ε2
6
18 · 6
ε2

15

5.2

Proof of Lemma 15

Denote u = u(β) = e−β/m and η = u(βmax ) = e−βmax /m . (The choice of η ∈ (0, 1) will be specified
later). We can write
z(β) = −β +

N
X

′

ln(ak + u)

z (β) = −

k=1

N
X

k=N

u
m(ak + u)

′′

z (β) =

N
X
k=1

ak u
2
m (ak +

u)2

(25)

N

1 X ak + 1
+
ln
η
ak + η
k=1

N 
1 X
η
1
′
′
ρ = |z (0) − z (βmax )| =
−
m
ak + 1 ak + η

(26)

q = z(0) − z(βmax ) = m ln

(27)

k=1

As for κ = maxβ∈R z ′′ (β), we will use the following bound.
Lemma 16. Suppose that a1 ≥ a2 ≥ . . . ≥ aN > 0. Then κ ≤ maxr∈[N −1] κr where we denoted
κr

" r
#
N
X
ak
1 X ar
+
=
m2
ak
ar+1
k=r+1

k=1

Proof. We need to show that
N
X
k=1

ak u
≤
m2 (ak + u)2

max κr

r∈[N −1]

∀u ∈ (0, +∞)
ak u
(ak +u)2

is increasing on [0, ak ] and
PN
decreasing on [ak , +∞) (with the maximum at u = ak ). Therefore, function g(u) =
k=1 gk (u)
attains a maximum at [aN , a1 ]. We can thus assume w.l.o.g. that u ∈ [aN , a1 ].
Let r ∈ [N − 1] be an index such that u ∈ [ar+1 , ar ]. For k ∈ [1, r] we have gk (u) ≤ gk (ar ) =
ar /ak
ak /ar+1
ak
. By summing
≤ aakr , and for k ∈ [r + 1, N ] we have gk (u) ≤ gk (ar+1 ) = (1+a
2 ≤ a
(1+ar /ak )2
r+1
k /ar+1 )
2
these inequalities we get that g(u) ≤ m κr .
By taking the derivative one can check that function gk (u) =

We can now prove Lemma 15. Define ak = 21−k and η = 21−N . For each k ∈ [N ] we have
+1
= ln a1k +ln(1+ak ) ≤ ln a1k +ak = (k−1) ln 2+21−k ,
≥ ln 2a1k = (k−2) ln 2 and ln aakk +η
< ln aka+1
ln
k
therefore
ak +1
ak +η



N
X
N
(k − 2) ln 2 =
m+
(N − 1) ln 2 − N ln 2
2
k=1


N h
i
X
N
1−k
(N − 1) ln 2 + 2
q < m(N − 1) ln 2 +
(k − 1) ln 2 + 2
<
m+
2
q > m(N − 1) ln 2 +

k=1

The following inequalities imply the last two claims of Lemma 15:
ρ >
κr <






N 
1
η
2N − 1
1 N
1 X
1 N
−
− N −1
−2
>
=
m
1 + 1 ak
m 2
2
m 2
k=1
" r
#


+∞
X
1 X ar
ak
4
1 2r − 1
+
+2 <
=
2
2
r−1
m
ak
ar+1
m
2
m2
k=1

k=r+1

16

∀r ∈ [N − 1]

Acknowledgements
I thank Laszlo Erdös and Alexander Zimin for useful discussions. In particular, the link [1] provided
by Alexander helped with the argument in Section 4.1. The author is supported by the European
Research Council under the European Unions Seventh Framework Programme (FP7/2007-2013)/ERC
grant agreement no 616160.

References
[1] http://math.stackexchange.com/questions/74454/paradox-of-a-poisson-process-on-mathbb-r.
[2] I. Bezáková, D. Štefankovič, V. V. Vazirani, and E. Vigoda. Accelerating simulated annealing for the
permanent and combinatorial counting problems. SIAM J. Comput., 37:1429–1454, 2008.
[3] G. P. Bhattacharjee. Algorithm AS 32: The incomplete gamma integral. Journal of the Royal Statistical
Society. Series C (Applied Statistics), 19(3):285–287, 1970.
[4] Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng, editors. Handbook of Markov chain
Monte Carlo. Chapman & Hall/CRC, 2011.
[5] M. Dyer and A. Frieze. Computing the volume of convex bodies: A case where randomness provably helps.
In Proceedings of AMS Symposium on Probabilistic Combinatorics and Its Applications 44, pages 123–170,
1991.
[6] J. A. Fill and M. L. Huber. Perfect simulation of Vervaat perpetuities. Electron. J. Probab., 15:96–109,
2010.
[7] G. S. Fishman. Choosing sample path length and number of sample paths when starting in the steady
state. Oper. Res. Lett., 16:209–219, 1994.
[8] Mark Huber. Perfect sampling using bounding chains. Annals of Applied Probability, 14(2):734–753, 2004.
[9] Mark Huber. Approximation algorithms for the normalizing constant of Gibbs distributions. The Annals
of Applied Probability, 25(2):974–985, 2015.
[10] Mark Huber and Sarah Schott. Using TPA for Bayesian inference. Bayesian Statistics 9, pages 257–282,
2010.
[11] Mark Huber and Sarah Schott. Random construction of interpolating sets for high-dimensional integration.
J. Appl. Prob., 51:92–105, 2014.
[12] M. Jerrum. A very simple algorithm for estimating the number of k-colourings of a low-degree graph.
Random Structures and Algorithms, 7:157–165, 1995.
[13] M. Jerrum and A. Sinclair. Polynomial-time approximation algorithms for the Ising model. SIAM J.
Comput., 22:1087–1116, 1993.
[14] Mark R. Jerrum, Leslie G. Valiant, and Vijay V. Vazirani. Random generation of combinatorial structures
from a uniform distribution. Theoret. Comput. Sci., 43(2-3):169–188, 1986.
[15] J. F. C. Kingman. Poisson Processes. Clarendon Press, 1992.
[16] James Matthews. Markov Chains for Sampling Matchings. PhD thesis, University of Edinburgh, School of
Informatics, 2008.
[17] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of state
calculation by fast computing machines. J. Chem. Phys., 21:1087–1092, 1953.
[18] James G. Propp and David B. Wilson. Exact sampling with coupled Markov chains and applications to
statistical mechanics. Random Structures and Algorithms, 9(1-2):223–252, 1996.
[19] Robert H. Swendsen and Jian-Sheng Wang. Replica Monte Carlo simulation of spin-glasses. Phys. Rev.
Lett., 57(21):2607–2609, 1986.
[20] J. P. Valleau and D. N. Card. Monte Carlo estimation of the free energy by multistage sampling. J. Chem.
Phys., 57:5457–5462, 1972.
[21] E. Vigoda. Improved bounds for sampling colorings. In FOCS, pages 51–59, 1999.
[22] D. Štefankovič, S. Vempala, and E. Vigoda. Adaptive simulated annealing: A near-optimal connection
between sampling and counting. J. of the ACM, 56(3):1–36, 2009.
[23] M. J. Wainwright and M. I. Jordan. Graphical models, exponential families, and variational inference.
Foundations and Trends in Machine Learning, 1(1-2):1–305, December 2008.

17

