Perturb-and-max-product: Sampling and learning in discrete energy-based models
Authors: Miguel Lazaro-Gredilla
Antoine Dedieu
Dileep George
Venue: NeurIPS 2021 Poster (2021)
Keywords: max-product
perturb-and-map
discrete graphical models
energy based models
belief revision
belief propagation

Abstract: Perturb-and-MAP offers an elegant approach to approximately sample from a energy-based model (EBM) by computing the maximum-a-posteriori (MAP) configuration of a perturbed version of the model. Sampling in turn enables learning. However, this line of research has been hindered by the general intractability of the MAP computation. Very few works venture outside tractable models, and when they do, they use linear programming approaches, which as we will show, have several limitations. In this work we present perturb-and-max-product (PMP), a parallel and scalable mechanism for sampling and learning in discrete EBMs. Models can be arbitrary as long as they are built using tractable factors. We show that (a) for Ising models, PMP is orders of magnitude faster than Gibbs and Gibbs-with-Gradients (GWG) at learning and generating samples of similar or better quality; (b) PMP is able to learn and sample from RBMs; (c) in a large, entangled graphical model in which Gibbs and GWG fail to mix, PMP succeeds.

Perturb-and-max--rroduct Sampling and learning in discrete energy-based models
Miguel Láraro-Gredilla, Antoine Dedieu, Dileep George
Vicarious Al SF Bay Area, CA
{miguel, antoine, d1llee))@ricarioussccmm
Abstract
Perturb-and-MAP offers an elegant approach to approximately sample from an energy-based model (EBM) by computing the maximum-a-posterinn figuration of perturbed version of the model. Sampling in turn enables learning. However, this line of research has been hindered by the general intractability of the MAP computation. Very few works venture outside tractable models, and when they do, they use linear programming approaches, which as we show, have several limitations. In this work. we present perturb-and aaarrr-andmaaxx max-muxx-roduct product (PMP), parallel and scalable mechanism for sampling and learning in discrete EBMs. Models can be arbitrary long as they are built using tractable factors. We show that (a) for Ising models, PMP orders of magnitude faster than Gibbs and Gibbs with- Gradients (GWG) at learning and generating samples of similar or better quality; (b) PMP able to learn and sample from RBMs; (C) in a large, entangled graphical model in which Gibbs and GWG fail to mix, PMP succeeds.
(MAP) con-
Introduction
In this work, concern ourselves with the problem of black-box parameter estimation for very general class of discrete energy-based models (EBMs), possibly with hidden variables. The golden measure of estimation quality is the likelihood of the parameters given the data, according to the model. However, the likelihood of EBMs is not computable in polynomial time due to the partition function, which in turn makes its optimization (parameter estimation) difficult.
This is an active research problem. with multiple approaches having been proposed over the last 14
38 19 These methods typically fall into one or more of the following categories (a) they depart from the maximum likelihood (ML. criterion to obtain a tractable objective, but in doing so they fail when the data does not match the model (e.g., pseudolikelihode [3): (b) they can only be applied to fully observed models (c.g., piecewise training (29]: (c) they
decades [37 29 33 2
13
12
37
17
fail catastrophically for simple models (e.g., Bethe approximation [2]): (d) they require additional designs (such as a custom encoder) for new models, i.e., they are not black-box (e.g. Advil 19)).
A general approach to discrete EBM learning that does not fall in any of those pitfalls is sampling: as we review in Section [2. sampling from an EBM results in efficient learning. Exact sampling is unfeasible in general, but approximate sampling works well enough in multiple models. Gibbs sampling is the most popular because simple, general, and can be applied in black-hox manner: the sampler can be generated directly from the model, even automatically [5] Its main limitation is that for some models, mixing can too slow be usable. Recently the Gibbs- tains-with-rree wrath-grddeetss (GWG) 8 sampler proposed, which significantly improves mixing speed, enabling tractability of more models. completely different approach to learning is Perturb-and-Ma. (22), but this approach has not seen widespread adoption because the general MAP problem could not be solved fast enough or with good enough quality, see Section 2.4.1
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

In this work, we show that max mrx-prrodcct, a relatively under-utilized algorithm, can be used to perform good enough MAP inference many non- non-tractable cases very fast. Thus, not only expands the applicability of perturb-and-MAP to new scenarios (we show that it can sample highly entangled and multimodal EBMs which the state-of-the-art GWG and traditional Perturb-and-MAN fail) but can also be much faster in practice than the state-of-the-art GWG.
2
Background
We review here MI. learning, Perturb-and-MAP, belief propagation, and how they connect.
2.1
Sampling and learning in discrete EBMs
Consider an EBM with variables x described by a set of factors {0a( 0.) = 0 oa (x())) and I unary terms {0,( f.) = 0 o, (,(x)))== Ia is the subset of variables used in factor a, 0a vector of factor parameters and du (Z.: 0a) is a vector of factor sufficient statistics. For the unary terms, the sufficient statistics (x;;;; are the indicator functions, which are set to I when the variable takes given value and to zero otherwise. The energy of the model can be expressed as E(x) E (a ; 0a) Σ (x;;; 02) or, collecting the parameters and sufficient statisties in corresponding vectors, Exx)
Notation:
-_ (x(x) 2,2,(2)) = - -xxxxx).
Given a discrete data set D (2x(m)) 1. where E {0, .O...., M - 1)2, the log- likeliidl of the parameters for a generic EBM in which some variables (Ir)) are hidden can be expressed as
D)
(
Σ log N
1
Σ 2hh
exp(O ( (()()) ,n)) log Σ exp(@Tan (z))) (1)
N
The learning problem can be described efficiently maximizing the key quantity Eq. ( However, Eq. (I) can only be computed for small models, due to the summation over I (the partition function), which involves exponential number of summands. In order to use gradient ascent methods, need the gradient of Eq. ① w.r.f. O, which is VeL(e)D)
N E (D(x(") xx)) p(aal)(m))))) (D(z)))(r=)) From this gradient we can make two standard observations: First, a maximum of Eq. [] the expectation of the sufficient statistics according to the data and the model match, since the gradient needs to be zero. Second, sampling enables efficient learning: we can efficiently approximate the expected sufficient statistics, with no bias, by sampling from the model (both unconditionally and conditioned on (*)). Exact sampling is also unfeasible, so we turn to approximate sampling.
2.2 Perturb-and-MAP
In [22 family of samplers, including the following convenient form, was proposed:
I = argmax (o) (.. + (9, + =()))
(2)
where (()) is random perturbation vector with as many entries as O (i.e., E, |oy]| entries, where Jey] is the cardinality of variable .y)). This means that to generate a sample, first we draw e(() from some distribution, use it to additively perturb the unary terms, and then we solve a standard maximum a posteriori (MAP) problem on the perturbed model.
While MAP inference is NP-complete, it is easier than partition function computation which is #P-complete. In fact, in log-supermodular models, such as attractive Ising models, the exact MAP can be computed using graph-cuts, but the partition function still needs exponential time.
The distribution of the samples (a) will of course depend on the perturbation density. We exclusively use independent, zero-mean densities
Gumbel(e; -c, 1), where x 0.5772 is the Euler-Mascheromi constant. With this choice, models with only unary terms result in the exact Gibbs distribution (*) x exp(OT T & ((x())) (see Lemma 1, (22). More realistic models containing non-unary terms, will in general not match the Gibbs density, but will approximate it.
The Perturb-and-MAP approach also allows for the computation of upper bound the partition function, log Σ, exp(@ θ ((x)) < (max, (0) & (x) + (Ø, + E) &,(x))) p(a) (see Corollary 1,

---

10). If we use this approximation in the log-likelihooe I we get much easier target for optimization
[(O1D) ma 2a ,x,)
xbh + (O, + +)) (x) xnn)
p(c)
max [ La (+: +E) 1((z)) + p(s)
(2))) p(r;))) where p(: O) is a (possibly conditional) sampler of the form of Eq. ( Note that this gradient has the same form as the gradient of the true likelihood, but with the exact sampler replaced by the approximate one. So this can be thought of as the approximate gradient of the true log-likelihood, or as the true gradient of the approximate log- likelihood.
whose gradient Vo((O))) = + Σ (()(n) Ih plenla
The expectations in Eq. ( and its gradient can be easily approximated, without bias, using the sample mean. Thus, can stochastic gradient descent (SGD) and minibatches, to get generally applicable technique for graphical models if we can solve the MAP problem efficiently. So far, the application of perturb-and-MAP has been restricted models whose structure allowed exact MAP solutions or allowed LP relaxations of the MAP problem to produce accurate results. As we show in Section 5.2
even simple Ising models violate these assumptions.
2.3
The Bethe approximation and belief propagation
To make perturb-ande MAP more widely useful, we need a more general (even if approximate) approach to MAP inference. Here we describe a way to achieve this leveraging the Bethe free energy, which uses a variational formulation and modified entropy to approximate the partition function.
Using an arbitrary variational distribution ((x:) we can lower bound the partition function using the expectation of the free energy. At the maximum with respect to this q(z)), obtain the equality:
B log E
exp(30 ((x)) max ((0) T (())))(()) + B --HH(((x)))
(4)
g(x)
where M(((z))) E ((x) log (z) is the entropy and 3 is an inverse temperature parameter. At 1, the l...s. of Eq. becomes the log-partition function. At 3 + do and with perturbing θ the problem posed by Eq. (4) matches the one in Eq. ( We keep this temperature parameter throughout
our derivation. The maximization in 4 still intractable, so we introduce two approximations and we obtain
Σ (de (xa; 0a)) (ala)) + Σ (1: (0: 0:)) + 8 H (xx))
max 9a (o())))n (x.) s.t.
=
ga(a) = q( [4(m)) 11 (x(xa) 2 0.
which is the Bethe approximation of the log-partition function (39 These approximations are (a)
the use of the Bethe free entropy H Hayy (xa), (xx))) Σ. H(qa (x.)) + E1 (x)))(((())) (which is exact for trees, approximate for loopy models) instead of the exact entropy; and (b) the use of locally consistent pseudo-marginal marginals qo (Ia) and (244), which might not correspond to any joint density g(z).
To convert this into unconstrained optimization problem, move to the dual. The dual variables are called "messages", since there one of them associated to each factor-variable pair. See [40] for details. Setting the gradient to O, we recover the belief propagation fixed-point equations
(x:) (1 0,)
(x;)(()))
benb())
ma-s((1)) = 1-111-10 log
(du kemboys (xa)))
where nb(-) denotes the neighbors of factor or variable. Updates must proceed sequentially, factor
by factor. Convergence to fixed point is not guaranteed for loopy models. We recover the primal quantities from the messages. particular, the unary marginals are (xx)) x exp(B(ox 0.) + Etenb() mx))( The recovered primal value in Eq. (5). for 3 1, is approximation of the rog-partiion function, but cannot be used for learning in general, see 23 [2] and Section 5.1
3

---

2.4 Max-product message passing
We are interested in the limit 3 + oo (zero temperature). The lh..s of Eq. (4). i.c. the quantity that we approximating, becomes max, θ ((z), exactly the MAP problem. The update for m re- mains unchanged, the other simplifies to ma-si = max Qa ( 0.) + keebbaa) my (T)))(TT))
The new "unary marginals" x; arg max. { { c; 0.) ) Ebenbca) me- c)} provide the maximizing configuration. This algorithm is known as max cux-prruuct or belief revision. A damping factor 0 S in the message updates can be used without affecting the fixed points but improving convergence. In fact, even parallel updates often provide satisfactory results with as little as a 2'
2.4.1 The connection with linear programming
When taking the limit 3 + oc, the astute reader will have probably noticed that, rather than modifying the fixed point updates, we could have removed the entropy in the r.h.s. of Eq. (5) and obtained a simple linear program (LP). Indeed, this is possible, and there is vast literature
how tackle
the resulting LP (a continuous relaxation of the MAP problem) 27 16 35 20 15 (4]. In fact, most approximate MAP inference research falls in this category, exploiting the convenience of the convexity of LPs.
It is common, although incorrect, to then assume that taking the limit 3 + co in the fixed point equations, i.e. max-pron max-product, is just some sort of coordinate-wir wise dual algorithm to optimize the same LP. This is not the case. The connection between the LP and max- max-produet has been studied in detail in [3.. Max Max-product is a minimax algorithm in the dual, and has different fixed points than the LP.
Here we argue that, for the purpose of sampling, max-product is superior to LP-based solvers. The perturbed MAP problems are highly multimodal, even for simple cases such as Ising models. The Bethe free energy is non-convex function with multiple minima. As 3 + oc, all those minima flatten out, but still exist [36), and become the fixed points of max-product, so the multimodality of the problem is preserved, and hopefully one of the modes is reached. In contrast, the LP approach merges all those minima into a single point or flat region, the global minimum, and the structure is lost. The relaxed solution of the LP on simple Ising models is often 0.5 for all or most variables, result that cannot be rounded to any useful discrete answer. We assess this in our experiments.
3 Perturb-and--mx--podd
(PMP) for sampling and learning
As explained in Section 1.4.1 LP-based solvers are not well-suited for perturb-and-Ma. inference, and fail even in simple cases. Furthermore, max-product is an anytime algorithm whose precision can be controlled by the number of times that the messages are updated, much like the update steps of Gibbs sampling, enabling an easy trade-off between compute time and accuracy.
Perturb-an and aed-max-prre max-product (PMP) learning is presented in Algorithm [ Minibatched PMP samples used to estimate the gradient of the approximate likelihood Eq. (3). The procedure is highly parallel: each time step compute the updates for all the factors (and for all the samples in the minibatch) at the same time. All the computation is local in the graph: updates for messages from a factor or variable only involve messages to that factor variable. In Algorithm[] we use θl to refer to the modificatior[ of the unary terms in θ that clamps the visible variables to the observed values r(()) The damping a is set to }.2 since this value seems to offer good trade-off between accuracy and speed in most cases.
A natural benchmark for this general-purposs algorithm is Gibbs sampling (GS). Both GS and PMP have computational complexity that scales as O(T E Jt), where Si is the number of factors variable is connected to, and T is the number of times we do a full sweep over all the variables. Sequentially, PMP is constant factor slower than GS, given the more involved updates. However, the high parallelism of PMP makes it up to orders of magnitude faster in practice. Recently, a strong competitor GS was presented: Gibbs with gradients (GWG) g. GWG essentially optimizes the
Observve that if visible variables x; we modify the unary parameters 0, so as to attribute score to the visible choice 2, = (() and -oo score to any other choice, the problem of finding the MAP for that new + (()))) is effectively the same as finding the MAP for the original O but clamping the visible variables to (n))

