Convergence of linear programming hierarchies for Gibbs states of spin systems
Authors: Hamza Fawzi
Omar Fawzi
Venue: Accepted by TMLR ()

Abstract: We consider the problem of computing expectation values of local functions under the Gibbs distribution of a spin system. In particular, we study two families of linear programming hierarchies for this problem. The first hierarchy imposes local spin flip equalities and has been considered in the bootstrap literature in high energy physics. For this hierarchy, we prove fast convergence under a spatial mixing (decay of correlations) condition. This condition is satisfied for example above the critical temperature for Ising models on a d-dimensional grid. The second hierarchy is based on a Markov chain having the Gibbs state as a fixed point and has been studied in the optimization literature and more recently in the bootstrap literature. For this hierarchy, we prove fast convergence provided the Markov chain mixes rapidly. Both hierarchies lead to an ε-approximation for local expectation values using a linear program of size quasi-polynomial in n/ε, where n is the total number of sites, provided the interactions can be embedded in a d-dimensional grid with constant d. Compared to standard Monte Carlo methods, an advantage of this approach is that it always (i.e., for any system) outputs rigorous upper and lower bounds on the expectation value of interest, without needing an a priori analysis of the convergence speed.

Published in Transactions on Machine Learning Research (11/2025)
Convergence of linear programming hierarchies for Gibbs states of spin systems
Hamza Fawzi and Omar Fawzi²
1AMTP, University of Cambridge, United Kingdom Inria, ENS de Lajon, UCBL, LIP, France
Reviewed on OpenReview: https: openreview. net/ pormfo id- mc1d
Abstract
We consider the problem of computing expectation values of local functions under the Gibbs distribution a spin system. In particular, we study two families of linear programming hicrarchies for this problem. The first hicrarchy imposes local spin flip equalitios and has been considered in the bootstrap literature in high energy physics. For this hierarchy, prove fast convergence under spatial mixing (decay of correlations) condition. This condition is satisfied for example above the critical temperature for Ising models dimensional grid. The second hierarchy is based on Markov chain having the Gibbs state as a fixed point and has been studied in the optimization literature and more recently in the bootstrap literature. For this hicrarchy, we prove fast convergence provided the Markov chain mixes rapidly. Both hierarchies lead to an s-approximation for local expectation values using a linear program of size quasi-polynomial in n/e, where n is the total number of sites, provided the interactions d. Compared to standard Monte Carlo methods, an advantage of this approach is that it always (i.e., for any system) outputs rigorous upper and lower bounds on the expectation value of interest, without needing an a priori analysis of the convergence speed.
d-
be embedded in d-dimensional grid with constant
1 Introduction
Consider a Gibbs distribution of the form
m(x) Hx((x
ex
(1)
where H is the potential function, Z En H(x), -((x) and S is the state space assumed finite. There are two fundamental computational questions that are associated to any such Gibbs distribution: (i) sampling from nt, and (ii) computing expectation values, i.e., to compute = Eres (xx) where f + R is a given function that we also call obscrvable. the vocabulary of graphical models, a Gibbs distribution is called a Markov random field and question (ii) one of the typical computational inference problems studied in the literature on graphical models MacKay (2003): Wainwright al. (2008).
An algorithm that generates samples from a be used to construct randomized algorithm to compute expectation values. More precisely, if we have an algorithm that runs in time T(e) and produces sample from a distribution y whose total variation distance from μ is then we have an algorithm that approximates p(f) to within & ++6 with running time proportional to T(s)/6² with probability, say, 2/3, assuming that f takes values in the interval []1 Existing algorithms for (i), and hence (ii), are mostly based on Markov Chain Monte Carlo (MCMC). The idea to run a Markov chain whose stationary distribution is μ until approximate convergence. The running time T(s) is then proportional to the mixing time of the Markov
The 1/62 factior comes from taking independent samples and ensuring variance O(82).

---

Published in Transactions on Machine Learning Research (11/2025)
chain. Despite its immense success, one drawback of most MCMC methods is the lack of simple way to certify convergence. 2 In general, one needs to analyze the mixing time of the Markov chain beforehand, which is hard in many cases.
In this paper we focus on the problem of computing expectation values u(f) and we are interested in certified algorithms. A certified algorithm will return, for each parameter quantities p and p such that piiin S u(f) S pa 'he parameter l controls the running time of the algorithm, and as f we have P p Pt min + 0. Importantly, a priori analysis of the convergence in l is not needed to get certified bounds p(f). By computing prtin and p we have achieved e-approximation of the desired quantity with E = p p
Spin models We focus here on spin models where the state space S is the Cartesian product S = ΣV where E = [-1,,,1]] and V is the vertex set of a graph G = (V, E). The set V is mostly assumed to be finito of size n, but we will highlight some of the results when they can equally be defined when V is infinite. The potential function H is of the form
H(x) 3
Σ he (11) (L))ER
(2)
where he: E² + R is a two-site interaction term which is assumed to be symmetric, i.e. Iz) he ho(z), m.). For example, the ferromagnetic (resp. antiferromagnetic) Ising model corresponds to ((2,,Ty) -zyry (resp. +2zzz). The parameter 8 2 0 plays the role of an inverse-tent teeess-temperaturee Note that when / = 0 then j is the uniform distribution on S, and when 3 = +00, p is supported on the states x E ES that minimize H (i.e., ground states).
Contributions In this paper we study certified algorithms to approximate v(f) when f is local function, i.e. depending only on variables in a set B C V of small size. Note that the antiferromagnetic model encodes the MAX CUT problem for 3 oo and such, computing p(J) in general is hard. This means we have to make assumptions of linear programs giving upper and lower bounds on p(f). The two hicrarchics are based on two different characterizations of the Gibbs measure:
H in order to obtain provable convergence guarantees. We study two hierarchies
The first one is based on imposing the so called DLR Sebbrushin-Lanfrrdd & Velenik (2017) locally on set A S V. The hicrarchy is thus parametrizod by sets satisfying B C A C V. For every such A, we define linear programs giving upper and lower bounds on the quantity of interest:
equations Friedli
Lpmin A,DLR (f) S ((f) S LP
A,DLR (f).
(3)
Note that the interval A.P. (f), A.aaa()))) monotone nonincreasing as a function of and for A V, we have equality in (3). The linear programs parametrized by A have a number of variables and constraints that are of order 20((A)), In Section 2, we prove that under spatial mixing condition the Gibbs distribution μ, have [LPmax A,SLa (f) LPPmi (f)l I [Iflloocie [Ifillsoce==oddis(BBB. "gd for some constants C2 > O (sce Theorem 2.4) and A V \ A. Here dist denotes the distance in the graph G. This implies that to achieve an approximation of s, it suffices to choose the set A so that dist B, A9) is sufficiently large while keeping |A| as small as possible. The spatial mixing condition is satisfied for example for Ising models on a d-dimensional grid provided B is below the critical temperattmp.
2. For the second hicrarchy which is described in Section 3, we consider a Markov chain whose unique stationary distribution is μ. Such Markov chains have been widely studied in the literature on spin systems, cf., Glauber dynamics Martinelli (1999). Tf P is the Markov kernel describing the Markov chain, then stationarity of μ corresponds to the equations p(Pg) = u(g) for all observables g. By imposing this condition on set of local functions g supported on some A V, we arrive at linear
2The technique of coupling from the past Propp Wilson (1996) notable exception where exact convergence can be pertified. is in general more costly that standard MCMC.
2

---

Published in Transactions on Machine Learning Research (11/2025)
programming relaxation the stationary distributions of P. This allows obtain upper and lower bounds on the expectation value p(f) for any local function f supported on B B
LPA, (J) S v(f) S LPP (J)
(4)
The size of the linear programs defining LP (f) and LPX,Mc are also of the order of 2 however they are slightly smaller than the one corresponding the first hierarchy. In fact, as shown in Remark 3.2, under the assumption that P is reversible, this hicrarchy is provably weaker than the first one in the sense that for any A, LP (f) S LP ALR (f) and A.DLR (f) < -P (f). Nevertheless, we show in Theorem 3.3 that when the Markov chain has a mixing time of On n): then we also have [IP A,3c A,Mc (f) X.X & |l -cgdlat(BBBB for some constants 6, 12 > 0.
The convergence result of the first hierarchy relies on the spatial mixing property, property that solely depends on the Gibbs distribution. The convergence result of the second hierarchy however, relies on the fast temporal mixing of the Markov chain P whose invariant distribution is μ.. It turns out that spatial mixing and fast temporal mixing of Glauber Markov chains are equivalent shown in Stroock & Zegarlinski (1992). In fact, our proof of convergence for the second hierarchy uses a result from Dyer et al. (2004) providing a combinatorial argument for the fact that fast temporal mixing implies spatial mixing.
Figure 1: The linear programming hierarchies studied in this paper compute expectation values, with respect to the Gibbs distribution μe, of observables f S + R with small support supp(f) = B. The hierarchies are indexed by subsets 2 B, and convergence expressed in the distance dist(B,A").
Related work The motivation for this paper came from recent work in the theoretical physics literature studying the bootstrap method, and in particular the papers Cho al. (2022); Cho & Sun (2023) where the two hierarchies considered here have been studied numerically, and where asymptotic convergence in the case of infinite systems have been proven. The numerical results of Cho et al. (2022); Cho & Sun (2023) show performance that is comparable to Markov chain Monte Carlo, in particular for systems that are away from criticality. A natural question left open by these works is to prove fast convergence of such certified methods for systems that are far from criticality. Our goal in this paper is to address this question. We note that the use of linear programming to compute moments invariant distributions of general Markov chains goes back to the late 1990s, see e.g., Hernández-Terma Lasserre, 2012, Chapter 12). The problems considered in this paper are also closely related to the literature e.g., Wainwright et al. (2008) and references therein. We note in particular the paper Risteski (2016) which analyzes convex programming hierarchy to estimate the free energy of Ising models in the case of dense graphs. Another recent paper closcly connected to our paper is Bach (2024) which proposes new algorithms based on surs of squares and quantum entropies to estimate log-partition functions. One difference between our paper and the two previously-cited papers is that we consider the problem of computing expectation values of observables directly, rather than the log-partition function. Finally we note that for quantum spin systems, certified algorithms for computing expectation values recently proposed in Fawzi et al. (2024).
Lasserre (1998) and (Hernández-Lemme
&
graphical models and variational inference,
3

---

Published in Transactions on Machine Learning Research (11/2025)
Notations For subset A V, the external boundary of A is defined
8 A = - 3i E A, {i.}} E}
(5)
and we let A = A A A. For spin configuration z E 5 = Σ and i E V we let x' be the configuration where the i'th spin is flipped. Given a function / R,, i V and I E we let
V.f(x) f(x')())))))))))) f(x).
(6)
The support of a function f:S - R is the set of spins that f depends on, i.e.,
supp(f) (iEV: f
(7)
Note that if sfpp(f) S A, then μ(f) = PA(f) where μa is the marginal μ on the sites and we slightly abused notation by using the same notation for f and its restriction to Σ>. Finally, given an arbitrary finitc set X we let. D(X) (p: X + R Exex
1} be the set of probability distributions on X.
2
Linear programming relaxation via the DLR equations
Note that the Gibbs distribution μ defined in (1) satisfies the following identity, for any S and i e V:
p(x) p x
v.H(())
In fact, j is the unique probability distribution satisfying these equalities. Now fix C V and assume that E By locality of the Hamiltonian, the support of V,H included in and so it easy to see that the oquation above is also true for the marginal fa on À A
(x(x) = M V.
Vx E ViEA.
(8)
These are known as the spin-flip equations, the DLR equations. This leads us to define
LP A. {u (D(EA):(( :((x) u(x')e VeH((a) Vo MEEAAiE
By construction, the marginal Mx belongs to A. Given a function f:S - supported on A, we define the following upper and lower bounds on the expectation value u(f):
LP A,DLR (f) = max(v(f) : E LP L
LP (f) = min : E LP p
(10)
Note that LP (f) and LPmin X,DLa (f) are linear programs having |Σ|||| variables and |A | equality constraints.
Relation to local Gibbs states The following proposition gives an alternative characterization of the polyhedron LP L. Before stating the proposition, we need to introduce the following important definition.
Definition 2.1. Let C V and 1 E S We define the local Gibbs state on A with boundary conditions 77 to be
1 27
(x)x) = exp( (x)(-11)(x))
Ya E
where
(](+= = Σ higlri, xz)) +
Σ (LJEE iEN
hig(Fe, 1z)
(KJ)EE
ijeA
is the Hamiltonian truncated to A with boundary conditions 7 and 7 = = exp(-HR(x)))
Proposition 2.2. Tf V E L.P BAR then there a probability distribution p(n) on A such that for any f f +R supported on A we have
x()= Σ (l)
4

---

Published in Transactions on Machine Learning Research (11/2025)
Proof. Let LPA. Fix E0 arbitrary spin configuration on d The constraint in (9) tells us that for any x E ΣA
((a,n) = exp(V,
Vi
This means that the conditional distribution of v on A given that the boundary values are 1 is precisely equal to the local Gibbs distribution un- Hence we let p(n) = Ezega v(z, 7)) be the marginal of v on the boundary, we get the desired equality.
Convergence
The main result in this section shows that if / has spatial mixing, then the linear programming based upper and lower bounds will converge exponentially fast (in dist(supp(f), (5)))))()))))))))))))))))))))))) to u(f). We adopt the following definition of spatial mixing, which corresponds to weak spatial mixing in (Martinelli, 1999, Definitions 2.3).
Definition 2.3. We say that probability measure Y has the spatial miring property if there exist
constants C1 C2 > 0 such that the following holds: for any two subsets B A V and any function f supported on B
sup T
1 - -T(I s C|
E exp(-C2distti,))). iEB.JE0=AA
(11)
For simplicity of expressions, we will be using the following implied inoquality
1x(00) xIOM | A/exp( e/xp((CCdditt nregorrA
sup
(12)
Infinite systems Tf the set of sites V is infinite, e.g-, the Ising model on the grid 2d, defining a Gibbs distribution is more subtle. The most standard definition is via the DLR equations: we say that a probability measure i on the set of configurations {-1,+1} is Gibbs state if for all finite A C V, the marginal of μ A satisfies (8). It turns out that for infinite systems, the Gibbs distribution might not be unique. Note that Definition 2.3 can be used as is for infinite systems as it never refers to Gibbs distributions of the whole system but rather to the distributions 1% which only depend on the Hamiltonian restricted to A. In fact, this condition is known to hold for Ising models on the infinite grid zd provided / is below the corresponding critical temperature (Martinelli, 1999, page 102). Note that when spatial mixing holds, it is simple to see that there is a unique Gibbs distribution that we denote μ.
Theorem 2.4. Assume that the Gibbs distribution u has the spatial miring properly of Desimition 2.3. Then there are constants C1, C2 such that for any B CACV and any f S - R supported on B:
LPmax A,DER (f) - LP A,DUR (f) 1I B| Al c/ C2d
(13)
As a result, if the gruph is d-dimensional grid, p(J) con be approzimated with addilive ervor e by choos- ing A to be a ball of radius ((log(BB\/E) around the sites in B which leads to linear program of size ((o(O()))(BB)))))
Note that this theorem holds even when V is infinite.
Proof. For any L.P p,pa, we know (from Proposition 2.2) that there is a probability distribution p(n) 17 D such that
v(f) = ((())
7
As such for any ʃ supported B C A, and any v LP A,RR, we have
w(f) - sup M.EENOEA
(X()) - xA( < (appl//rrAex(((CAaAAAa (o0(-Cddist(A.) A²),
which proves the desired statement.
5

---

Published in Transactions on Machine Learning Research (11/2025)
3 Linear programming relaxation via a Markov chain
We now consider a second linear programming hierarchy compute expectation values of Gibbs state p.. This hierarchy is based the characterization of the Gibbs distribution μ as the unique stationary measure of some local Markov chain. We start by introducing these local Markov chains which are also known Glauber dynamics. Recall that a Markov chain on S is described via the transition probabilities Pa. of moving from I E to x' E.
Definition 3.1. We say that Markov chain on S with transition matrix P (Pa, )a, is local if it satisfics the following:
1. Single-site update: Pa.a = 0 if x' differs from on more than one site
2. Locality: For any S and E V, P2,x' depends only spins zj for in the neighborhood of
Boundedness: Pa.x' S for all a + ES
We can interpret such a Markov chain as follows: each time step, a site i E is chosen uniformly at random, and the i'th spin is flipped with certain probability e(i, [0, In this case the probability transition matrix is given by: Pa. xi,x)/n, PZæ 1 - Eiev xi,,))nn and Pr,x' 0 otherwise. To satisfy the assumption of locality the function c(t,x) should only depend on the spins 2j where j is in the neighborhood of i.
We will interested local Markov chains whose stationary distribution is the Gibbs distribution p. A sulficient condition for y be stalionary for P is that reversible for μl ie., u(x))a.. = p(x') For local Markov chains, this reversibility condition simplifies to:
Pa Pe.)) = Par ((((")
ViE ZES
which, using the definition of μ is equivalent to
Pa,x* = Pat Pat
(14)
where
v,H())=H(z))-
The simplest concrete example of a Markov chain satislying the above assumptions the heat-bath dynamics, where a given site V is flipped with probability
(pp(-H(2') c(i,z) (x)(-H(x)) + exp( Losp(-H((p H(x)))
e ((H)
(15)
V.H-v
Note that (x) only depends on the spins zj for sites j that are neighbors of i, and so c(i, is local as desired.
Observables Assume (X1)220 is Markov chain with transition matrix P. Il / EV + R is an observable, then Pf(x)= = Ewess Pa (x') the expectation value of f(X1) conditioned on Xo = x. More generally, (P*))(z) the expoctation of f(X,) conditioncd on Xo æ. Using the locality of P, observe that
= Des f())) + Σ d.0 Vif(Ø).
(16)
iEV
iEV
where we used the fact that Pe,x 1 Pa.a. Note that v.f(x) = for i & supp(f). Since Pr.a only depends on the spins zj in the neighborhood of i, the above shows that supp(Pf) S supp(f), where we recall that supp(J) denotes the union of supp(f) with its external boundary.

---

Published in Transactions on Machine Learning Research (11/2025)
LP hierarchy based on Glauber dynamics Assume P is local Markov chain for which is invariant distribution, i.e,, μP μ.. The latter equality can be written equivalently as: p(Pg) p(g) for all observables + R. If we restrict to be supported on a certain A S V, then note that s((g) and ((g) only depend on the marginal of on A. Given a subset A S V we can thus consider the following set of distributions a on A:
MP = {u E D(EA)
{
v(1) =1
(Normalization)
(17)
((h) 2 0 vh 2 0, supp(h) S A (Positivity)
PPg) v(g) 4g supp(g) S (Stationarity) }
By construction, the marginal of μ on A belongs to LPA.CC A.M.. Hence if is any observable supported on B S A we can obtain lower and upper bounds on s(f) by solving the following linear programs:
LPRin AMM ( min(v(() E A.AMC
LP A (f) m max(u(f) : A.A.M.)-
These are linear programs in |S] |Al variables and at most |E||A] equality constraints.
Remark 3.2 (Relation with DLR equations). We show this remark that reversible, then I.P., S LP A.M- Let us rewrite the stationarity condition of (17) more explicitly. supp(g) S A, and a distribution on DA, one can check that
v(Pg) v(g)= = Σ ((P)(()))())) - g(x))
ZEER
Σ v(x)
P PR. Vig(z)
ZE
Σ
Σ
Pa.i Vig(z)
IEA
Σ
Σ
(v(x') -v()) Pa Pa )9(x),
ZEDA EEA
where the last equality we used the fact that EvErs ((>)V = EmerÄ v.f(z)g(z). By taking g to be an indicator function of σ E D we see that the equation v(Pg) = v(g) in the definition of LP A is equivalent to the following linear equation in v:
Σ ZENA YEA
Σ
D((d)) - v(D) = 0.
(18)
AA 7A
The outer summation over æ is effectively a summation over spin values on the boundary A= = A A.
Now using the assumption that P reversible with respect to μ. i.e., that (14) is satisfied, condition (18)
for all a e A simplifies to
Pe. v(o)e (.H()) = Vo EA
(19)
ZEEA SEA
ZA=OA
We see that these equation are implied by the DLR equations (9); more precisely we that the DLR equations impose that each term v(x') =(z)e-vV.ife v.H(2) of the sum is equal to zero. In other words, we have LPAmE LPA.c..
Convergence We are now ready to state our main convergence theorem.
Theorem 3.3. Assume is graph of bounded degree with n = [V] vertices. Let P be a local Markov Choin according to Definition 3.1 and assume that there are constants ci, 02 > 0 such that for all 2 0,
sup |p plll1 Scnn(1 c
(20)

---

Published in Transactions on Machine Learning Research (11/2025)
Let f: .VV - I supported on B C CV. Then there constants C1, C2 > 0 that depend only on A, |B], and C1, such that the following holds: for ony A 2 B
LP
M.M.c( (f) - LP (M) S Cinllf Crill//loae-Coddss..... CCddst((BA")
Remark 3.4. To get accuracy E, the above theorem means that need to take A such that ne- ne-dist(B,A") = O(e), i..., dist(B, A²) ((log((/EE). If the graph is d-dimensional grid, and A is chosen to be ball of radius r ((log(n/e)) centered at B then = 0(r)) and the LP has size 2 So the running time for accuracy is O(exp(Clog' (n/=)))
Remark 3.5. The condition (20) is known to hold for the heat bath dynamics defined in (15) on Ising models at sufficiently small inverse-eemperature B. More specifically, it holds for any Ising model on a graph with degree bounded by A provided A (B) < (Levin & Peres, 2017, Theorem 15.1)
We first need the following important lemma about the speed of propagation of information local Markov chains. We have seen that if f supported on B C V, then f'f is supported on B' := {i e V dist(i, B) s t). The lemma shows that Pt can actually be well approximated by a function supported on Bat/n for constant >> O. This rosult is well-known, see e.g., (Martinelli, 1999, Lemma 3.2). Here, adapt the proof of (Dyer et al., 2004, Lemma 3.1).
Lemma 3.6. Assume the graph G has n = [V] nodes and is of marimum degree A > 1, and let P be a local Markov chain on G. Let f : + R be an observable supported on B, and let A B. Coll dist (iist(,,A A*). Assume z, x' are two spin configurations such that H, z for all A. Then for any integer 2 0,
- Pt P('')) S 0//
where [| and c 1)
Proof. This is a direct consequence of (Dyer et al., 2004, Lemma 3.1) on the speed of propagation of infor- mation of the Glauber dynamics. We reproduce the argument here for convenience. Let (X1)e20 and (Y)) be two copies of the Glauber dynamics with initial configurations respectively x and x' coupled as follows. At each time step t 2 1, the vertex i, chosen uniformly variable. Then X is obtained by flipping the spin in of X, when U S c(be, X2) and otherwise X = X Similarly, Y,+1 is obtained by flipping the spin i, of Y, when U S c(t,, Note that if X, and Y, coincide on the neighborhood of the vertex 4, then Xe+1l&]- = Y++1 [ie], where we use the notation X(i] to denote the i-th spin of the configuration X. More generally for a set B C V,, we let X[B] be the string of spins in the set B. From the interpretation mentioned above, setting Xo = æ and Yo- = x', we have
random in and U [O, 1] is a uniform random
1 -- = (x) f(ri]l
S PB # Y/
S Σ Pr[X, + YD
JEB
Note that in order to have X,(]] # ViDi, there has to be a path of disagreements from A to J. i.e., Up, Ul,..., bt with up E A and vi- = j. Note that for a fixed path Ut, the probability that there exist times 0 < 1 </ ... < 4 S t such that in = for all is at most (9) (1)* The number of simple paths of length I going from A to j is at most A(A - 1)--1 Rocall that dist A B) which implics that I 2 r. Thus,
Pr[MI # YbII < Σ AIA 1)-1() ()
E(A - 1)(t/n)
<
8

---

Published in Transactions on Machine Learning Research (11/2025)
r 2 (A 1)(t/m), bound
(1 ((A-))(/))
Pr[X.[] # KbI] <
In this case, this leads
(e (c(A-1)(t/i) 1)((tm))
(P'f(z) 201
S 8| B exp(-rlogg logr) S |l c/xpppttt///// r).
In the case r (<<(AA-1)(////)
1)(t/n). then we simply the bound
I |lf/l/cceppett/n-
We can now prove Theorem 3.3.
Proof of Theorem 3.3. It is clear that MP (f (f-)(f) = APA (f)- so we may assume in our analysis that u(f) = In this case by the fast mixing assumption, we know that | flloo S |llllooc (cc//1))))))))) c2/n)) for all t 0.
Let f 0 to be fixed later and define g -E P'f. Note that
Pg = g f - pl f.
(21)
Given a function h: 1-1,,+1 (- +R, lot be the function obtained from by fixing all the spins outside A to some fixed value (say all +1) so that supp(h) C By the constraints the linear program, in particular the stationarity equalities in (17), we can then write for any feasible
|n()) = w(f - (PT 1000-99) ğ))l
< || -PP - 5 lxx
(22)
||Pg g g-PP g) pl fllx
< [I(P -g 3o pt1 flloo
where the first inequality follows from the positivity condition of and the following equality follows from (21).
It remains bound each term above individually. The first. term can be bounded using Lemma 3.6. Indeed,
using the latter we have, for any +-1,+13
4-1
1 - [(x) <
Σ t=
(23)
where fist(B,A)) For the first term of (22), note that (P-1)(o-))(a)) = = z] < 2||g øo1oo _ 2t |I ce /i/n-- As result,
|I(- ax S 2t ce
The second term of (22) can be bounded by the mixing time assumption. Putting things together we get
|n()) < 26, fe + [1fllleei(1--op] c2/n)*
Taking t nr/2e gives
I(fi S line Car
for some appropriate constants C1 C2 0, desired.

---

Published in Transactions on Machine Learning Research (11/2025)
4 Outlook
We considered certified algorithms for computing local expectation values of Gibbs distributions using linear programming relaxations. In practice, it is possible to strengthen the hierarchy using additional constraints, for example semidefinit programming (SDP) constraints such as the so-called reflection positivity constraints for the Ising model as done Cho & Sun (2023).
Another variant of the hierarchy we considered in this paper is the following sequence of SDP relaxations (which can be interpreted searching for aum-of-squares certificates) indexed by l N:
={ SDP; = {ue E Dze:
v11) =
(Normalization)
(24)
v v )) 2 vh deg(h) Se
(Positivity)
v ((Pg) v(g) Vg deg(g) S 28 -
A
(Stationarity))
Here, D2t is the set of linear forms the space of degree 26 polynomials in (II,..., kn). Given low-degree
observable f. we let
SDP min(v(f): SDP)),
SDP (ʃ) max(v(f) : v SDP2)
(25)
It is obvious to check that SDPPI (f) and SDPE (f) are lower and upper bounds respectively on μ(f). An open question is to study the convergence rate of SDP (f) SDPP (f) as + oo under appropriate conditions on the mixing time of the Markov chain.
Another interesting direction for future research to use techniques from variational inference improve the performance of the convex relaxations. Variational inference methods Wainwright et al. (2008) are alternatives to MCMC that can be used to compute expectation values of Gibbs distributions (among other things). Such methods are not certified (except for very specific instances), but have very good performance in practice in many settings of interest.
Acknowledgements
HF would like to thank Oisín Fanst and Shengding Sun for discussions that helped simplify the proof of Theorem 3.3, and Minjae Cho for discussions on the bootstrap. HF acknowledges funding from UK Research and Innovation (UKRI) under the UK goverument's 's Horizon Europe funding guarantee EP EP/X032000 /X032051/1, OF acknowledges funding from the European Research Council (ERC Grant AlgoQIP, Agreement No. 851716).
References
Francis Bach. Sum-of-squares relaxations for information theory and variational inference. Foundations of Compulational Mathemalics, pp. 1 1-39, 2024.
Minjae Cho and Xin Sun. Bootstrap, Markov Chain Monte Carlo, and LP/SDP hierarchy for the lattice Tsing model. Journal of High Eneryy Physics, 2023(11):1-22 -2, 2023.
Minjac Cho, Barak Gabai, Ying-Hsuan Hsuan Lin, Victor A Rodriguez, Joshua Sandor, and Xi Yin. Bootstrapping the Ising model the lattice. arXiv preprint Xrxi::220 12538 2022.
Martin Dyer, Alistair Sinclair, Eric Vigoda, and Dror Weitz. Mixing in time and space for lattice spin systems: A combinatorial view. Random Structures & Algorithmss, 2444):461-479 2004.
Hamza Fawzi, Omar Fawzi, and Samuel O Scalet. Certified algorithms for equilibrium states of local quantum hamiltonians. Nature Communications, 15(1):7394, 2024.
Sacha Friedli and Yvan Velenik. Statistioal mechanics of latlice systems: a concrele mathematical indroduc- tion. Cambridge University Press, 2017.
Onésimo Hernández-Lerma and Jean B Lasserre. Linear programming approximations for markov control
processes in metric spaces. Acta Applicandae Mathematica, 51:123-1,: 139, 1998.
10

---

Published in Transactions on Machine Learning Research (11/2025)
Onésimo Hernández-Lerma and Jean B Lasserre. Markow choins and invariant probabilitics, volume 211. Birkhäuser, 2012.
David A Levin and Yuval Peres. Markov chains and miring times, volume 107. American Mathematical Soc., 2017.
David JC MacKay. Information theory, inference and learning algorithms. Cambridge university press, 2003.
Fabio Martinelli.
Lectures glauber dynamics for discrete spin models. Lectures on probability theory and statistics (Saint-Flour, 1997), 1717:93-191, 1999.
James Gary Propp and David Bruce Wilson. Fxact sampling with coupled markov chains and applications statistical mechanics. Random Structures & Algorithms, 9(1-2):223 252, 1996.
Andrej Risteski. How to calculate partition functions using convex programming hierarchies: provable bounds for variational methods. In Conference on Learning Theory, pp. 1402 1416. PMLR, 2016.
Daniel W Stroock and Boguslaw Zegarlinski. The logarithmic sobolev inequality for discrete spin systems on lattice. Communications in Mathematical Physics, 149:175-193, 1992.
Martin J Wainwright, Michael I Jordan, et al. Graphical models, exponential families, and variational inference. Foundations and Trends& in Mochine Learning, 1(1-2):1 (115)):13005, 2008.
11


Reviews:
[{"review_id": "Fndgwk4R1Y", "reviewer": "Reviewer_pptv", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors investigate two hierarchical linear programming (LP) relaxations for computing the expectation of a local function \\$ f \\$ under the Gibbs distribution of a spin system:\n\n1. **DLR-based Hierarchy:**\nThe first approach leverages the Dobrushin-Lanford-Ruelle (DLR) equations, resulting in LPs of size \\$ 2^{\\Theta(|\\Lambda|)} \\$, where \\$ B \\subseteq \\Lambda \\subseteq V \\$. Here, \\$ V \\$ is the set of nodes in the graph, and \\$ B \\$ contains the nodes on which the function \\$ f \\$ depends.\n2. **Glauber Dynamics Hierarchy:**\nThe second approach constructs a hierarchy grounded in Glauber dynamics, again leading to LPs of order \\$ 2^{\\Theta(|\\Lambda|)} \\$.\n\nThe central contribution of the paper is the development of methods that enable the computation of the expectation of local functions, with provable convergence guarantees.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The authors formally propose their methods and provide convergence guarantees in Sections 2 and 3.", "audience_yn": "Yes", "audience_explanation": "I am not an expert in the area and feel that the topic may not be very relevant to a general ML audience. It is mostly relevant to researchers in theoretical physics. However, I do expect a potential overlap with some ML researchers (especially with those who work with Ising models).", "requested_changes": "The paper is generally well written, and the proposed methods are underpinned by solid formal foundations and convergence guarantees. The novelty primarily lies in the introduction of these new LP hierarchies. The proofs draw carefully on results from prior literature.\n\nA notable area for improvement could be the inclusion of empirical comparisons with MCMC-based baseline methods, which would significantly enhance the impact and practical relevance of the study.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "KLpD6mtHtZ", "reviewer": "Reviewer_rRjm", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies two families of linear programming (LP) hierarchies for approximating local expectation values under Gibbs distributions of spin systems. The first hierarchy is based on the DLR (spin-flip) equations, while the second is derived from local Markov chains (such as Glauber dynamics) having the Gibbs state as their stationary distribution. In both constructions, one selects a finite region $\\Lambda \\subseteq V$ (with $V$ the set of sites of the system) that contains the support $B$ of the observable of interest. The LP constraints are then imposed on $\\Lambda$: spin-flip equalities for the DLR hierarchy, or local stationarity conditions for the Markov chain hierarchy. For any such finite $\\Lambda$, the resulting LP provides valid upper and lower bounds on the observable. As $\\Lambda$ grows, these bounds converge to the true expectation value. The authors show that under additional assumptions, this convergence is exponentially fast in the distance between $B$ and the complement $\\Lambda^c$: spatial mixing of the Gibbs distribution suffices for the DLR-based hierarchy (Theorem~2.4), while a temporal mixing bound on the Markov chain (Eq.~(20)) suffices for the MC-based hierarchy (Theorem~3.3). On a $d$-dimensional grid and at sufficiently high temperature, the results imply that $\\varepsilon$-accuracy can be obtained with LPs of quasipolynomial size in $n/\\varepsilon$, specifically $\\exp(O(\\log^d(n/\\varepsilon)))$.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The paper is clearly written and technically correct. The analysis builds on classical results about spatial and temporal mixing of spin systems, combined with standard LP formulations. The main contribution is to present these ingredients in a unified framework and to highlight that LP relaxations provide {certified} bounds on observables. The mathematics is elementary but the exposition is careful and accessible.", "audience_yn": "Yes", "audience_explanation": "The scope of the contribution is somewhat specialized: the focus is on spin systems and Gibbs measures, and the computational guarantees, while quasipolynomial, are unlikely to be practical at scale. Nevertheless, the perspective of using LP hierarchies for certified inference may be of interest to researchers in theoretical machine learning who are concerned with rigorous methods for inference in graphical models. For all these reasons, I am inclined to recommend this paper for publication in TMLR.", "requested_changes": "\\begin{itemize}\n\\item The paper briefly mentions connections to variational inference and to sum-of-squares relaxations (e.g., in the related work section), but this point could be emphasized more strongly. Situating the LP hierarchies between variational inference (widely used but typically uncertified) and the sum-of-squares hierarchy (more powerful but less tractable) would help clarify where this contribution fits within contemporary ML theory.\n\\item How do the certified bounds established in the present paper extend to more general spin systems, such as the Potts model, hard-core models, and proper coloring models? It would be interesting to know if\t the more modern (and powerful) results by Lubetzky and Sly on general spin systems sampling (see, e.g., \\textquotedblleft Cutoff for General Spin Systems with Arbitrary Boundary Conditions\\textquotedblright, CPAM, 2014) could be employed in this direction.\n\\item {p.\\ 1.} When deriving the $\\delta^{-2}$ dependence, it may help to comment explicitly that this comes from a variance bound (e.g., Hoeffding’s inequality applied under approximate independence of samples).\n\\item {p.\\ 2.} The hierarchy based on the Markov chain is claimed to be ``provably weaker.'' If I understand correctly, the justification comes from the discussion around equation (19), and this implication appears to rely on reversibility of the Markov chain. It would be useful to make this step transparent from the very the start, and refer to a later paragraph for the rigorous justification.\n\\item {p.\\ 9.} In equation (22), the first line makes use of the third condition in (17) (the stationarity condition). Pointing this out would improve clarity.\n\\item {Typos.}  \np.~2. ``corresponds'' $\\to$ ``corresponding'';\np.~7. ``imposes'' $\\to$ ``impose''; \np.~8. ``$X[B]$'' should read ``$X[B]$ be''.\n\\end{itemize}", "broader_impact_concerns": "\\begin{itemize}\n\\item The paper briefly mentions connections to variational inference and to sum-of-squares relaxations (e.g., in the related work section), but this point could be emphasized more strongly. Situating the LP hierarchies between variational inference (widely used but typically uncertified) and the sum-of-squares hierarchy (more powerful but less tractable) would help clarify where this contribution fits within contemporary ML theory.\n\\item How do the certified bounds established in the present paper extend to more general spin systems, such as the Potts model, hard-core models, and proper coloring models? It would be interesting to know if\t the more modern (and powerful) results by Lubetzky and Sly on general spin systems sampling (see, e.g., \\textquotedblleft Cutoff for General Spin Systems with Arbitrary Boundary Conditions\\textquotedblright, CPAM, 2014) could be employed in this direction.\n\\item {p.\\ 1.} When deriving the $\\delta^{-2}$ dependence, it may help to comment explicitly that this comes from a variance bound (e.g., Hoeffding’s inequality applied under approximate independence of samples).\n\\item {p.\\ 2.} The hierarchy based on the Markov chain is claimed to be ``provably weaker.'' If I understand correctly, the justification comes from the discussion around equation (19), and this implication appears to rely on reversibility of the Markov chain. It would be useful to make this step transparent from the very the start, and refer to a later paragraph for the rigorous justification.\n\\item {p.\\ 9.} In equation (22), the first line makes use of the third condition in (17) (the stationarity condition). Pointing this out would improve clarity.\n\\item {Typos.}  \np.~2. ``corresponds'' $\\to$ ``corresponding'';\np.~7. ``imposes'' $\\to$ ``impose''; \np.~8. ``$X[B]$'' should read ``$X[B]$ be''.\n\\end{itemize}", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Wecmswhjsn", "reviewer": "Reviewer_6A1y", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work focus on the theoretical convergence of two methods to obtain $\\epsilon$-approximation of the expectation of local function under a Gibbs distribution of a spin system. The two methods are based on linear programs which allow to give lower and upper bounds on the integral. The main theoretical results are upper bounds on the difference of these bounds, which can be decreased to get an $\\epsilon$-approximation.", "claims_and_evidence_yn": "Yes", "claims_explanation": "This paper introduces the necessary background to understand the problem, and all the results are supported by proofs which seem correct.", "audience_yn": "Yes", "audience_explanation": "This paper provides a theoretical result, justifying the use of 2 LP based methods to approximate integrals under a Gibbs distribution of a spin system, which have been used numerically in previous works. I am not very familiar with this literature, but it seems that it could interest some audience (but not sure for the TMLR's audience).", "requested_changes": "The paper is very theoretical, and it is not clear to me what are the motivations to approximate this type of integral. I am aware that there is a literature that studies this problem, but for a TMLR paper, I believe that it would be important to better motivate the problem in the introduction.\n\nIt is also briefly mentioned in the Related work section that the two studied hierarchies have been used numerically in previous works. But it is not clear if these methods work well or not compared to other methods, which would then further motivate their theoretical study.\n\nIt is proved that the DLR bounds are sharper than MC bounds. Does it mean that the DLR method provides better results than the MC one?\n\n\nTypos:\n- Beginning of page 3: in the inequality, $\\mathrm{LP}\\_{\\Lambda,\\mathrm{DLR}}$ instead of $\\mathrm{LP}_{\\Lambda, \\mathrm{MC}}$.\n- Equation (13): misses DLR?", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)