{'fraction_non_alphanumeric': 0.05204365868910238, 'fraction_numerical': 0.01193792280143255, 'mean_word_length': 4.281146802761933, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Decomposition Bounds for Marginal MAP

Marginal MAP inference involves making MAP predictions in systems defined with latent variables or missing information. It is significantly more difficult than pure marginalization and MAP tasks, for which a large class of efficient and convergent variational algorithms, such as dual decomposition, exist. In this work, we generalize dual decomposition to a generic power sum inference task, which includes marginal MAP, along with pure marginalization and MAP, as special cases. Our method is based on a block coordinate descent algorithm on a new convex decomposition bound, that is guaranteed to converge monotonically, and can be parallelized efficiently. We demonstrate our approach on marginal MAP queries defined on real-world problems from the UAI approximate inference challenge, showing that our framework is faster and more reliable than previous methods.


Introduction
Probabilistic graphical models such as Bayesian networks and Markov random fields provide a useful framework and powerful tools for machine learning. Given a graphical model, inference refers to answering probabilistic queries about the model. There are three common types of inference tasks. The first are max-inference or maximum a posteriori (MAP) tasks, which aim to find the most probable state of the joint probability; exact and approximate MAP inference is widely used in structured prediction [26]. Sum-inference tasks include calculating marginal probabilities and the normalization constant of the distribution, and play a central role in many learning tasks (e.g., maximum likelihood). Finally, marginal MAP tasks are "mixed" inference problems, which generalize the first two types by marginalizing a subset of variables (e.g., hidden variables) before optimizing over the remainder. 1 These tasks arise in latent variable models [e.g., 29,25] and many decision-making problems [e.g., 13]. All three inference types are generally intractable; as a result, approximate inference, particularly convex relaxations or upper bounding methods, are of great interest.
Decomposition methods provide a useful and computationally efficient class of bounds on inference problems. For example, dual decomposition methods for MAP [e.g., 31] give a class of easy-toevaluate upper bounds which can be directly optimized using coordinate descent [37,6], subgradient updates [14], or other methods [e.g., 22]. It is easy to ensure both convergence, and that the objective is monotonically decreasing (so that more computation always provides a better bound). The resulting bounds can be used either as stand-alone approximation methods [6,14], or as a component of search [11]. In summation problems, a notable decomposition bound is tree-reweighted BP (TRW), which bounds the partition function with a combination of trees [e.g., 34,21,12,3]. These bounds are useful in joint inference and learning (or "inferning") frameworks, allowing learning with approximate inference to be framed as a joint optimization over the model parameters and decomposition bound, often leading to more efficient learning [e.g., 23]. However, far fewer methods have been developed for marginal MAP problems.
In this work, we deveop a decomposition bound that has a number of desirable properties: (1) Generality: our bound is sufficiently general to be applied easily to marginal MAP. (2) Any-time: it yields a bound at any point during the optimization (not just at convergence), so it can be used in an anytime way. (3) Monotonic and convergent: more computational effort gives strictly tighter bounds; note that (2) and (3) are particularly important for high-width approximations, which are expensive to represent and update. (4) Allows optimization over all parameters, including the "weights", or fractional counting numbers, of the approximation; these parameters often have a significant effect on the tightness of the resulting bound. (5) Compact representation: within a given class of bounds, using fewer parameters to express the bound reduces memory and typically speeds up optimization.
We organize the rest of the paper as follows. Section 2 gives some background and notation, followed by connections to related work in Section 3. We derive our decomposed bound in Section 4, and present a (block) coordinate descent algorithm for monotonically tightening it in Section 5. We report experimental results in Section 6 and conclude the paper in Section 7.

Background
Here, we review some background on graphical models and inference tasks. A Markov random field (MRF) on discrete random variables x = [x 1 , . . . , x n ] ∈ X n is a probability distribution, where F is a set of subsets of the variables, each associated with a factor θ α , and Φ(θ) is the log partition function. We associate an undirected graph G = (V, E) with p(x) by mapping each x i to a node i ∈ V , and adding an edge ij ∈ E iff there exists α ∈ F such that {i, j} ⊆ α. We say node i and j are neighbors if ij ∈ E. Then, F is the subset of cliques (fully connected subgraphs) of G.
The use and evaluation of a given MRF often involves different types of inference tasks. Marginalization, or sum-inference tasks perform a sum over the configurations to calculate the log partition function Φ in (1), marginal probabilities, or the probability of some observed evidence. On the other hand, the maximum a posteriori (MAP), or max-inference tasks perform joint maximization to find configurations with the highest probability, that is, Φ 0 (θ) = max x α∈F θ α (x α ).
A generalization of max-and sum-inference is marginal MAP, or mixed-inference, in which we are interested in first marginalizing a subset A of variables (e.g., hidden variables), and then maximizing the remaining variables B (whose values are of direct interest), that is, where A ∪ B = V (all the variables) and A ∩ B = ∅. Obviously, both sum-and max-inference are special cases of marginal MAP when A = V and B = V , respectively.
It will be useful to define an even more general inference task, based on a power sum operator: where f (x i ) is any non-negative function and τ i is a temperature or weight parameter. The power sum reduces to a standard sum when τ i = 1, and approaches max x f (x) when τ i → 0 + , so that we define the power sum with τ i = 0 to equal the max operator.
The power sum is helpful for unifying max-and sum-inference [e.g., 36], as well as marginal MAP [15]. Specifically, we can apply power sums with different weights τ i to each variable x i along a predefined elimination order (e.g., [x 1 , . . . , x n ]), to define the weighted log partition function: where we note that the value of (3) depends on the elimination order unless all the weights are equal. Obviously, (3) includes marginal MAP (2) as a special case by setting weights τ A = 1 and τ B = 0. This representation provides a useful tool for understanding and deriving new algorithms for general inference tasks, especially marginal MAP, for which relatively few efficient algorithms exist.

Related Work
Variational upper bounds on MAP and the partition function, along with algorithms for providing fast, convergent optimization, have been widely studied in the last decade. In MAP, dual decomposition and linear programming methods have become a dominating approach, with numerous optimization techniques [37,6,32,14,38,30,22], and methods to tighten the approximations [33,14].
For summation problems, most upper bounds are derived from the tree-reweighted (TRW) family of convex bounds [34], or more generally conditional entropy decompositions [5]. TRW bounds can be framed as optimizing over a convex combination of tree-structured models, or in a dual representation as a message-passing, TRW belief propagation algorithm. This illustrates a basic tension in the resulting bounds: in its primal form 2 (combination of trees), TRW is inefficient: it maintains a weight and O(|V |) parameters for each tree, and a large number of trees may be required to obtain a tight bound; this uses memory and makes optimization slow. On the other hand, the dual, or free energy, form uses only O(|E|) parameters (the TRW messages) to optimize over the set of all possible spanning trees -but, the resulting optimization is only guaranteed to be a bound at convergence, 3 making it difficult to use in an anytime fashion. Similarly, the gradient of the weights is only correct at convergence, making it difficult to optimize over these parameters; most implementations [e.g., 24] simply adopt fixed weights.
Thus, most algorithms do not satisfy all the desirable properties listed in the introduction. For example, many works have developed convergent message-passing algorithms for convex free energies [e.g., 9,10]. However, by optimizing the dual they do not provide a bound until convergence, and the representation and constraints on the counting numbers do not facilitate optimizing the bound over these parameters. To optimize counting numbers, [8] adopt a more restrictive free energy form requiring positive counting numbers on the entropies; but this cannot represent marginal MAP, whose free energy involves conditional entropies (equivalent to the difference between two entropy terms).
On the other hand, working in the primal domain ensures a bound, but usually at the cost of enumerating a large number of trees. [12] heuristically select a small number of trees to avoid being too inefficient, while [21] focus on trying to speed up the updates on a given collection of trees. Another primal bound is weighted mini-bucket (WMB, [16]), which can represent a large collection of trees compactly and is easily applied to marginal MAP using the weighted log partition function viewpoint [15,18]; however, existing optimization algorithms for WMB are non-monotonic, and often fail to converge, especially on marginal MAP tasks.
While our focus is on variational bounds [16,17], there are many non-variational approaches for marginal MAP as well. [27,40] provide upper bounds on marginal MAP by reordering the order in which variables are eliminated, and using exact inference in the reordered join-tree; however, this is exponential in the size of the (unconstrained) treewidth, and can easily become intractable. [20] give an approximation closely related to mini-bucket [2] to bound the marginal MAP; however, unlike (weighted) mini-bucket, these bounds cannot be improved iteratively. The same is true for the algorithm of [19], which also has a strong dependence on treewidth

Fully Decomposed Upper Bound
In this section, we develop a new general form of upper bound and provide an efficient, monotonically convergent optimization algorithm. Our new bound is based on fully decomposing the graph into disconnected cliques, allowing very efficient local computation, but can still be as tight as WMB or the TRW bound with a large collection of spanning trees once the weights and shifting variables are chosen or optimized properly. Our bound reduces to dual decomposition for MAP inference, but is applicable to more general mixed-inference settings.
Our main result is based on the following generalization of the classical Hölder's inequality [7]: For a given graphical model p(x; θ) in (1) with cliques F = {α} and a set of nonnegative weights τ = {τ i ≥ 0, i ∈ V }, we define a set of "split weights" w α = {w α i ≥ 0, i ∈ α} on each variable-clique pair (i, α), that satisfies α|α i w α i = τ i . Then we have where the left-hand side is the powered-sum along order [x 1 , . . . , x n ] as defined in (3), and the right-hand side is the product of the powered-sums on subvector x α with weights w α along the same elimination order; that is, . , x kc ] should be ranked with increasing index, consisting with the elimination order [x 1 , . . . , x n ] as used in the left-hand side.
Proof details can be found in Section E of the supplement. A key advantage of the bound (4) is that it decomposes the joint power sum on x into a product of independent power sums over smaller cliques x α , which significantly reduces computational complexity and enables parallel computation.

Including Cost-shifting Variables
In order to increase the flexibility of the upper bound, we introduce a set of cost-shifting or reparameterization variables δ = {δ α i (x i ) | ∀(i, α), i ∈ α} on each variable-factor pair (i, α), which can be optimized to provide a much tighter upper bound. Note that Φ τ (θ) can be rewritten as, where N i = {α | α i} is the set of cliques incident to i. Applying inequality (4), we have that where the nodes i ∈ V are also treated as cliques within inequality (4), and a new weight w i is introduced on each variable i; the new weights The bound L(δ, w) is convex w.r.t. the cost-shifting variables δ and weights w, enabling an efficient optimization algorithm that we present in Section 5. As we will discuss in Section 5.1, these shifting variables correspond to Lagrange multipliers that enforce a moment matching condition.

Dual Form and Connection With Existing Bounds
It is straightforward to see that our bound in (5) reduces to dual decomposition [31] when applied on MAP inference with all τ i = 0, and hence w i = w α i = 0. On the other hand, its connection with sum-inference bounds such as WMB and TRW is seen more clearly via a dual representation of (5): The tightest upper bound obtainable by (5), that is, , i ∈ α} is a set of pseudo-marginals (or beliefs) defined on the singleton variables and the cliques, and L is the corresponding local consistency polytope defined by Here, H(·) are their corresponding marginal or conditional entropies, and pa α i is the set of variables in α that rank later than i, that is, for the global elimination order [x 1 , . . . , x n ], pa α i = {j ∈ α | j i}. The proof details can be found in Section F of the supplement. It is useful to compare Theorem 4.2 with other dual representations. As the sum of non-negatively weighted conditional entropies, the bound is clearly convex and within the general class of conditional entropy decompositions (CED) [5], but unlike generic CED it has a simple and efficient primal form (5). 4 Comparing  (5) further splits the graph into small cliques (here, edges), introducing additional cost-shifting variables but allowing for easier, monotonic optimization. (d) Primal TRW splits the graph into many spanning trees, requiring even more cost-shifting variables. Note that all three bounds attain the same tightness after optimization.
to the dual form of WMB in Theorem 4.2 of [16], our bound is as tight as WMB, and hence the class of TRW / CED bounds attainable by WMB [16]. Most duality-based forms [e.g., 9,10] are expressed in terms of joint entropies, θ, b + β c β H(b β ), rather than conditional entropies; while the two can be converted, the resulting counting numbers c β will be differences of weights {w α i }, 5 which obfuscates its convexity, makes it harder to maintain the relative constraints on the counting numbers during optimization, and makes some counting numbers negative (rendering some methods inapplicable [8]). Finally, like most variational bounds in dual form, the RHS of (7) has a inner maximization and hence guaranteed to bound Φ τ (θ) only at its optimum.
In contrast, our Eq. (5) is a primal bound (hence, a bound for any δ). It is similar to the primal form of TRW, except that (1) the individual regions are single cliques, rather than spanning trees of the graph, 6 and (2) the fraction weights w α associated with each region are vectors, rather than a single scalar. The representation's efficiency can be seen with an example in Figure 1, which shows a 3 × 3 grid model and three relaxations that achieve the same bound. Assuming d states per variable and ignoring the equality constraints, our decomposition in Figure 1(c) uses 24d cost-shifting parameters (δ), and 24 weights. WMB (Figure 1(b)) is slightly more efficient, with only 8d parameters for δ and and 8 weights, but its lack of decomposition makes parallel and monotonic updates difficult. On the other hand, the equivalent primal TRW uses 16 spanning trees, shown in Figure 1(d), for 16 · 8 · d 2 parameters, and 16 weights. The increased dimensionality of the optimization slows convergence, and updates are non-local, requiring full message-passing sweeps on the involved trees (although this cost can be amortized in some cases [21]).