---

order in which variables in the model are updated, thus providing faster mixing for the same number of updates (and with 2-3x times more compute), so we also consider it in our comparisons.
Intuitively, Gibbs sampling propagates less information than PMP, since rather than real-valued messages, it only propagates discrete values. To illustrate this effect, consider simple chain EBM: Gibbs sampling can be exponentially slow [ whereas PMP finds the optimal configuration in linear time. The trade-off is that PMP not unbiased (even asymptotically) and a bit less general than GS: it can only be applied when all the factors in the model are max marginaaanale (the max-product updates for the factor are tractable). Fortunately, most small factors and some not-so- not-so-small one are, and they can be combined without limitation into rich EBMs. Neural functions are typically not easy to max-manee maxammaaiallzze and therefore using a whole neural network (NN) to define the energy not compatible with PMP. A NN that is broken into max- maaginnanal pieces, however, can be used.
One important distinction when using PMP for learning is not to confuse the learned θрмP (which are the parameters that the PMP sampler uses to try match the data distribution) with the parameters of Gibbs distribution. As we show in Section 5.1 the learned θрMP can be very far from the Gibbs distribution parameterization that produces the same density.
Algorithm 1: Learning and sampling with perturb and-max product (PMP)
Input: Sufficient statistic (-)), data (z:(") IR=1, step size y), minibatch size S, sampling steps T Output: A learned set of parameters θ ~ (0, 0.01) for I + to L do for < 1 to S do
// Init
Uniform(1,N)
// Choose sample from the training // Get sample from the posterior // Get a sample from the prior
pmp (6) [m). 2)
pmp (O, T)
end
- θ +
3/2 Een p(vl") - - -
// Stochastic gradient ascent
end
return
Procedure pmp (0, (9, T) ~ Gumbel(-en -, 1) Vi; m 10) (x.)) + i,,i, m (x.) + 0 Va, 1,
// Init
for tort- - to I' do
((-1((x)) maa (m:) +- 6s( Og + +() + Edenbiora m (as) Va, i,
11 Max Max-product
m ( t- 1
(1-1]] (Ix) + maxxmr (Da (D((Ea;0)) 0n) + +k mi ( Va, i,
end
f arg max, { (0)(=== = 0 + + Encu(s) mo-i(zi =) Vi
// Decode
return x
4 Related work
has been proposed in different variants in recent literature. I 25). with the common theme being that either only tractable MAP inference was considered LP solvers (including graph-cuts in this category) were used. As we argue Section
Learning using
22
10 6 26 25
1. this approach does not work on multimodal problems (such as those invariant to permutations
like in Section (.6) and performs poorly even simple settings, as we show in Sections 5.2 and [.3
In the seminal work (22), the authors suggest using QPBO [24] for the MAP problem, which, as convex method, also fails even in small Ising models. The approaches in [0 use the max-product LP (MPLP) 7 to solve the MAP problem, whose poor performance in non-attractive problems we
show in Section 5.2 The recent work (5 is limited to tractable models. The structured learning of
26 also relies on Sug-supermo potentials that make MAP inference tractable.
For instance, n- n-wise factors require (CC") compute (where C' is the cardinality of the variables), logical factors (such as the OR of variables) only need Onn) time [20] cardinality potentials limiting k out of variables to be ON need (kk) 30), other higher-order potentials accept efficient computation [31], etc.
5

---

A different approach is proposed in [6]: it absorbs the parameters of the model in the perturbation distribution, and projects each data point to a concrete perturbation value (the MAP the perturbation given the data point). The parameters of the perturbation are then slightly adjusted towards fitting the projected data, and the process repeated, in the form of hard expectation maximization. The optimizations are "inverse" problems in the continuous space of perturbations, typically solved using quadratic programming. This approach has the advantage of being able to handle any noise distribution, but it is very computationally demanding and does not support Gumbel-distributed perturbations.
Experiments
We now study the empirical performance of PMP for sampling and learning. We measure performance using the maximum mean discrepancy (MMD) 9 between a set of model samples x and set of ground truth samples z':
N
N
MMD MMD3(x, (xx NN
k(x x(0i) + ,((())x(0)) - 2k(x x ((()))
where k(x(9, ,0(0)) exp - (b Ea # 21) is the standard average Hamming distance
kernel. MMD is also used in the original GWG paper 8 Annealed importance sampling, while
popular, is inadequate for this use case, see Appendix B.4 details.
We use the term "full sweep" to refer to a full MCMC update for all variables in GS and GWG, and single update of all the messages for PMP.
Our experiments are run GeForce RTX 2080 Ti. All compared methods are coded on JAX and run with changes only in the samplers for maximum comparability of their runtimes. Code to reproduce these experiments can be found at htttt::/// github. htmp/vicaaroossnnpprttuuuu and_max product/ The specific max-produce updates used in the presented models are derived in Code automating max-product inference in general models can be found in [1].
Appendix C
5.1 Learning and sampling from the "wrong" model
We want to clarify that during learning PMP selects the parameters ÔрмP so to match the expected sufficient statistics To drive this point home, and also to show that PMP can learn and sample from a distribution that cannot be learned with the Bethe approximation, we use the toy model from (23).
data, when sampling. Those are not the parameters of Gibbs distribution.
Consider an EBM with variables, fully connected, with energy function E(x) = -i T,Ij, where x ([1,,+14 This model not learnable using the Bethe approx- imation for > 0.316 [23] We set 0.5 and run PMP learning on infinite data sampled from this model (this is feasible, since the expected statistics can be computed exactly). We use Adam for 200 iterations with 0.01 learning rate: each iteration considers 100 chains and run 100 full sweeps. The result ÔPμP 0.331. One can be tempted to compare this value with the original = 0.5 of the generating Gibbs distribution to assess the quality of the PMP approxi- mation, but that would be incorrect. θPмP is the parameter of the PMP sampler, not of Gibbs distribution. To see this, note that the KI. divergence from data to (a) the Gibbs distribution KL(p(z|0 = 0.5) || p(r(OpиP = 0.331)) = 0.119, is much worse than to (b) the learned sam- pling distribution (x((p==00 = 0.5) PPMP 0PMP = 0.331)) )) & 0.008, which produces an almost perfect match. One way to look this that we are learning the "wrong" model, but it is wrong in such way that compensates for how PMP sampling is also inexact, thus resulting in learned sampling distribution that closer to that of the data.
5.2 MPLPP Max-product is not broken
In the seminal paper "Fixing max-product" 7 Globerson et al. show how the I.P in Eq. (5) (for
+ oo) can be solved via a convergent sequence of message passing updates, which they call max-product LP (MPLP). We show that MPLP can perform very poorly even in simple cases. We

---

first consider a 2D eyclical square lattice with energy E(z)= - - Ei 2, and side sizes {5, 10, 15). We then generate 100 fully-connected Ising model with dimensions {5, 10, 15) and energy [xx) = Li.eeeds WijR,Jj - E, b where Way is uniformly sampled in [-2, and 6; in --0.1,0.1,, with E (-1,11111111 +1}. For both problems, we calculate the exact partition function and estimations using PMP and MPLP (see Section ()) We run PMP for 200 full sweeps, and MPLP until the message updates change by less than 10 (which corresponds to more than 200 full sweeps in practice). The same perturbations are used for both methods, so they both solve the same exact MAP problem. The estimations returned are upper bounds, so the approximation error should be positive. Figure [ we see that this is the case for the 2D lattice, where both approaches perform similarly. However, for Ising models, neither method provides an actual upper bound (since the MAP inference and the expectation are only computed approximately), but MPLP is significantly worse, see Figure [ center]. center). Furthermore, MPLP is sequential algorithm (for which the damping trick used to turn max-product into parallel algorithm does not work practice): for the 15 15 2D lattice sampling from MPLP uses 500x more compute time than PMP, making it impractical for sampling.
Log: partition on eyc. lattice
Log partition lsing model
Ising trained on 'zeros'
[Z (2 sppros) logiZeret)
MPLP (500x slower)
-10
Gibbs 11.0x slowert
PMP
DING
Gibbs-ruset slower?
-20
GWG reset (3 Ox slower)
-40
MPLP (500x slower)
Side of cyclic lattice
Dimensions of Ising model
Full sweeps (sampling steps)
Figure 1: [Left and center] Log-partition estimation error using PMP and MPLP. Closer to zero is better. Shaded area represents 1 std deviation. [Right] See section 5.3 Sampling quality for
different methods, keeping training fixed and varying the number full sweeps at sampling time.
5.3 Ising model: MNIST zeros
In this experiment, we train fully connected Ising model on dataset of contours derived from the Os of the MNIST dataset [18] (see top row of Figure bottoom left)). For learning, we initialize all the parameters to 0 and use Adam for 1000 iterations with 0.001 learning rate. At every iteration, we consider 100 chains and run 50 full sweeps. The MCMC chains are persistent, and non-persistent
version which resets on every learning iteration is marked with the suffix "reset". After training, we produce new samples following each model and method for a varying number of sampling steps (see Figure rigght).
As expected, GWG outperforms Gibbs. The reset versions are truncated, biased samplers, during learning, so they perform better when sampling for a similar number of full sweeps. Despite the prevalence of persistent contrastive divergence (PCD), short-run reset samplers seem best when we want to generate a high quality sample with a truncated MCMC chain. PMP outperforms all the other methods, and does so with much smaller compute time, as indicated in the legend. PMP actually performing around 3x more operations than Gibbs, but its parallelism more than makes up for it.
We also (see Appendix D This approach reaches a log( (MMD") of ---343 -4.33 (0..20). LP performance is poor and the approach is unpractical: took three orders of magnitude more compute. See the LP samples at the bottom row of Figure bottom, left] and how they compare with the other methods.
LP solver (instead of max max- product) with Pertub-and-Man
learn the Ising model
5.4 Structured models: 2D lattices and Erdös-Renyi graphs
Here we reproduce the experiments from 8 with structured models, following precisely the definition, sampling, and learning described. See details in Appendix [ The authors used the -el-menn--uuree error (RMSE) of the estimated parameters (which we include for reference) their performance
7

---

Cyc. cettice
Cyc. 2D lattice
Erdos-Renyl model
-12
PMP
PMP
Gilbs
-15
PMP
GWG-rese (800s)
GWS
GWG
log(AMSE)
-2.0-
-43
-3.0
-444
-2.-
-1.0
-30-
-47
-.5
-3.5
20
40
60
69
100
20
40
100
40
20
100
Learning and sampling, full sweeps
Learning, # MCMC steps equiv
Learning. MCMC steps of equiv
Cyr. lattice
Cyc. lattice
Erdos-Renyl model
-7.5
-7.6
Gibbs (150k)
PMP
-45
PMP
Gibbs
GWG-reset (800x)
GWG
log(MMD*>
-5.5
PMP
-7.9
550
Gibbs
5.5
-7.0
-82
-32
-8.0
20
4.
60
60
100
20
60 60
100
20
40
60 80
1.0
Learning and sampling, # full sweeps
Learning. # MCMC steps equiv
Learning. # MCMC steps equiv
Figure 2: Performance on the lattice model (0 = 0011, side 25) and the Erdös-Renyi graphs from [ [Left column] x-axis determines the number of full sweeps for training and sampling, for all methods. [Center and right columns] x-axis determines the number of individual variable updates for Gibbs and GWG at training time. PMP is trained for a number of full sweeps that result in less training time than GWG at the same x-axis location. To sample, all methods are run for 50 full sweeps.
metric, but this metric is not as meaningful for PMP, since as we saw in Section 5.1 a large difference in the parameters does not imply low quality samples.
Non- non-persistent training: In the left column of Figure 2 we show the performance of training in a 2D lattice from [ with 20 chains and a varying number of full sweeps both for learning and sampling. The chains are reset in between learning iterations. All methods perform well, with PMP taking the lead (even in RMSE), and also being substantially faster (see legend).
Persistent training: This is the actual setting from 8 shown in the remaining plots of Figure 2 Less than a full sweep is run in between iterations and chains are not reset. Samples are produced by fixed 50 full sweeps for all methods. The x-axis measures here the number of individual variable updates in between learning iterations, so all the points in the x-axis correspond to less than a full sweep. A persistent version of PMP is outside the scope of this work (although we provide possible approach in Appendix A). Just as the non con-persistent MCMC samplers, non non-persistent PMP cannot produce meaningful samples in less than full sweep, it cannot compete by that metric. Instead, we adjust the number of full sweeps of PMP to take equal or smaller compute time than GWG at the same x-axis location. This allows a few full sweeps of PMP, sufficient significantly outperform GWG with the same compute budget. Also, notice that once again the persistent regime not the most conducive to the highest sampling quality for any of the methods tested. Hence, despite being faster and more commonly used, this regime might not be the best choice for all applications.
5.5 Restricted Boltzmann machines: learning and sampling
Here we train a restricted Boltzmann machine (RBM) on MNIST data and sample from using PMP. Its energy is defined as F(I., Rx) = =- -E., WayTeyThy - Σ, byTnj E, CaPei with x, (0. 1)N,, I6 E {0, N)., To the best of our knowledge, PMP is the first reported instance of a Berturb-and-MAP method able to train an RBM and sample recognizable MNIST-like digits from it.
Indeed, we ran the author's code from 33 and the samples did not look like MNIST digits by a far margin. Note that no samples are reported in that work. See Appendix B.4 for discussion.
8

---

PCD1 PCD100 G-nosst PMP Data
3
4 6
8
1
0
9 8
-4-
log(
4
6
8 1
S 3
7
2
& / 24 9
3 ( 4 K -; 4 (2 0
3 7
6
1 O
PMP
PCD-100 faster)
Gibbs-reset (3x faster))
PCD-1 (300x faster)
-0-
909900 0- 3 41
101
10a
10a
# Full sweeps (sampling steps)
Figure 3: [Left] Sampling quality for an RBM trained once and sampled for an increasing number of full sweeps. [Right] Samples per method at 100 full sweeps. See Appendix B.2 for more samples.
We train an RBM with 500 hidden units on MNIST using 100 non-persistent full sweeps of PMP and block Gibbs (Gibbs-reset) and also single-sweep (PCD-1) and 100-sweep (PCD-100) persistent block Gibbs variants. We initialize W (O, 0.01) and 6, c NN 1). We use stochastic gradient descent with learning rate 0.01 and train for 200 epochs using minibatches of size 100.
Note that here it does not make sense efficient than a black-box selection of the next variable to sample. The block Gibbs sampler is informed by the structure of the model. However, PMP is not. It is processing all factors in parallel, unaware of the conditional independence in the RBM.
GWG, since the block Gibbs sampler is much more
After learning, the MMD for each model as the number of full sweeps increases is reported in Figure leett] It is interesting to see the marked peak of Gibbs-reset at the 100 full sweeps mark. Indeed, we have learned biased sampler tuned for 100 full sweeps, more MCMC steps make that model worse. Again, PMP outperforms the competition in sample quality. However, in this particular case, the block Gibbs sampler specifically designed for RBMs can leverage as much parallelism as PMP, and having a smaller multiplicative constant due to its simplicity, it is faster than PMP. The trade-off is not being black-box and providing worse MMD figures.
Figure 3 rightt] shows sampling results for each method after 100 full sweeps. Samples for other numbers of full sweeps and a brief discussion can be found in Appendix B.2] Similar to Section 5.3 we tried to use an LP solver instead of max- produut We were able to get it to work to some extent on smaller dataset, as reported in Appendix B.3 but the IP approach was still outperformed by PMP.
5.6
2D blind deconvolution
Finally, we consider the problem of sampling from an EBM with 40680 binary variables that are intricately connected, and which are invariant to permutations.
A binary 2D convolution generative model combines two sets of binary variables W, and S to form a set of binary images X. We illustrate in the top row of Figure 4 for the generation of 2 independent images. Each binary entry of W and S is modeled with an independent Bernoulli prior. The fact that the pixels in W are contiguous and form lines is not captured by the model, and is done only to produce shapes that humans can parse more easily.
W (size: nea x fea x featwian) contains several 2D binary features. S (size: Vimages Tifeat xhx w) is a set of binary indicator variables representing whether each feature is present at each possible image location. The two are combined by convolution, placing the features defined by W at the locations specified by S in order to form an image. The image is binary, so feature overlap acts as an OR, not a sum. The example in the top row of Figure 4 shows how the features W on the left are placed at the locations marked by S, with the top row of S producing the top image in X and the bottom row of S producing the bottom image in X. Each column of S corresponds to the locations of each of the 5 available features in W.
Observe that, unlike S, W is shared by all images. This makes all variables highly interdependent and the posterior highly multimodal: permuting the first dimension of W and the second dimension of S in the same manner does not change X, so this naturally results in multiple equivalent modes. The

---

W
S
X
Data
GWS
G-nesot
pMlppp
Figure 4: Top row: Binary convolution example, with features (W), locations (S) and resulting images (X). Bottom row: [Left] Data and samples from MNIST zeros Ising model after 25 full sweeps of sampling using each tested method. [Center] examples of blind deconvolution dataset (full dataset X in the Appendix B). [Right] samples of the posterior using PMP.
proposed problem Bayesian version of the task in 27, where the weights have been converted into binary variables. Note that Boolean matrix factorization emerges as particular case.
directly turn to do not need to train it. Instead, convolution), we inference: given X (which contains 100 images generated according to a binary 2D sample from the joint posterior p (W, X) Out of those 100 images, are shown in Figure 4 bottom row, center]. On its right, 5 samples (one per row) from the posterior are obtained by PMP after 1000 full sweeps. Each sample comes from consecutive random seed, showing that the model is consistently good at finding samples from the posterior. The dimensions of W used for the generation of X were 4 scenario in which we do not know their ground truth values. The samples from the posterior used that additional feature to include extra copy of the symbols used for generation, but otherwise looked correct. The size of S was set to 100 x 5 x 9 x 9, so that the convolution output matches the size of X. GWG was much slower per full sweep, and were not able to produce any meaningful posterior samples in reasonable time. GWG did work when S or W were clamped to the true values, but not when trying to sample from both simultaneously. All 100 images in X are provided in Appendix B
Since this model has free parameters,
but during inference they were set 5 x 6 to simulate more realistic
6 Discussion
In this work, we present PMP, method for sampling and learning in EBMs. PMP improves over standard, LP-based Perturb-and-Man
that offers MAP solver that fast, parallel, anytime, and produces high quality solutions in problems where I.P-based MAP solvers fail (e.g.. simple non-attractive Ising models and challenging models with symmetries, such as the one in Section (.6).
Compared with Gibbs sampling and GWG, PMP offers faster mixing for a given number of full sweeps, and its parallelism helps complete those full sweeps in substantially less time. PMP also has the ability to solve some problems exponentially faster than Gibbs sampling, both theoretically in the case of chains and experimentally in cases like Section [.6) Its main limitation that is only applicable in EBMs consisting of max-marginn marginaalzzablee factors, although this class is very rich, contains many well-known models, and can be made arbitrarily complex.
A different way to pose this problem is to assume that W a set of parameters of the generative model, and infer them via maximum likelihood. Because we place prior on them and infer their posterior instead, we are doing full Bayesian inference.
10

---

Acknowledgments and Disclosure of Funding
We thank Wolfgang Lehrach and Guangyao Zhou for useful discussions during the preparation of this manuscript.
This work was supported by Vicarious Al and ONR grant N00014-19-1-2366
References
[1] Adrian Barbu and Song-Chun Zhu. Generalizing swendsen- wang to sampling arbitrary posterior
probabilities. IEEE Trunsactions on Pattern Analysis and Machine Intelligence, 27(8): 1239- 1239 1253, 2005.
[2] Julian Besag. Statistical analysis of non-lattice data. Journal of the Royal Statistical Society: Series D (The Statistician), 1793)1179-195, 1975.
[3] Christopher M Bishop. Pattern recognition. Machine learning, 128(9), 2006.
[4] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman- Milne, and Qiao Zhang. JAX: composable transformations of Python+NummPy programs, 2018.
[5] Sarah Depaoli, James P Clifton, and Patrice R Cobb. Just another gibbs sampler (jags) flexible software for meme implementation. Journaa of Educational and Behavioral Statistics, (6): 649, 2016.
[6] Andreea Gane, Tamir Hazan, and Tommi Jaakkola. Learning with maximum a-posteriori perturbation models. In Artificial Intelligence and Statistics, pages 247-256. PMLR, 2014.
[7] Amir Globerson and Tommi Jaakkola. Fixing max-product: Convergent message passing algorithms for map Ip-relaxations. Advances in neural information processing systems, 20:553- 560, 2007.
[8] Will Grathwohl, Kevin Swersky, Milad Hashemi, David Duvenaud, and Chris J Maddison. Oops took a gradient: Scalable sampling for discrete distributions. arXiv preprint arXiv:2102.04509, 2021.
[9] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard SchöIkopf, and Alexander Smola. A kernel two-sample test. The Journal of Machine Learning Research, 13(1):723-773, 2012.
[10] Tamir Hazan and Tommi Jaakkola. On the partition function and random maximum a-posterioni perturbations. arXiv preprint 1rXiv:1206.6410. 2012.
[11] Tamir Hazan, Francesco Orabona, Anand D Sarwate, Subhransu Maji, and Tommi S Jaakkola. High dimensional inference with random maximum a-posteriori perturbations. IEEE Transac- fions on Information Theory, 65(10))6539-6560 2019.
[12] Uri Heinemann and Amir Globerson. What cannot be learned with bethe approximations. arXiv preprint arXiv: 2012.
[13] Aapo Hyvärinen. Consistency of pseudolikelihood estimation of fully visible boltzmann machines. Neural Computation, 18(10):2233-229) 2006.
[14] Aapo Hyvärinen. Some extensions score matching. Computational statistics & data analysis, 51(5):2499-255 2512, 2007.
[15] Vladimir Jojic, Stephen Gould, and Daphne Koller. Accelerated dual decomposition for map inference. In ICML, 2010.
[16] Nikos Komodakis, Nikos Paragios, and Georgios Tziritas. Mrf energy minimization and beyond via dual decomposition. IEEE transactions on pattern analysis and machine intelligence, 3353):531-552, 2010.
11

---

[17] Volodymyr Kuleshov and Stefano Ermon. Neural variational inference and learning in undirected graphical models. In Advances in Neural Information Processing Systems, pages 7734-6743, 2017.
[18] Yann LeCun, Corinna Cortes, and CJ Burges. Mnist handwritten digit database, 2010.
[19] Chongxuan Li, Chao Du, Kun Xu, Max Welling, Jun Zhu, and Bo Zhang. To relieve your headache training an mrf, take advil. arXiv preprint arXiv: 1000::901000000 .84000, 2019.
[20] André FT Martins, Mário AT Figuciredo, Pedro MQ Aguiar, Noah A Smith, and Eric P Xing. Ad3: Alternating directions dual decomposition for map inference in graphical models. The Journal of Machine Learning Research, 16(1):495-545, 2015.
[21] Ofer Meshi, Mehrdad Mahdavi, and Alexander G Schwing. Smooth and strong: Map inference with linear convergence. In NIPS, pages 298-306. Citeseer, 2015.
[22] George Papandreou and Alan L Yuille. Perturb-am aarrb-and-mapp random fields: Using discrete optimiza- tion to learn and sample from energy models. In 2011 International Conference on Computer Vision, pages 193-200. IEEE, 2011.
[23] Zachary Pitkow, Yashar Ahmadian, and Ken Miller. Learning unbelievable probabilities. In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Q Weinberger, editors, Advances in Neural Information Processing Systems, volume 24. Curran Associates, Inc., 2011.
[24] Carsten Rother, Vladimir Kolmogorov, Victor Lempitsky, and Martin Szummer. Optimizing binary mrfs via extended roof duality. 2007 IEEE Conference on Computer Vision and Puttern Recognition, pages 1-8. IEEE, 2007.
[25] Tatiana Shpakova. On Parameter Learning for Perturb-and-MAN Models. PhD thesis, PSL Research University, 2019.
[26] Tatiana Shpakova, Francis Bach, and Anton Osokin. Marginal weighted maximum log- likelihood for efficient learning of perturb-and-map models. arXiv preprint arXiv: 1,0725,, 2018.
[27] Tomáš Singliar and Miloš Hauskrecht. Noisy-or component analysis and its application to link analysis. The Journal of Machine Learning Research, 7:2189-2213 2006.
[28] David Alexander Sontag. Approximate inference graphical models using LP relaxations. PhD thesis, Massachusetts Institute of Technology, 2010.
[29] Charles Sutton and Andrew McCallum. Piecewise training for undirected models. In Proceed- ings of the Twenry-Finn Conference on Uncertainty in Artificial Intelligence, pages AUAI Press, 2005.
[30] Kevin Swersky, Daniel Tarlow, Ilya Sutskever, Ruslan Salakhutdinov, Richard Zemel, and Ryan Prescott Adams. Cardinality restricted boltzmann machines. Advances in neural informa- tion processing systems, 2012.
[31] Daniel Tarlow, Inmar Givoni, and Richard Zemel. Hop-map: Efficient message passing with high order potentials. In Proceedings of the Thirteenth International Conference Artificial Intelligence and Statistics, pages 812-819. JMLR Workshop and Conference Proceedings, 2010.
[32] Jakub M Tomezak, Szymon Zaręba, Siamak Ravanbakhsh, and Russell Greiner. Low- dimensional perturb-and-map approach for learning restricted boltzmann machines. Neural Processing Letters, 50(2):440 2019.
[33] Martin J Wainwright, Tommi Jaakkola, and Alan S Willsky, Tree-reweighted belief propa- gation algorithms and approximate ml estimation by pseudo-moment matching. In AISTATS. volume 3, page 3, 2003.
[34] Martin J Wainwright, Tommi S Jaakkola, and Alan S Willsky. Map estimation via agreement on trees: gressage-passingg and linear programming. IEEE transactions on information theory, 51(11)369773717, 2005.
12

---

[35] Huayan Wang and Koller Daphne. Subproblem-tee tree calibration: A unified approach to max- product message passing, In International Conference on Machine Learning, pages 190-198. PMLR, 2013.
[36] Yair Weiss, Chen Yanover, and Talya Meltzer. Map estimation, linear programming and belief propagation with convex free energies. arXiv preprint 120.52100..5286, 2012.
[37] Max Welling and Charles A Sutton. Learning Markov random fields with contrastive free energies. In A/STATS, 2005.
[38] Sam Wiseman and Yoon Kim. Amortized bethe free energy minimization for learning mrfs, In Advances in Neural Information Processing Systems, pages 15546-15557, 2019.
[39] Jonathan S Yedidia, William T Freeman, and Yair Weiss. Bethe free energy, kikuchi approxima- tions, and belief propagation algorithms. Advances in neural information processing systems, 13, 2001.
[40] Jonathan S Yedidia, William T Freeman, Yair Weiss, al. Generalized belief propagation. In NIPS, volume 13, pages 689-695, 2000.
[41] Guangyao Zhou", Nishanth Kumar", Miguel Lázaro-Gredilla and Dileep George. PGMax: Fac- tor graph discrete variables and hardwaar--aceelerateed differentiable loopy belief propagation in JAX. http:// / htte//giitan github. hot..//////aaaarriooii////caaaa. 2021.
13


Reviews:
[{"review_id": "xudvT7lFbZ", "reviewer": "Reviewer_MkBQ", "summary": "This paper proposes to use max-product belief propagation to provide an approximate solution to the MAP problem encountered in the perturb-and-MAP sampling method. The paper argues that, despite its lack of convergence guarantees, max-product is faster (due in large part to its parallelizability) and better for sampling than linear-programming based methods (the standard approach) in problems with multiple modes. Experiments in a number of synthetic tasks show that the proposed approach, perturb-and-max-product (PMP), is an efficient and effective sampling method that can also be used for parameter estimation of discrete energy-based models in a black-box manner. ", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "# Overall\nIn summary, I found this paper to be a straightforward and pleasant read. The paper very clearly lays out the problem being tackled, the proposed approach, and how this proposed approach fits in to existing work in this area. The proposed method is, to my knowledge, novel. I was somewhat surprised by this, as using max-product for perturb-and-MAP seems like something that would have been tried before; however, I do not know of and did not find existing work taking this approach so it indeed seems novel, and the fact that the method appears obvious in hindsight takes nothing away from the paper. More importantly, the method outperforms the compared baselines quite consistently. While the experimental domains are all synthetic and straightforward, they are common domains for this problem and do demonstrate the benefits of the method, both for sampling and parameter estimation. That said, there are some minor issues with clarity and explanation of one of the experiments (blind deconvolution), and a fairly major issue with a missing baseline in the experiments, making the comparisons seem unfair.\n\n# Quality\nThe submission is sound and claims are mostly well-supported (see below). The paper does a great job of justifying the decisions made throughout and providing supporting evidence for them. The authors are careful and honest about the method and its evaluation in relation to what is contained in the paper. I enjoyed reading the paper until I reached the experimental section where it seems the most natural PMP baseline---that of perturb-and-MAP with a LP solver---is not included (indeed, this is even the standard method reference in the abstract). This seems like a fairly significant issue with the experiments, and I do not see why it would not be included. Even if it takes much longer to run, then run it for a reasonable amount of time and show the results. I am aware that the authors do compare to MPLP in section 5.2 and state that it is too slow for sampling, but it could presumably at least be used for 1-10 sweeps (instead of 200) to get the timescale down to a reasonable amount. \n\nBeyond this, there are two other areas where the authors could further strengthen the paper. The first would be to provide a better understanding and quantification of the relationship between the learned parameters $\\Theta_{PMP}$ and the parameters of the Gibbs distribution (as discussed in Section 5.1), and what role, if any, this plays when using the learned model. The second would be to show the benefit of PMP in a more challenging experiment on real-world data. I do not think that either of these are necessary but do think their inclusion would strengthen the paper.\n\n\n# Clarity\nAs stated above, the paper is overall quite clear, well-organized, and informative to the reader. However, one area in which I found the paper confusing and had to re-read sections to properly understand is in relation to the hidden variables $x_h$. These are introduced in eq (1) clearly enough but then disappear in the approximation, eq (3). When reading through, I interpreted this as implying that the remainder of the paper would not concern itself with hidden variables. But then the update rules in Algorithm 1 and the text on lines 146-147 indicate that the observations $x^{(n)}$ are in fact partial and thus there are hidden variables but the hidden variables are never mentioned explicitly after eq (1). I would encourage the authors to be more consistent and specific in their use of the hidden variables in eq (3) and beyond, and to mention them explicitly in the text beyond section 2.1.\n\nThe second area that I found difficult to understand was the exact setup of the blind deconvolution experiment. Some of my questions about this section are as follows\n1. It's unclear to me what the effect of adding a feature and \"making them all bigger\" (lines 300-302) for W when generating samples is. \n2. Why not sample W and S jointly and produce samples that way? \n3. Also if W is being provided for generation, how is this different than the case with GWG discussed on line 306? \n4. What is a \"reasonable time\"? (line 305) What do samples look like if run for an \"unreasonable\" amount of time? \nOverall, I think the presentation of this experiment could be greatly clarified in the main body and expanded if necessary in the supplement (e.g., even if it's a very simple generative model, please write it down for reference). Finally, this experiment shows that PMP can recover some of the shapes in the dataset but does not show that the samples match the dataset at all (images of single larger shapes vs of multiple smaller overlapping shapes). Is this the goal? Unclear from the text. \n\nOther minor issues with clarity:\n- Line 43 references \"ML learning\" without defining the acronym.\n- Related to the hidden variables $x_h$, the definition of the dataset on line 52 does not admit the inclusion of missing values unless $D < |\\text{domain}(x)|$ unless ${0,1}^D \\subset \\text{domain}(x)$, but then domain($x$) is not defined.\n- The notation $ \\Theta {|{x^{(n)}}} $ to represent clamping of variables is confusing, because $\\Theta$ is the set of parameters, not a function of the variables. If I understand the notation correctly, this should instead be $\\Phi|_{x^{(n)}}$ since $\\Phi$ is a function of $x$ and $\\Theta$ is not.\n- Can you speak to the poor performance of PMP on MMD for the Erdos-Renyi model in Fig 2?\n\n# Significance\nThe method is novel and interesting and I think will prove useful for other researchers but issues with the experiments hold the paper back.\n\n----------------------------------------------------\n### After author response\n\nThank you for your thorough response. I'm satisfied with the response (as much as possible without seeing an updated version of the paper and results) and have increased my score above the acceptance threshold. I think the paper will be much stronger after the information provided in the author response has been included.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes, the discussion in section 3 and elsewhere provides a nice summary of the limitations and tradeoffs of the method. Societal impact is not discussed but the presented method does not require such a discussion.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "vn1f0M-w3Pr", "reviewer": "Reviewer_E7Us", "summary": "The paper presents the perturb-and-max-product (PMP) algorithm for sampling and learning in discrete energy-based models (EMB). The algorithm is scalable and parallel. Empirical evaluation is conducted on Ising models, 2D lattices and Erodos-Renyi graphs, and restricted Boltzmann machines. PMP is compared against Gibbs and Gibbs-with-Gradients under various settings.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I do not have any major comments to make. All-in-all, I think the algorithm is clear, and the paper is well written and executed. The empirical evaluation demonstrates the competitiveness of PMP to the previous mentioned comparable methods.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations are addressed and potential negative societal impact is N/A.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "JxfXjEV4rqI", "reviewer": "Reviewer_hSgj", "summary": "Perturb & Max is an established method for sampling (and therefore learning) a certain class of energy-based models, namely latent-variable models whose joint distribution over visibles & latents is exponential family. However, the max operation poses computational difficulties and existing approaches based on LP solvers are rather limited. This paper proposes an alternative approach, called PMP, that is based on max-product messaging passing. PMP shows favourable empirical performance compared to LP-based solvers, as well as a recent state-of-the-art Gibbs-with-Gradients.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Overall verdict**\n\nThis paper was well-communicated, had a couple of surprising insights and strong empirical results. I enjoyed reading it, and recommend it for acceptance. My only complaint is that I felt that certain limitations of the method (and advantages of competing methods) were understated. \n\n**Originality & Quality**\n\nPROS\n\n1) One of the author's key observations is that the LP & max-product solutions to Equation 4 yield different results. This is an interesting observation and, crucially, enables the authors to improve upon existing LP solvers.\n\n2) The exprimental results are generally strong. In particular, the RBM results are quite remarkable: PMP outperforms block-gibbs (visually & in terms of MMD) *without* knowledge of the model's conditional independencies (although you do pay for it in terms of computational cost).\n\n3) The parallelizability of PMP is very nice, and its good to see that there is a JAX implementation.\n\n4) Whilst experiment 5.6 was not easy to understand without delving into the appendix, I apppreciated that the authors had gone to the effort of constructing a high-dim complex task which appears to stretch the limits of state-of-the-art samplers like GWG. I expect this experimental setup could be useful in future work for testing the efficacy of novel discrete samplers.\n\nCONS\n\n1) The authors claim that the \"[energy-based] models can be arbitrary as long as they are built using tractable factors\", and their method handles \"black-box parameter estimation for a very general class of models\". Whilst neither claim is false, I think they give a misleading impression. I would prefer more direct statements of the form: \"We assume a joint distribution, over observed and latent variables, is exponential family (i.e. the energy is linear in the parameters)\". Whilst such models are very useful and popular, they do have limitations. For instance, neural energy functions are not included, whereas competing methods such as Gibbs & Gibbs-with-Gradients are applicable here - see Section 8 in [1].\n\n2) In fact, there are even greater parameteric restrictions than the previous point suggests. The authors say in line 161 that the factors in the model must be \"max-marginalizable\". I wasn't sure exactly what was meant by this phrase (please clarify) but I inferred it to mean that the message-passing updates are tractable to compute. This does seem like a fairly big limitation.\n\n2) The question of which perturbation distribution to use seems pretty crucial, as it controls the amount of bias in PMP. I would have liked to have seen experiments that varied the choice of perturbation distribution (and the parameters of these distributions) to see how much impact it has.\n\n3) I would have liked to see other metrics beyond MMD. Annealed Importance Sampling can be used to approximate the partition function and therefore the average log-likelihood (and for small RBMs, exact computation is feasible).\n\n4) As explained in Section 5.1, the parameters estimated by PMP will not, in general, match ground truth parameters (even for identifiable models). This implication is that PMP isn't exactly a method for learning energy-based models, but rather a method for learning a sampler with an EBM-based parameterisation. This stands in tension with the claim (in the title & abstract) that PMP is a method for sampling *and* learning of EBMS. That said, the RMSE results in Figure 2 do at least provide some empirical support to the claim that PMP can accurately learn the energy function. I would be even more convinced if there were AIS calculations for the RBM experiments, as this would shed some more light on how well the energy function has been estimated.\n\n**Clarity**\n\nIn section 2.3, the author's take an 'apparent detour'. I don't like this writing style; it requires the readers to hold too much in memory, and then re-read parts of section 2.2 after finishing 2.4 in order to make sense of everything. In particular, it should be possible to introduce Equation 4 in a way that directly relates to Equation 3.\n\nEquation 1 (and its gradient in line 57) include latent variables. But Equation 3 and its gradient in line 80 do not. I think Equation 3 should be altered to include latent variables, since this would help in understanding all the steps in Algorithm 1.\n\n**Significance**\n\nI think this paper has a decent chance of being widely read and cited in the EBM literature. Moreover, many communities outside of ML use energy-based models (e.g. Potts models in proteomics) and the algorithm presented here may be very useful to them.\n\n**References**\n\n[1] Grathwohl, Will, et al. \"Oops I Took A Gradient: Scalable Sampling for Discrete Distributions.\" arXiv preprint arXiv:2102.04509 (2021).\n\n**minor points/questions**\n\n- The expression for the energy in line 50 looks incorrect. Where are the parameters \\theta_i?\n- Could PMP be tractably applied to deep Boltzmann machines?\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I have made recommendations regarding technical limitations above. I am happy with the societal impact section.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4.5"}, "extra_text": {}}, {"review_id": "8AfBXyiMzXn", "reviewer": "Reviewer_aLeX", "summary": "The manuscript is clear, well written, and presents some interesting results in the use \nof the Belief Propagation (BP) algorithm to train discrete energy-based models. \n\nThe authors investigate the use of BP in its zero-temperature/max-product version to \napproximately solve the MAP phase of the Perturb-and-MAP approach, as opposed to costly and\nless effective linear programming relaxations or other heuristics. ", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Relations to works in the area are carefully discussed. The authors provide a good coverage\nof experiments and datasets, and compare to standard Gibbs sampling and to the recently introduced\nGibbs-with-gradients.\nUnfortunately, there is little novelty in the technique proposed, but the surprisingly good performance of BP in the proposed approach is a nice discovery. \n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "\nIn my opinion, the major limitations of the paper are the following:\n\n- The approach is not suitable for energy models that are not highly factorized, \n  as in the notable case of an energy function expressed by a multi-layer perceptron.\n  \n- On loopy graphical models, BP is an uncontrolled approximation. Comfortably though, \n  their experiments on loopy systems show good results.  \n \n- The framework doesn't allow the learning of an explicit energy model, but of an implicit sampler where the epsilon noise distribution plays an important role. Therefore, the possibility of composing energies, one of the few advantages of EBM over other generative models are lost. Also, it cannot be used to solve inverse problems (i.e. estimate the energy parameters).\n\nThe paper marks good progress in the field of discrete energy-based models and in \nthe application of message passing algorithms, I'm not sure it's up to the NeurIPS\ninnovation standards though.\n\nMinor comments: \n\n- Some of the hyperparameters used in the experiments are not entirely clear. E.g. from line 229\n  \"100 chains, 50 full sweeps\" I infer that the perturbation epsilon is sampled 100 times,\n  and each time the algorithm is run for 50 iterations to find the approximate MAP. This could be clarified a bit. Also, I could not find the number of gradient steps or epochs used. \n\n- The definition of MMD could be reported to facilitate the reading.\n\n- Fig. 2: not clear how Zapprox is computed", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors present an interesting framework in which they combine the "perturb-and-MAP" strategy for generating samples from a Gibbs distribution with the perspective that incorrect inference (such as substituting the simple and efficient max-product method) can be compensated by learning the model parameters using a matching approximation, giving a mechanism for (approximately) generating new samples from the training data distribution.

The paper is well-written and clear.  Novelty is somewhat limited, since the work is a straightforward combination of several ideas in the literature.  Experimental validation is reasonable / sufficient, although several reviewers requested some extensions, such as including e.g. off the shelf or limited-budget LP solvers, and expanding on the authors' observations about the performance of such methods.  The authors' response clarified and addressed some of the reviewers' concerns.  Reviewers felt the impact of the work was borderline for NeurIPS, and suggested that some of the authors' claims be toned down or moderated by clear discussion of the method's cons, but were generally in favor of acceptance if possible.

(I also should add that there are *many* method for generating samples from discrete distributions beyond Gibbs-like methods, including tempering methods and "discontinuous" Hamiltonian Monte Carlo for Gibbs distributions, and on the "learning to sample" side, discrete (normalizing) flow models, just to name only a few.)

Decision: Accept (Poster)