Monotonically Tightening the Bound
In this section, we propose a block coordinate descent algorithm (Algorithm 1) to minimize the upper bound L(δ, w) in (5) w.r.t. the shifting variables δ and weights w. Our algorithm has a monotonic convergence property, and allows efficient, distributable local computation due to the full decomposition of our bound. Our framework allows generic powered-sum inference, including max-, sum-, or mixed-inference as special cases by setting different weights.

Moment Matching and Entropy Matching
We start with deriving the gradient of L(δ, w) w.r.t. δ and w. We show that the zero-gradient equation w.r.t. δ has a simple form of moment matching that enforces a consistency between the singleton beliefs with their related clique beliefs, and that of weights w enforces a consistency of marginal and conditional entropies.
can be interpreted as a singleton belief on x i , and µ α (x α ) can be viewed as clique belief on x α , defined with a chain rule (assuming where Z i is the partial powered-sum up to x 1:i on the clique, that is, where the summation order should be consistent with the global elimination order o = [x 1 , . . . , x n ]. (2) The gradients of L(δ, w) w.r.t. the weights {w i , w α i } are marginal and conditional entropies defined on the beliefs {µ i , µ α }, respectively, Therefore, the optimal weights should satisfy the following KKT condition The proof details can be found in Section G of the supplement. The matching condition (8) enforces that µ = {µ i , µ α | ∀(i, α)} belong to the local consistency polytope L as defined in Theorem 4.2; similar moment matching results appear commonly in variational inference algorithms [e.g., 34]. [34] also derive a gradient of the weights, but it is based on the free energy form and is correct only after optimization; our form holds at any point, enabling efficient joint optimization of δ and w.

Block Coordinate Descent
We derive a block coordinate descent method in Algorithm 1 to minimize our bound, in which we sweep through all the nodes i and update each block δ Ni = {δ α i (x i ) | ∀α ∈ N i } and w Ni = {w i , w α i | ∀α ∈ N i } with the neighborhood parameters fixed. Our algorithm applies two update types, depending on whether the variables have zero weight: (1) For nodes with τ i = 0 (corresponding to max nodes i ∈ B in marginal MAP), we derive a closed-form coordinate descent rule for the associated shifting variables δ Ni ; these nodes do not require to optimize w Ni since it is fixed to be zero. (2) For nodes with τ i = 0 (e.g., sum nodes i ∈ A in marginal MAP), we lack a closed form update for δ Ni and w Ni , and optimize by local gradient descent combined with line search.
The lack of a closed form coordinate update for nodes τ i = 0 is mainly because the order of power sums with different weights cannot be exchanged. However, the gradient descent inner loop is still efficient, because each gradient evaluation only involves the local variables in clique α.
Closed-form Update. For any node i with τ i = 0 (i.e., max nodes i ∈ B in marginal MAP), and its associated δ Ni = {δ α i (x i ) | ∀α ∈ N i }, the following update gives a closed form solution for the zero (sub-)gradient equation in (8) (keeping the other {δ α j |j = i, ∀α ∈ N i } fixed): where |N i | is the number of neighborhood cliques, and γ α i (x i ) = log . Note that the update in (11) works regardless of the weights of nodes {τ j | ∀j ∈ α, ∀α ∈ N i } in the neighborhood cliques; when all the neighboring nodes also have zero weight (τ j = 0 for ∀j ∈ α, ∀α ∈ N i ), it is analogous to the "star" update of dual decomposition for MAP [31]. The detailed derivation is shown in Proposition H.1 and H.2 in the supplement.
The update in (11) can be calculated with a cost of only O(|N i | · d |α| ), where d is the number of states of x i , and |α| is the clique size, by computing and saving all the shared {γ α i (x i )} before updating δ Ni . Furthermore, the updates of δ Ni for different nodes i are independent if they are not directly connected by some clique α; this makes it easy to parallelize the coordinate descent process by partitioning the graph into independent sets, and parallelizing the updates within each set. To handle the normalization constraint (6) on w Ni , we use an exponential gradient descent: let ; taking the gradient w.r.t. v i and v α i and transforming back gives the following update where η is the step size and pa α i = {j ∈ α | j i}. In our implementation, we find that a few gradient steps (e.g., 5) with a backtracking line search using the Armijo rule works well in practice. Other more advanced optimization methods, such as L-BFGS and Newton's method are also applicable.

Experiments
In this section, we demonstrate our algorithm on a set of real-world graphical models from recent UAI inference challenges, including two diagnostic Bayesian networks with 203 and 359 variables and max domain sizes 7 and 6, respectively, and several MRFs for pedigree analysis with up to 1289 variables, max domain size of 7 and clique size 5. 7 We construct marginal MAP problems on these models by randomly selecting half of the variables to be max nodes, and the rest as sum nodes.
We implement several algorithms that optimize the same primal marginal MAP bound, including our GDD (Algorithm 1), the WMB algorithm in [16] with ibound = 1, which uses the same cliques and a fixed point heuristic for optimization, and an off-the-shelf L-BFGS implementation that directly optimizes our decomposed bound. For comparison, we also computed several related primal bounds, including standard mini-bucket [2] and elimination reordering [27,40], limited to the same computational limits (ibound = 1). We also tried MAS [20] but found its bounds extremely loose. 8 Decoding (finding a configurationx B ) is more difficult in marginal MAP than in joint MAP. We use the same local decoding procedure that is standard in dual decomposition [31]. However, evaluating the objective Q(x B ) involves a potentially difficult sum over x A , making it hard to score each decoding. For this reason, we evaluate the score of each decoding, but show the most recent decoding rather than the best (as is standard in MAP) to simulate behavior in practice. Figure 3 compare the convergence of the different algorithms, where we define the iteration of each algorithm to correspond to a full sweep over the graph, with the same order of time complexity: one iteration for GDD is defined in Algorithm 1; for WMB is a full forward and backward message pass, as in Algorithm 2 of [16]; and for L-BFGS is a joint quasi-Newton step on all variables. The elimination order that we use is obtained by a weighted-min-fill heuristic [1] constrained to eliminate the sum nodes first.

Figure 2 and
Diagnostic Bayesian Networks. . We randomly select half the nodes as max-nodes in these models. We tune the damping rate of WMB from 0.01 to 0.05.
experimented different damping ratios for WMB, and found that it is slower than GDD even with the best damping ratio found (e.g., in Figure 2(a), WMB works best with damping ratio 0.035 (WMB-0.035), but is still significantly slower than GDD). Our GDD also gives better decoded marginal MAP solution x B (obtained by rounding the singleton beliefs). Both WMB and our GDD provide a much tighter bound than the non-iterative mini-bucket elimination (MBE) [2] or reordered elimination [27,40] methods.
Genetic Pedigree Instances. Figure 3 shows similar results on a set of pedigree instances. Again, GDD outperforms WMB even with the best possible damping, and out-performs the non-iterative bounds after only one iteration (pass through the graph).

Conclusion
In this work, we propose a new class of decomposition bounds for general powered-sum inference, which is capable of representing a large class of primal variational bounds but is much more computationally efficient. Unlike previous primal sum bounds, our bound decomposes into computations on small, local cliques, increasing efficiency and enabling parallel and monotonic optimization. We derive a block coordinate descent algorithm for optimizing our bound over both the cost-shifting parameters (reparameterization) and weights (fractional counting numbers), which generalizes dual decomposition and enjoy similar monotonic convergence property. Taking the advantage of its monotonic convergence, our new algorithm can be widely applied as a building block for improved heuristic construction in search, or more efficient learning algorithms.

Supplement A Experiment on Ising grid
Our GDD directly optimizes a primal bound, and is thus guaranteed to be an upper bound of the partition function even before the algorithm converges, enabling a desirable "any-time" property. In contrast, typical implementations of tree reweighted (TRW) belief propagation optimize the dual free energy function [34], and are not guaranteed to be a bound before convergence. We illustrate this point using an experiment on a toy 5×5 Ising grid, with parameters generated by normal ditribution N (0, 2) and half nodes selected as max-nodes for marginal MAP. Figure 4(a)-(b) shows the TRW free energy objective and GDD, WMB upper bounds across iterations; we can see that TRW does violate the upper bound property before convergence, while GDD and WMB always give valid upper bounds.

B More Results on Diagnostic Bayesian Networks
In addtion to the marginal MAP results on BN-1 and BN-2 in main text, we vary the percentage of max-nodes when generating the marginal MAP problems; the reported results in Figure 5(a)-(b) are the best bound obtained by the different algorithms with the first 20 iterations. In all cases, GDD's results are as good or better than WMB. WMB-0.5 (WMB with damping ratio 0.5) appears to work well on sum-only and max-only (MAP) problems, i.e., when the percentage of max-nodes equals 0% and 100% respectively, but performs very poorly on intermediate settings. The far more heavily damped WMB-0.04 or WMB-0.02 work better on average, but have much slower convergence.

C More Results on Pedigree Linkage Analysis
We test our algorithm on additional 6 models of pedigree linkage analysis from the UAI08 inference challenge. We construct marginal MAP problems by randomly selected 50% of nodes to be max-nodes, and report all the results in Figure 6. We find that our algorithm consistently outperforms WMB with the best possible damping ratio.

D Extensions to Junction Graph
Our bound (5) in the main text uses a standard "factor graph" representation in which the cost-shifts {δ α i } are defined for each variable-factor pair (i, α), and are functions of single variables xi. We can extend our bound to use more general shifting parameters using a junction graph representation; this allows us to exploit higher order clique structures, leading to better performance.
Let (C, S) be a junction graph of p(x; θ) where C = {c | c ⊂ V } is the set of clusters, and S = {s = c k ∩ c l | c k , c l ∈ C} is the set of separators. Assume p(x; θ) can be reparameterized into the form,  Figure 6: Marginal MAP inference on additional pedigree linkage analysis models. We randomly selected 50% of nodes as max-nodes in these models. We tune the damping rate of WMB from 0.01 to 0.06, but we omit WMB-0.06 in the plot if WMB-0.05 is already diverged. and the weighted log partition function is rewritten as Φτ (θ) = log τ x exp c∈C θc(xc) . Similar to the derivation of bound (5) in the main text, we can apply Theorem 4.1, but with a set of more general cost-shifting variables δ c s , defined on each adjacent separator-cluster pair (s, c); this gives the more general upper bound, where we introduce the set of non-negative weights w s = {w s i | i ∈ s} on each separator and w c = {w c i | i ∈ c} on each cluster, which should satisfy s∈N se i w s i + c∈N c i w c i = τi, where N se i = {s | i ∈ s} are all the separators that include node i, and N c i = {c | i ∈ c} are all the clusters that include node i. Obviously, our earlier bound (5) in the main text can be viewed as a special case of (14) with a special junction graph whose separators consist of only single variables, that is, S = V .
A block coordinate descent algorithm similar to Algorithm 1 can be derived to optimize the junction graph bound. In this case, we sweep through all the separators s and perform block coordinate update on all {δ c s |∀c ⊇ s} at each iteration. Similarly to Algorithm 1, we can derive a close form update for separators with all-zero weights (that is, τi = 0, ∀i ∈ s, corresponding to s ⊆ B in marginal MAP), and perform local gradient descent otherwise.

E Proof of Thereom 4.1
Proof. Note the Hölder's inequality is where {fj(x)} are arbitrary positive functions, and {ξj} are non-negative numbers that satisfy j ξj = ξ0. Note we extend the inequality by defining power sum with ξj = 0 to equal the max operator. Our result follows by applying Hölder's inequality on each xi sequentially along the elimination order [x1, x2, · · · , xn].

F.1 Background
The log-partition function Φ(θ) has the following variational (dual) form where M(G) is the marginal polytope [35]. Then, for any scalar ε > 0 (including ε → 0 + ), we have As stated in [15,16], we can further generalize the variational form of above scalar-weighted log partition function to the vector-weighted log partition function (3) in the main text, where H(xi|xi+1:n; b) is the conditional entropy on b(x), and is defined as H(xi|xi+1:n; b) = − x b(x) log(b(xi|xi+1:n)). See more details of its derivation in Theorem 4.1 within [15].
A notable special case of (15) is the dual form of marginal MAP by setting weights τ A = 1 and τ B = 0.

F.2 Proof of Thereom 4.2
We now prove the following dual representation of our bound, where is the local consistency polytope, and pa α i = {j ∈ α|j i}. Thereom 4.2 follows directly from (17).
By applying the dual form of the powered sum (15) on each node and clique respectively, we have where pa α i is the set of variables in α that are summed out later than i, M(Gi) and M(Gα) are the marginal polytopes on singleton node i and clique α respectively, which enforce {bi, bα} to be properly normalized. The above equation can be more compactly rewritten as, where the order of min and max are commuted according to the strong duality (it's convex with θ, and concave with b).
The inner minimization min θ θ, b is a linear program, and it turns out can be solved analytically. To see this, given the relationship between θ and δ, we rewrite the linear program as Then, it is easy to observe that the linear program is either equal to θ, b only if b satisfy the marginalization constraint x α\i bα(xα) = bi(xi) for ∀(i, α), or it will become negative infinity. Considering the outer maximization, we have where L(G) is the local consistency polytope that is obtained by enforcing both M and the marginalization constraint.

F.3 Connection with Existing Free Energy Forms
Most variational forms are expresssed in the following linear combination of local entropies [39, 10], where β refers the region, c β refers the general counting number, b β (x β ) is the local belief.
We can rewrite our dual representations (17) as, where pa α i is the set of variables in α that rank later than i. Without loss of generality, assuming xα = [x1, · · · , xi, xj, · · · xc], i.e. xi and xj are adjacent in the order, we can get One can view (19) in terms of (18), by selecting the region β ∈ {i ∈ V } ∪ {α ∈ F} ∪ {pa α i | ∀(i, α)}; some counting numbers c β will be the differences of weights w α j − w α i .

F.4 Matching Our Bound to WMB
After the weights are optimized, our GDD bound matches to WMB bound with optimum weights. A simple weight initialization method matches our bound to WMB with uniform weights on each mini-bucket, which often gives satisfactory result; a similar procedure can be used to match the bound with more general weights as in Section D. We first set wi = 0 for all nodes i. We then visit the nodes xi along the elimination order o = [x1, x2, · · · , xn], and divide xi's neighborhood cliques Ni = {α|α i} into two groups: (1) the children cliques in which all x α\i have already been eliminated, that is, N ch i = {α | ∀j ∈ α\i, j ≺ i in o}; (2) the other, parent cliques N pa i = {α | ∃j ∈ α\i, j i in o}. We set w α i = 0 for all the children cliques (α ∈ N ch i ), and uniformly split the weights, that is, w α i = τi/|N pa i |, across the parent cliques.

G Proof of Therom 5.1
Proof. For each δ α i (xi), the involved terms in L(δ, w) are L α i (δ) = Φw i (δ) + Φwα (δ), where Our result follows by showing that The gradient of Φw i (δ) is straightforward to calculate, wi log The gradient of Φwα (δ) is more involved; see Proposition I.1 for a detailed derivation.
Given the gradients, the moment matching condition (8) in Therom 5.1 obviously holds. We now prove the entropy matching condition in (10). The constraint optimization is Note, when τi = 0, the optimization is trival, so we simply assume τi > 0. We frame the Lagrangian as Note g ≤ 0 (dual feasibility), otherwise max g,λ K(w, λ, g) will approach infinity. The KKT conditions are stationarity: complementary slackness: giwi = 0, g α i w α i = 0.
We plug it into (23) and (24), and obtain the entropy matching condition (10) in the Therom 5.1.

H Derivations of Closed-form Update
We first derive the closed-form update rule for δ α i (xi) in Proposition H.1. We derive the closed-form update rule for the block δN i = {δ α i (xi) | ∀α ∈ Ni} in Proposition H.2.
Proposition H.1. Given max node i in marginal MAP (i.e., τi = 0 ) and one clique α i (i.e. i ∈ α), keeping all δ fixed except δ α i (xi), there is a closed-form update rule, where x α\i = {xj : j ∈ α, j = i}, w α \i = {w α j : j ∈ α, j = i}, and Ni = {α|α i} is the set of all clique factors in the neighborhood of node i. Futhermore, this update will monotonically decrease the upper bound.
Proof. The terms within the bound L(δ, w) that depend on δ α i (xi) are, The sub-gradient of (26) w.r.t. δ α i (xi) equal to zero if and only if, which is "argmax" matching. One sufficient condition of this matching is, which impllies matching of "pseudo marginals". Then, one can pull δ α i (xi) outside from the operator log To prove monotonicity, we substitute above update equation of δ α i (xi) into (26); then we get, Clearly, (27) ≤ (26) by using the fact that maxx[f (x) + g(x)] ≤ maxx f (x) + maxx g(x).
It is easy to verify AδN i (xi) = γ i (xi). from (30). Since A is invertable, one can solve δN i (xi) = A −1 γ i (xi). Then, one can read out the closed-form update rule (28). The monotonicity holds directly by noticing that the update rule (28) are solutions which jointly satisfy equation (25).