Convergence Rates of Smooth Message Passing with Rounding in Entropy-Regularized MAP Inference

TL;DR

This paper analyzes the convergence rates of smooth message passing algorithms with rounding for entropy-regularized MAP inference in graphical models, providing theoretical guarantees on iteration complexity for recovering the true MAP solution.

Contribution

It offers the first theoretical analysis of convergence rates for entropy-regularized message passing algorithms in MAP inference, including conditions for exact recovery.

Findings

Convergence rates depend on regularization parameters and problem structure.

Under certain conditions, the algorithm guarantees recovery of the true MAP solution.

Provides bounds on the number of iterations needed for $ ext{epsilon}$-optimality.

Abstract

Maximum a posteriori (MAP) inference is a fundamental computational paradigm for statistical inference. In the setting of graphical models, MAP inference entails solving a combinatorial optimization problem to find the most likely configuration of the discrete-valued model. Linear programming (LP) relaxations in the Sherali-Adams hierarchy are widely used to attempt to solve this problem, and smooth message passing algorithms have been proposed to solve regularized versions of these LPs with great success. This paper leverages recent work in entropy-regularized LPs to analyze convergence rates of a class of edge-based smooth message passing algorithms to $\epsilon$-optimality in the relaxation. With an appropriately chosen regularization constant, we present a theoretical guarantee on the number of iterations sufficient to recover the true integral MAP solution when the LP is tight and the solution is unique.

Equations

$${\max }_{x_{\mathcal{V}}\in {\chi }^n}{\sum }_{i\in \mathcal{V}}{\theta }_i(x_i)+{\sum }_{ij\in \mathcal{E}}{\theta }_{ij}(x_i,x_j).$$

$$\displaystyle{\mathcal{M}\stackrel{{\scriptstyle\text{def.}}}{{=}}\left\{\bm{\mu}\ :\ \exists\ \mathbb{P},\begin{array}[]{lr}\mathbb{P}_{X_{i}}(x_{i})=\bm{\mu}_{i}(x_{i}),\ \forall i,x_{i}\\ \mathbb{P}_{X_{i},X_{j}}(x_{i},x_{j})=\bm{\mu}_{ij}(x_{i},x_{j}),\\ \quad\forall ij,x_{i},x_{j}\\ \end{array}\right\}}$$

$$\max ⟨\theta ,\mathbf{\mu }⟩\text{ s.t. }\mathbf{\mu }\in \mathcal{M}$$

$$\displaystyle\mathbb{L}_{2}\stackrel{{\scriptstyle\text{def.}}}{{=}}\left\{\bm{\mu}\geq 0\ :\ \begin{array}[]{lr}\bm{\mu}_{i}\in\Sigma_{d}&\forall i\in\mathcal{V}\\ \bm{\mu}_{ij}\in\mathcal{U}_{d}(\bm{\mu}_{i},\bm{\mu}_{j})&\forall ij\in\mathcal{E}\\ \end{array}\right\}$$

$$\max ⟨\theta ,\mathbf{\mu }⟩\text{ s.t. }\mathbf{\mu }\in {\mathbb{L}}_m,$$

$$\min ⟨C,\mathbf{\mu }⟩-\frac{1}{\eta }H(\mathbf{\mu })\text{ s.t. }\mathbf{\mu }\in {\mathbb{L}}_m,$$

$${\left(\mathrm{round}({\mathbf{\mu }}_{\eta }^{\ast })\right)}_i:=\underset{x\in \chi }{\mathrm{argmax}}({\mathbf{\mu }}_{\eta }^{\ast }{)}_i(x)$$

$$\min {\mathcal{D}}_{\Phi }\left(\mathbf{\mu },\exp (-\eta C)\right)\text{ s.t. }\mathbf{\mu }\in {\mathbb{L}}_m,$$

$${\mathcal{X}}_{ij\to i}=\{\mathbf{\mu }:{\mathbf{\mu }}_{ij}\text{mathbbm}1={\mathbf{\mu }}_i\}$$

$${\mathcal{X}}_{ij,i}=\{\mathbf{\mu }:{\mathbf{\mu }}_i^{\top }\text{mathbbm}1=1,\text{mathbbm}{1}^{\top }{\mathbf{\mu }}_{ij}\text{mathbbm}1=1\}$$

$${\mathcal{X}}_{ij\to j}=\{\mathbf{\mu }:{\mathbf{\mu }}_{ij}^{\top }\text{mathbbm}1={\mathbf{\mu }}_j\}$$

$${\mathcal{X}}_{ij,j}=\{\mathbf{\mu }:{\mathbf{\mu }}_j^{\top }\text{mathbbm}1=1,\text{mathbbm}{1}^{\top }{\mathbf{\mu }}_{ij}\text{mathbbm}1=1\}.$$

$$\displaystyle{\textsc{GreedyEdge}(\bm{\mu})=\operatorname*{arg\,max}_{ij\in\mathcal{E}}\left\{\begin{array}[]{lr}\max\{\|\bm{\mu}_{ij}^{(k)}\mathbbm{1}-\bm{\mu}_{i}^{(k)}\|_{1},\\ \|(\bm{\mu}_{ij}^{(k)})^{\top}\mathbbm{1}-\bm{\mu}_{j}^{(k)}\|_{1}\}\end{array}\right\}}$$

$$\begin{array}{rl}L(\lambda ,\xi )& =-{\sum }_{ij\in \mathcal{E}}{\sum }_{x_i,x_j\in \chi }\exp \left(-\eta C_{ij}(x_i,x_j)-{\lambda }_{ij}(x_i)-{\lambda }_{ji}(x_j)-{\xi }_{ij}\right)\\ & -{\sum }_{i\in \mathcal{V}}{\sum }_{x\in \chi }\exp \left(-\eta C_i(x)-{\xi }_i+{\sum }_{j\in N_r(i)}{\lambda }_{ij}(x)+{\sum }_{j\in N_c(i)}{\lambda }_{ji}(x)\right)\\ & -{\sum }_{ij\in \mathcal{E}}{\xi }_{ij}-{\sum }_{i\in \mathcal{V}}{\xi }_i+{\sum }_{ij\in \mathcal{E}}{\sum }_{x_i,x_j\in \chi }\exp (-\eta C_{ij}(x_i,x_j))+{\sum }_i{\sum }_{x\in \chi }\exp (-\eta C_i(x))\end{array}$$

$$S$$

$$+\sum _{i\in \mathcal{V}}\left[\log \sum _{x\in \chi }{e}^{-\eta C_i(x)}+\sum _{x\in \chi }\frac{\eta }{d}{C}_i(x)\right].$$

$$L({\lambda }^{\ast },{\xi }^{\ast })-L({\lambda }^{(1)},{\xi }^{(1)})\le {\mathcal{S}}_0.$$

$$\frac{{ϵ}^2}{4}\le \frac{1}{4}\Vert {\mathbf{\mu }}_{ij}\text{mathbbm}1-{\mathbf{\mu }}_i{\Vert }_1^2\le 2h^2({\mathbf{\mu }}_{ij}\text{mathbbm}1,{\mathbf{\mu }}_i),$$

$$C_{ij}(x_i,x_j)$$

$${\mathcal{D}}_{\Phi }(\mathbf{\mu },\text{mathbbm}1)$$

$$=-H(\mathbf{\mu })$$

$$\min \eta ⟨C,\mathbf{\mu }⟩+{\mathcal{D}}_{\Phi }\left(\mathbf{\mu },\text{mathbbm}1\right)\text{ s.t. }\mathbf{\mu }\in {\mathbb{L}}_m.$$

$$\nabla \Phi ({\mathbf{\mu }}^{\prime })=\nabla \Phi (\text{mathbbm}1)-\eta C,$$

$$\min {\mathcal{D}}_{\Phi }\left(\mathbf{\mu },\exp (-\eta C)\right)\text{ s.t. }\mathbf{\mu }\in {\mathbb{L}}_m$$

$${\mathcal{X}}_{ij\to i}$$

$$\mathcal{L}({\mathbf{\mu }}^{\prime },\alpha )$$

$$=\sum _{x_i,x_j}{\mathbf{\mu }}_{ij}^{\prime }(x_i,x_j)\left(\log \frac{{\mathbf{\mu }}_{ij}^{\prime }(x_i,x_j)}{{\mathbf{\mu }}_{ij}(x_i,x_j)}-1+\alpha (x_i)\right)+\sum _{x_i}{\mathbf{\mu }}_i^{\prime }(x_i)\left(\log \frac{{\mathbf{\mu }}_i^{\prime }(x_i)}{{\mathbf{\mu }}_i(x_i)}-1-\alpha (x_i)\right).$$

$${\mathbf{\mu }}_{ij}^{\prime }(x_i,x_j)$$

$${\mathbf{\mu }}_i^{\prime }(x_i)$$

$$\mathcal{L}(\alpha )=-\sum _{x_i,x_j}{\mathbf{\mu }}_{ij}(x_i,x_j)e^{-\alpha (x_i)}-\sum _{x_i}{\mathbf{\mu }}_i(x_i)e^{\alpha (x_i)}.$$

$${\alpha }^{\ast }(x_i)=\frac{1}{2}\log \frac{{\sum }_{x_j}{\mathbf{\mu }}_{ij}(x_i,x_j)}{{\mathbf{\mu }}_i(x_i)}.$$

$${\mathcal{X}}_{ij,i}$$

$${\mathbf{\mu }}_{ij}^{\prime }(x_i,x_j)$$

$${\mathbf{\mu }}_i^{\prime }(x_i)$$

Full text

1 INTRODUCTION

Undirected graphical models are a central modeling formalism in machine learning, providing a compact and powerful way to model dependencies between variables. Here we focus on the important class of discrete-valued pairwise models. Inference in discrete-valued graphical models has applications in many areas including computer vision, statistical physics, information theory, and genome research (Antonucci et al., 2014; Wainwright and Jordan, 2008; Mezard and Montanari, 2009).

We focus on the problem of identifying a configuration of all variables that has highest probability, termed maximum a posteriori (MAP) inference. This problem has an extensive literature across multiple communities, where it is described by various names, including energy minimization (Kappes et al., 2013) and constraint satisfaction (Schiex et al., 1995). In the binary case, the MAP problem is sometimes described as quadratic-pseudo Boolean optimization (Hammer et al., 1984) and it is known to be NP-hard to compute exactly (Kolmogorov and Zabin, 2004; Cooper, 1990) or even to approximate (Dagum and Luby, 1993). Consequently, much work has attempted to identify settings where polynomial-time methods are feasible. We call such settings “tractable” and the methods “efficient.” A general framework for obtaining tractable methodology involves “relaxation”—the MAP problem is formulated as an integer linear program (ILP) and is then relaxed to a linear program (LP). If the vertex at which the LP achieves optimality is integral, then it provides an exact solution to the original problem. In this case we say that the LP is tight. If the LP is performed over the convex hull of all integral assignments, otherwise known as the marginal polytope $\mathcal{M}$, then it will always be tight. Inference over the marginal polytope is generally intractable because it requires exponentially many constraints to enforce global consistency.

A popular workaround is to relax the marginal polytope to the local polytope $\mathbb{L}_{2}$ (Wainwright and Jordan, 2008). Instead of enforcing global consistency, the local polytope enforces consistency only over pairs of variables, thus yielding pseudo-marginals which are pairwise consistent but may not correspond to any true global distribution. The number of constraints needed to specify the local polytope is linear in the number of edges. More generally, Sherali and Adams (1990) introduced a series of successively tighter relaxations of the marginal polytope, or convex hull, while retaining control on the number of constraints. However, even with these relaxations, it has been observed that standard LP solvers do not scale well (Yanover et al., 2006), motivating the study of solvers that exploit the structure of the problem, such as message passing algorithms.

Of particular interest to this paper are smooth message passing algorithms, i.e. algorithms derived from regularized versions of the relaxed LP (Meshi et al., 2012; Savchynskyy et al., 2011, 2012; Hazan and Shashua, 2008; Ravikumar et al., 2010). These regularized LPs conduce to efficient optimization in practice and have the special property that their fixed points are unique and optimal; however, this comes at the cost of solving an approximation of the true MAP problem and, without rounding, they do not recover integral solutions in general. Non-asymptotic convergence rates to the optimal regularized function value have been studied (Meshi et al., 2012), but guarantees on the number of iterations sufficient to recover the optimal integral assignment of the true MAP problem have not been considered to our knowledge.

In this work we provide a sharp analysis of the entropy-regularized MAP inference problem with Sherali-Adams relaxations. We first characterize the approximation error of the regularized LP in $l_{1}$ distance, based on new results on entropy-regularized LPs (Weed, 2018). We then analyze an edge-based smooth message passing algorithm, modified from the algorithms described in Werner (2007) and Ravikumar et al. (2010). We prove a $O(1/\epsilon^{2})$ rate of convergence of iterates in $l_{1}$ distance. Combining the approximation error and convergence results, we present a guarantee on the number of iterations sufficient to recover of the true integral MAP assignment using a standard vertex rounding scheme when the LP relaxation is tight and the solution is unique.

2 RELATED WORK

The idea of entropy regularization to aid optimization in inference problems is well studied. It is well known that solving a scaled and entropy-regularized linear program over the marginal polytope yields the scaled Gibbs free energy, intimately related to the log partition function, when the temperature parameter equals one (Wainwright and Jordan, 2008). As the temperature parameter is driven to zero, the calculation of the free energy reduces to the value of the MAP problem. However, this problem is intractable due to the difficulty of both computing the exact entropy and characterizing the marginal polytope (Deza and Laurent, 2009). Therefore, there has been much work in trying to turn this observation into tractable inference algorithms. The standard Bethe approximation instead minimizes an approximation of the true entropy (Bethe, 1935). It was show by Yedidia et al. (2003) that fixed points of the loopy belief propagation correspond to its stationary points, but still the optimization problem resulting from this approximation is non-convex and convergence is not always guaranteed.

To alleviate convergence issues, much work has considered convexifying the free energy problem leading to classes of convergent convex belief propagation often derived directly from convex regularizers (Meshi et al., 2009; Heskes, 2006; Hazan and Shashua, 2008; Johnson and Willsky, 2008; Savchynskyy et al., 2012). For instance, Weiss et al. (2007) proposed a general convexified belief propagation and explored some sufficient conditions that enable heuristically recovering the MAP solution of the LP via a convex sum-product variant. However, the approximation error was still unclear and non-asymptotic convergence rates were not considered. A number of algorithms have also been proposed to directly optimize the unregularized LP relaxation often with only asymptotic convergence guarantees such as block-coordinate methods (Werner, 2007; Globerson and Jaakkola, 2008; Kovalevsky and Koval, 1975; Tourani et al., 2018; Kappes et al., 2013) and tree-reweighted message passing (Wainwright et al., 2005; Kolmogorov, 2006). The relationship between the regularized and unregularized problems can equivalently be viewed as applying a soft-max to the dual objective typically considered in the latter to recover that of the former (Nesterov, 2005; Sontag et al., 2011). Many other convergent methods exist such as augmented Lagrangian (Martins et al., 2011; Meshi and Globerson, 2011), bundle (Kappes et al., 2012), and steepest descent (Schwing et al., 2012, 2014) approaches, but again they are difficult to compare without rates.

Most closely related to our work is recent work in convergence analysis of certain smoothed message passing algorithms that aim to solve the regularized LP objective. Savchynskyy et al. (2011) proposed an accelerated gradient method that achieves $O(1/\epsilon)$ convergence to the optimal regularized dual objective value. Convergence of the primal iterates was only shown asymptotically. Meshi et al. (2012) considered a general dual coordinate minimization algorithm based on the entropy-regularized MAP objective. They proved upper bounds on the rate of convergence to the optimal regularized dual objective value; however, closeness to the true MAP assignment was not formally characterized. Furthermore, convergence in the dual objective value again does not make it easy to determine when the true MAP assignment can be recovered. Meshi et al. (2015) later studied the benefits of adding a quadratic term to the LP objective instead and proved similar guarantees. Ravikumar et al. (2010) also considered entropic and quadratic regularization, using a proximal minimization scheme with inner and outer loops. They additionally provided rounding guarantees to recover true primal solutions. However, as noted by the authors, the inexact calculation of the inner loop prevents a convergence rate analysis once combined with the outer loop. Additionally, rates on the inner loop convergence were not addressed.

The approach of this paper can be understood as the bridging the gap between Meshi et al. (2012) and Ravikumar et al. (2010). Our first contribution is a characterization of the approximation error of the entropy-regularized MAP inference problem. We then study an edge-based message passing algorithm that solves the regularized LP, which is essentially a smoothed max-sum diffusion (Werner, 2007) or the inner loop of the proximal steps of Ravikumar et al. (2010). For our main contribution, we provide non-asymptotic guarantees to the integral MAP assignment for this message passing algorithm when the LP is tight and the solution is unique. To our knowledge, this is the first analysis with rates guaranteeing recovery of the true MAP assignment for smooth methods.

3 BACKGROUND

We denote the $d$-dimensional probability simplex as $\Sigma_{d}\stackrel{{\scriptstyle\text{def.}}}{{=}}\left\{p\in\mathbb{R}_{+}^{d}\>:\>\sum_{i}p_{i}=1\right\}$. The set of joint distributions which give rise to $p,q\in\Sigma_{d}$ is defined as $\mathcal{U}_{d}(p,q)\stackrel{{\scriptstyle\text{def.}}}{{=}}\left\{P\in\mathbb{R}^{d\times d}_{+}\>:\>P\mathbbm{1}=p,\ P^{\top}\mathbbm{1}=q\right\}.$ For any two vectors or matrices $p$ and $q$ having the same number of elements, we use $\langle p,q\rangle$ to denote the dot product, i.e. elementwise multiplication then sum over all elements. We use $\|p\|_{1}$ to denote the sum of absolute values of the elements of $p$. The Bregman divergence between $p,q\in\mathbb{R}_{+}^{d}$ with respect to a strictly convex function $\Phi:\mathbb{R}_{+}^{d}\mapsto\mathbb{R}$ is $\mathcal{D}_{\Phi}(p,q)\stackrel{{\scriptstyle\text{def.}}}{{=}}\Phi(p)-\Phi(q)-\langle\nabla\Phi(q),p-q\rangle.$ We will consider the Bregman divergence with respect to the negative entropy $\Phi(p)=-H(p)\stackrel{{\scriptstyle\text{def.}}}{{=}}\sum_{i}p_{i}(\log p_{i}-1)$, where $p$ need not be a distribution. When $p$ is a distribution, this corresponds to the Kullback-Leibler (KL) divergence. The Bregman projection with respect to $\Phi$ of $q\in\mathbb{R}_{+}^{d}$ onto the set $\mathcal{X}$ is defined as $\mathcal{P}_{\mathcal{X}}\left(q\right)\stackrel{{\scriptstyle\text{def.}}}{{=}}\operatorname*{arg\,min}_{p\in\mathcal{X}}\ \mathcal{D}_{\Phi}(p,q)$. The Hellinger distance between $p,q\in\Sigma_{n}$ is defined as $h(p,q)\stackrel{{\scriptstyle\text{def.}}}{{=}}\frac{1}{\sqrt{2}}\|\sqrt{p}-\sqrt{q}\|_{2}$, where $\|\cdot\|_{2}$ is the $l_{2}$-norm. We denote the square of the Hellinger distance by $h^{2}(p,q)$. We will often deal with marginal vectors which are ordered collections of joint and marginal distributions in the form of matrices and vectors, respectively.

3.1 Pairwise Models

For a set of vertices, $\mathcal{V}=\{1,\ldots,n\}$, and edges $\mathcal{E}$, a pairwise graphical model, $\mathcal{G}\stackrel{{\scriptstyle\text{def.}}}{{=}}\{\mathcal{V},\mathcal{E}\}$, is a Markov random field that represents the joint distribution of variables $X_{\mathcal{V}}\stackrel{{\scriptstyle\text{def.}}}{{=}}\left(X_{i}\right)_{i\in\mathcal{V}}$, taking on values from the set of states $\chi=\{0,\ldots,d-1\}$. We assume that each vertex has at least one edge. For pairwise models, the joint distribution can be written as a function of doubletons and singletons: $\textstyle p_{\theta}(x_{\mathcal{V}})\propto\exp\left(\sum_{i\in\mathcal{V}}\theta_{i}(x_{i})+\sum_{ij\in\mathcal{E}}\theta_{ij}(x_{i},x_{j})\right).$ We wish to find maximum a posteriori (MAP) estimates of this model. That is, we consider the integer program:

[TABLE]

The maximization in (Int) can be written as a linear program by defining a marginal vector $\bm{\mu}$ over variable vertices $\{\bm{\mu}_{i}\}_{i\in\mathcal{V}}$ and variable edges $\{\bm{\mu}_{ij}\}_{ij\in\mathcal{E}}$. The vector $\bm{\mu}_{i}\in\mathbb{R}_{+}^{d}$ represents the marginal distribution probabilities on vertex $i$ while the matrix $\bm{\mu}_{ij}\in\mathbb{R}_{+}^{d\times d}$ represents the joint distribution probabilities shared between vertices $i$ and $j$. We follow the notation of Globerson and Jaakkola (2008) and denote indexing into the vector and matrix variables with parentheses, e.g. $\bm{\mu}_{ij}(x_{i},x_{j})$ for $x_{i},x_{j}\in\chi$. The set of marginal vectors that are valid probability distributions is known as the marginal polytope and is defined as

[TABLE]

We can think of $\mathcal{M}$ as the set of mean parameters of the model for which there exists a globally consistent distribution $\mathbb{P}$. We abuse notation slightly and dually view $\theta$ as a potential “vector.” The edge matrix $\theta_{ij}\in\mathbb{R}^{d\times d}$ is indexed as $\theta_{ij}(x_{i},x_{j})$, indicating the element at the $x_{i}$th row and $x_{j}$th column. The vertex vector $\theta_{i}$ is indexed as $\theta_{i}(x_{i})$, indicating the $x_{i}$th element. The MAP problem in (Int) can be shown to be equivalent to the following LP (Wainwright and Jordan, 2008):

[TABLE]

where $\langle\theta,\bm{\mu}\rangle=\sum_{i\in\mathcal{V}}\sum_{x_{i}}\theta_{i}(x_{i})\bm{\mu}_{i}(x_{i})+\sum_{ij\in\mathcal{E}}\sum_{x_{i},x_{j}}\theta_{ij}(x_{i},x_{j})\bm{\mu}_{ij}(x_{i},x_{j})$.

3.2 Sherali-Adams Relaxations

The number of constraints in $\mathcal{M}$ is unfortunately superpolynomial (Sontag, 2010). This motivates considering relaxations of the marginal polytope to outer polytopes that involve fewer constraints. For example, the local outer polytope is obtained by enforcing consistency only on edges and vertices:

[TABLE]

Relaxations of higher orders have also been studied, in particular by Sherali and Adams (1990) who introduced a hierarchy of polytopes by enforcing consistency on joint distributions of increasing order up to $n$: $\mathbb{L}_{2}\supseteq\mathbb{L}_{3}\supseteq\ldots\supseteq\mathbb{L}_{n}\equiv\mathcal{M}$. The corresponding Sherali-Adams LP relaxation of order $m$ is then

[TABLE]

where $1\leq m\leq n$. Because $\mathbb{L}_{m}$ is an outer polytope of $\mathcal{M}$, we no longer have that the solution to (LP) recovers the true MAP solution of (Int) in general. However if the solution to (LP) is integral, then $x_{i}=\operatorname*{arg\,max}_{x}\bm{\mu}_{i}(x)$ recovers the optimal solution of the true MAP problem. In this case, we say $\mathbb{L}_{m}$ is tight.

4 ENTROPY-REGULARIZED MAP

In this section, we present our first main technical contribution, characterizing the approximation error in the entropy-regularized MAP problem for Sherali-Adams relaxations. In contrast to solving the exact (LP), we aim to solve the entropy-regularized LP:

[TABLE]

where $C\stackrel{{\scriptstyle\text{def.}}}{{=}}-\theta$ and $H(\bm{\mu})=\langle\bm{\mu},-\log\bm{\mu}+\mathbbm{1}\rangle$. The hyperparameter $\eta$ adjusts the level of regularization. Denote by $\bm{\mu}^{*}_{\eta}$ the solution of (Reg) where we omit the reference to $m$ to alleviate notation. In addition to their extensive history in inference problems, entropy-regularized LPs have arisen in a number of other fields to aid optimization when standard LP solvers are insufficient. For example, recent work in optimal transport has relied on entropy regularization to derive alternating projection algorithms (Cuturi, 2013; Benamou et al., 2015) which admit almost linear time convergence guarantees in the size of the cost matrix (Altschuler et al., 2017). Some of our theoretical results draw inspiration from these works.

4.1 Approximation Error

When $\mathbb{L}_{m}$ is tight and the solution is unique, we show that approximate solutions from solving (Reg) are not necessarily detrimental because we can apply standard vertex rounding schemes to yield consistent integral solutions. It was shown by Cominetti and San Martín (1994), and later refined by Weed (2018), that the approximation error of general entropy-regularized linear programs converges to zero at an exponential rate in $\eta$. Furthermore, it is possible to determine how large $\eta$ should be chosen in order for rounding to exactly recover the optimal solution to (Int). The result is summarized in the following extension of Theorem 1 of Weed (2018)111The entropy is defined without the linear offset in Weed (2018)..

Theorem 1**.**

Let $\mathcal{R}_{1}=\max_{\bm{\mu}\in\mathbb{L}_{m}}\|\bm{\mu}\|_{1}$, $\mathcal{R}_{H}=\max_{\bm{\mu},\bm{\mu}^{\prime}\in\mathbb{L}_{m}}H(\bm{\mu})-H(\bm{\mu}^{\prime})$, $\mathcal{V}_{m}$ be the set of vertices of $\mathbb{L}_{m}$, and $\mathcal{V}^{*}_{m}\subseteq\mathcal{V}_{m}$ the set of optimal vertices with respect to $C$. Let $\Delta=\min_{V_{1}\in\mathcal{V}_{m}\backslash\mathcal{V}^{*}_{m},V_{2}\in\mathcal{V}_{m}^{*}}\langle C,V_{1}\rangle-\langle C,V_{2}\rangle$ be the smallest gap in objective value between an optimal vertex and any suboptimal vertex of $\mathbb{L}_{m}$. Suppose $\mathbb{L}_{m}$ is tight and $|\mathcal{V}^{*}_{m}|=1$. If $\eta\geq\frac{2\mathcal{R}_{1}\log{64\mathcal{R}_{1}}+2\mathcal{R}_{1}+2\mathcal{R}_{H}}{\Delta}$, the following rounded solution is a MAP assignment:

[TABLE]

Proof.

Define $\widetilde{C}=C+\mathbf{1}\frac{1}{\eta}$, where $\mathbf{1}$ denotes an all-ones vector with the same dimensions as $C$. If $\eta\geq\frac{4\mathcal{R}_{1}}{\Delta}$ then $\widetilde{\mathcal{V}}_{m}^{*}$, the set of optimal vertices of $\mathbb{L}_{m}$ with respect to $\widetilde{C}$, satisfies $\widetilde{\mathcal{V}}_{m}^{*}=\mathcal{V}_{m}^{*}$ and $\min_{V_{1}\in\mathcal{V}_{m}\backslash\widetilde{\mathcal{V}}^{*}_{m},V_{2}\in\widetilde{\mathcal{V}}_{m}^{*}}\langle C,V_{1}\rangle-\langle C,V_{2}\rangle\geq\frac{\Delta}{2}$. If $V\in\widetilde{\mathcal{V}}_{m}^{*}$; and $V^{\prime}\in\mathcal{V}_{m}\backslash\widetilde{\mathcal{V}}_{m}^{*}$, then $\langle\widetilde{C},V^{\prime}\rangle-\langle\widetilde{C},V\rangle\geq\Delta-\frac{1}{\eta}\|V^{\prime}-V\|_{1}\geq\frac{\Delta}{2}$. Let $\widetilde{\Delta}=\frac{\Delta}{2}$. If $\eta\geq\frac{\mathcal{R}_{1}\log{64\mathcal{R}_{1}}+\mathcal{R}_{1}+\mathcal{R}_{H}}{\widetilde{\Delta}}$, and $|\widetilde{\mathcal{V}}_{m}^{*}|=1$ then $2\mathcal{R}_{1}\exp\left(-\eta\frac{\widetilde{\Delta}}{\mathcal{R}_{1}}+\frac{\mathcal{R}_{1}+\mathcal{R}_{H}}{\mathcal{R}_{1}}\right)\leq\frac{1}{32}$. And therefore, by Corollary 9 of Weed (2018) $\min_{\bm{\mu}\in\mathcal{V}^{*}_{m}}\|\bm{\mu}-\bm{\mu}_{\eta}^{*}\|_{1}\leq\frac{1}{32}$. Since $\mathbb{L}_{m}$ is assumed to be tight and $\widetilde{\mathcal{V}}_{m}^{*}=\mathcal{V}_{m}^{*}$ contains a single integral vertex $\bm{\mu}^{*}$, the last equation implies $\mathrm{round}(\bm{\mu}_{\eta}^{*})=\bm{\mu}^{*}$. ∎

Consequently, since $\mathcal{R}_{1}\leq\sum_{j=1}^{m}\binom{n}{j}d^{j}$ and $\mathcal{R}_{H}\leq\sum_{j=1}^{m}\binom{n}{j}\log(d^{j})$222For $m=2$ we can get tighter bounds corresponding to the number of edges in the graph $\mathcal{G}$., we have:

Corollary 1**.**

If $\mathbb{L}_{m}$ is tight, $|\mathcal{V}_{m}^{*}|=1$, and $\eta\geq\frac{\log(8mn^{m}d^{m})+2mn^{m}d^{m}}{\Delta}$, the rounded solution $\mathrm{round}(\bm{\mu}_{\eta}^{*})$ is a MAP assignment.

In general the dependence of $\Delta$ on $\eta$ suggested by Theorem 1 is not improvable (Weed, 2018). Nevertheless, when $m=2$ and $d=2$, since all vertices in $\mathcal{V}_{2}$ have entries equal to either $0,\frac{1}{2}$ or $1$—see Padberg (1989) or Theorem 3 of Weller et al. (2016)—if the entries of $C$ are all integral, we have $\Delta\geq\frac{1}{2}$, thus yielding a more concrete guarantee. The disadvantage of choosing exorbitantly large $\eta$ is that efficient computation of solutions often becomes more difficult in practice (Weed, 2018; Benamou et al., 2015; Altschuler et al., 2017). Thus, in practice, there exists a trade-off between computation time and approximation error that is controlled by $\eta$. We will provide a precise theoretical characterization of the trade-off in Section 6. In our guarantees, multiplying $C$ by a constant $a$ (and therefore multiplying $\Delta$ by $a$) is equivalent to multiplying $\eta$ by the same value.

4.2 Equivalent Bregman Projection

The objective (Reg) can be interpreted as a Bregman projection. This interpretation has been explored by Ravikumar et al. (2010) as a basis for proximal updates and also Benamou et al. (2015) for the optimal transport problem. The objective is equivalent to

[TABLE]

where $\Phi:=-H$. The derivation, based on a mirror descent step can be found in the appendix. The projection, however, cannot be computed in closed form in general due to the complex geometry of $\mathbb{L}_{m}$.

Ravikumar et al. (2010) proposed using the Bregman method (Bregman, 1966), which has been applied in many fields to solve difficult constrained problems (Benamou et al., 2015; Goldstein and Osher, 2009; Osher et al., 2005, 2010), to compute $\mathcal{P}_{\mathbb{L}_{m}}(\exp({-\eta C}))$ for the inner loop calculation of their proximal algorithm. While the outer loop proximal algorithm can be shown to converge at least linearly, the inner loop rate was not analyzed and the constants (possibly dependent on dimension) were not made clear. Furthermore, the Bregman method is in general inexact, which makes the approximation and the effect on the outer loop unclear (Liu and Ihler, 2013).

5 SMOOTH MESSAGE PASSING

We are interested in analyzing a class of algorithms closely inspired by max-sum diffusion (MSD) as presented by Werner (2007) and the proximal updates of Ravikumar et al. (2010) to solve (Proj) over the $\mathbb{L}_{2}$ polytope. We describe it in detail here, with a few minor modifications and variations to facilitate theoretical analysis. In $\mathbb{L}_{2}$, the constraints occur only over edges between vertices333Written explicitly, the constraints actually occur between any pair of vertices, but these variables play no role in the objective or constraints.. Given an edge $ij\in\mathcal{E}$, we must enforce the constraints prescribed by (7), which is the intersection of the following sets:

[TABLE]

The normalization of the joint distribution $\bm{\mu}_{ij}$ in (b) and (d) is actually a redundant constraint, but it facilitates analysis as we demonstrate in Section 6. For each of these affine constraints, we can compute the Bregman projections in closed form with simple multiplicative updates.

Proposition 1**.**

For a given edge ${ij}\in\mathcal{E}$, the closed-form solutions of the Bregman projections for each of the above individual constraints are given below.

1. (a)

Left consistency: If $\bm{\mu}^{\prime}=\mathcal{P}_{\mathcal{X}_{ij\rightarrow i}}(\bm{\mu})$, then for all $x_{i},x_{j}\in\chi$, $\bm{\mu}_{ij}^{\prime}(x_{i},x_{j})\leftarrow\bm{\mu}_{ij}(x_{i},x_{j})\sqrt{\frac{\bm{\mu}_{i}(x_{i})}{\sum_{x}\bm{\mu}_{ij}(x_{i},x)}}$ and $\bm{\mu}_{i}^{\prime}(x_{i})\leftarrow\bm{\mu}_{i}(x_{i})\sqrt{\frac{\sum_{x}\bm{\mu}_{ij}(x_{i},x)}{\bm{\mu}_{i}(x_{i})}}$. 2. (b)

Left normalization: If $\bm{\mu}^{\prime}=\mathcal{P}_{\mathcal{X}_{ij,i}}(\bm{\mu})$, then for all $x_{i}\in\chi$, $\bm{\mu}_{i}^{\prime}\leftarrow\frac{\bm{\mu}_{i}}{\sum_{x}\bm{\mu}_{i}(x)}$ and $\bm{\mu}_{ij}^{\prime}\leftarrow\frac{\bm{\mu}_{ij}}{\sum_{x_{i},x_{j}}\bm{\mu}_{ij}(x_{i},x_{j})}$. 3. (c)

Right consistency: If $\bm{\mu}^{\prime}=\mathcal{P}_{\mathcal{X}_{ij\rightarrow j}}(\bm{\mu})$, then for all $x_{i},x_{j}\in\chi$, $\bm{\mu}_{ij}^{\prime}(x_{i},x_{j})\leftarrow\bm{\mu}_{ij}(x_{i},x_{j})\sqrt{\frac{\bm{\mu}_{j}(x_{j})}{\sum_{x}\bm{\mu}_{ij}(x,x_{j})}}$ and $\bm{\mu}_{j}^{\prime}(x_{j})\leftarrow\bm{\mu}_{j}(x_{j})\sqrt{\frac{\sum_{x}\bm{\mu}_{ij}(x,x_{j})}{\bm{\mu}_{j}(x_{j})}}$. 4. (d)

Right normalization: If $\bm{\mu}^{\prime}=\mathcal{P}_{\mathcal{X}_{ij,j}}(\bm{\mu})$, then for all $x_{j}\in\chi$, $\bm{\mu}_{j}^{\prime}\leftarrow\frac{\bm{\mu}_{j}}{\sum_{x}\bm{\mu}_{j}(x)}$ and $\bm{\mu}_{ij}^{\prime}\leftarrow\frac{\bm{\mu}_{ij}}{\sum_{x_{i},x_{j}}\bm{\mu}_{ij}(x_{i},x_{j})}$.

These update rules are similar to a number of algorithms throughout the literature on LP relaxations. Notably, they can be viewed as a smoothed version of MSD (Werner, 2007; Kovalevsky and Koval, 1975) in that the updates enforce agreement between variables on the edges and vertices. Nearly identical smoothed updates were also initially proposed by Ravikumar et al. (2010). As in MSD, it is common for message passing schemes derived from LP relaxations to operate on dual objective instead. We presented the primal view here as the Bregman projections lend semantic meaning to the updates and ultimately the stopping conditions in the algorithms. An equivalent dual view is presented in Appendix C.1.

Based on these update rules, we formally outline the algorithms we wish to analyze, which we call edge-based message passing (EMP) for convenience. We consider two variants: EMP-cyclic (Algorithmic 1), which cyclically applies the updates to each edge in each iteration and EMP-greedy (Algorithmic 2), which applies a single projection update to only the edge with the greatest constraint violation in each iteration. We emphasize that these algorithms are not fundamentally new, but our analysis in the next section is our main contribution. EMP-cyclic is the Bregman method, almost exactly the inner loop proposed by Ravikumar et al. (2010). In both variants, $\bm{\mu}^{(1)}$ is defined as the normalized value of $\exp(-\eta C)$. The GreedyEdge operation in EMP-greedy is defined as

[TABLE]

These procedures are then repeated again until the stopping criterion is met, which is that $\bm{\mu}^{(k)}$ is $\epsilon$-close to satisfying the constraint that the joint distributions sum to the marginals for all edges. Both algorithms also conclude with a rounding operation. Any fixed point of EMP must correspond to an optimal $\bm{\mu}_{\eta}^{*}$ (see details in appendix). Computationally, EMP-greedy requires a search over the edges to identify the greatest constraint violation, which can be efficiently implemented using a max-heap (Nutini et al., 2015).

6 THEORETICAL ANALYSIS

We now present our main contribution, a theoretical analysis of EMP-cyclic and EMP-greedy. This result combines two aspects. First, we present a convergence guarantee on the number of iterations sufficient to solve (Proj), satisfying the $\mathbb{L}_{2}$ constraints with $\epsilon>0$ error in $l_{1}$ distance. We note that, in finite iterations, the pseudo-marginals of EMP are not primal feasible in general due to this $\epsilon$-error. We then combine this result with our guarantee on the approximation error in Theorem 1 to show a bound on the number of iterations sufficient to recover the true integral MAP assignment by rounding, assuming the LP is tight and the solution is unique. This holds with sufficient iterations and a sufficiently large regularization constant even though the pseudo-marginals may not be primal feasible. We emphasize that these theorems are a departure from usual convergence rates in the literature (Meshi et al., 2012, 2015). Prior work has guaranteed convergence in objective value to the optimum of the regularized objective (Proj), making it unclear whether the optimal MAP assignment can be recovered, e.g. by rounding. We address this ambiguity in our results.

We begin with the upper bound iterations to obtain $\epsilon$-close solutions, which is the result of two facts which we show. The first is that the updates in Proposition 1 monotonically improve a Lyapunov (potential) function by an amount proportional to the constraint violation as measured via the Hellinger distance. The second is that the difference between the initial and optimal values of the Lyapunov function is bounded.

Let $\deg({\mathcal{G}})$ denote the maximum degree of graph $\mathcal{G}$ and define:

[TABLE]

Theorem 2**.**

For any $\epsilon>0$, EMP is guaranteed to satisfy $\|\bm{\mu}_{ij}\mathbbm{1}-\bm{\mu}_{i}\|_{1}<\epsilon$ and $\|\bm{\mu}_{ij}^{\top}\mathbbm{1}-\bm{\mu}_{j}\|_{1}<\epsilon$ for all $ij\in\mathcal{E}$ in $\lceil\frac{4\mathcal{S}_{0}(\deg(\mathcal{G})+1)}{\epsilon^{2}}\rceil$ iterations for EMP-cyclic and $\lceil\frac{4\mathcal{S}_{0}}{\epsilon^{2}}\rceil$ iterations for EMP-greedy.

Here, $\mathcal{S}_{0}=\min(\|\eta C/d+\exp({-\eta C})\|_{1},S)$. In this theorem, we give our guarantee in terms of $l_{1}$ distance rather than function value convergence. As we will see, this is significant, allowing us to relate this result to Theorem 1 in order to derive the main result. The proof is similar in style to Altschuler et al. (2017). We leave the full proof for EMP-cyclic for the appendix due to a need to handle tedious edge cases, but we state several intermediate results and sketch the proof for EMP-greedy for intuition as it reveals possibly how similar message passing algorithms can be analyzed. We first introduce a Lyapunov function written in terms of dual variables $(\lambda,\xi)$, indexed by the edges and vertices to which they belong in $\mathbb{L}_{2}$. We denote the iteration-indexed dual variables as $(\lambda^{(k)},\xi^{(k)})$. For a given edge $ij\in\mathcal{E}$, constraints enforcing row and column consistency correspond to $\lambda_{ij},\lambda_{ji}\in\mathbb{R}^{m}$, respectively. Normalizing constraints correspond to $\xi_{i},\xi_{j},\xi_{ij}\in\mathbb{R}$. The Lyapunov function, $L(\lambda,\xi)$, is shown in Figure 3.

We note that maximizing $L$ over $(\lambda,\xi)$ satisfies all constraints and yields the solution to (Proj) by first-order optimality conditions. We now present a result that establishes the monotone improvement in $L$ due to the updates in Proposition 1.

Lemma 1**.**

For a given edge $ij\in\mathcal{E}$, let $\bm{\mu}^{\prime}$ and $(\lambda^{\prime},\xi^{\prime})$ denote the updated primal and dual variables after a projection from one of (a)–(d) in Proposition 1. We have the following improvements on $L$. If $\bm{\mu}^{\prime}$ is equal to:

1. (a)

$\mathcal{P}_{\mathcal{X}_{ij\rightarrow i}}(\bm{\mu})$, then $L(\lambda^{\prime},\xi^{\prime})-L(\lambda,\xi)=2h^{2}(\bm{\mu}_{ij}\mathbbm{1},\bm{\mu}_{i})$ 2. (b)

$\mathcal{P}_{\mathcal{X}_{ij,i}}(\bm{\mu})$, then $L(\lambda^{\prime},\xi^{\prime})-L(\lambda,\xi)\geq 0$ 3. (c)

$\mathcal{P}_{\mathcal{X}_{ij\rightarrow j}}(\bm{\mu})$, then $L(\lambda^{\prime},\xi^{\prime})-L(\lambda,\xi)=2h^{2}(\bm{\mu}_{ij}^{\top}\mathbbm{1},\bm{\mu}_{j})$ 4. (d)

$\mathcal{P}_{\mathcal{X}_{ij,j}}(\bm{\mu})$, then $L(\lambda^{\prime},\xi^{\prime})-L(\lambda,\xi)\geq 0$.

This result shows that $L$ improves monotonically after each of the four updates in Proposition 1. Furthermore, at every update, $L$ improves by twice the squared Hellinger distance of the constraint violation between the joint and the marginals.

Lemma 2**.**

Let $\lambda^{*}$, $\xi^{*}$ denote the maximizers of $L$. The difference in function value between the optimal value of $L$ and the first iteration value is upper bounded

[TABLE]

Turning to Theorem 2, the result is obtained by observing that as long as the constraints are violated by an amount $\epsilon>0$ (i.e., the algorithm has not terminated), then the Lyapunov function must improve by a known positive amount at each iteration. We provide a proof sketch for EMP-greedy.

Proof Sketch of Theorem 2 for EMP-greedy.

We now show how to combine the results of Lemma 1 and Lemma 2 to obtain Theorem 2. Let $k^{*}$ be the first iteration such that the termination condition in Algorithm 2 holds with respect to some $\epsilon>0$. Then, for any $k$ satisfying $1\leq k<k^{*}$, we have that $\textsc{GreedyEdge}(\bm{\mu})$ selects $ij$ such that either $\|\bm{\mu}_{ij}\mathbbm{1}-\bm{\mu}_{i}\|_{1}\geq\epsilon$ or $\|\bm{\mu}_{ij}^{\top}\mathbbm{1}-\bm{\mu}_{j}\|_{1}\geq\epsilon$.

Without loss of generality, suppose $\|\bm{\mu}_{ij}\mathbbm{1}-\bm{\mu}_{i}\|_{1}\geq\|\bm{\mu}_{ij}^{\top}\mathbbm{1}-\bm{\mu}_{j}\|_{1}$. Therefore, we have

[TABLE]

where again $h^{2}(\bm{\mu}_{ij}\mathbbm{1},\bm{\mu}_{i})$ denotes the squared Hellinger distance and the last inequality is the Hellinger inequality. Since $\bm{\mu}_{ij}$ and $\bm{\mu}_{i}$ are normalized for each iteration, this inequality is valid. Thus, $L$ improves by $2h^{2}(\bm{\mu}_{ij}\mathbbm{1},\bm{\mu}_{i})$ when $\mathcal{P}_{\mathcal{X}_{ij\rightarrow i}}$ occurs and by a non-negative amount when $\mathcal{P}_{\mathcal{X}_{ij,i}}$ occurs by Lemma 1. Therefore, we can guarantee improvement of at least $\frac{\epsilon^{2}}{4}$ each iteration. Since the optimality gap is at most $\mathcal{S}_{0}$ by Lemma 2, this means the algorithm must terminate in $\lceil\frac{4\mathcal{S}_{0}}{\epsilon^{2}}\rceil$ iterations.∎

We now turn to our main theoretical result. We combine our approximation and iteration convergence guarantees to fully characterize the convergence of EMP for $\mathbb{L}_{2}$ to the optimal MAP assignment when the relaxation is tight and the solution is unique.

Theorem 3**.**

Let $\eta\geq\frac{2\log(16n^{2}d^{2})+16|\mathcal{E}|d^{2}}{\min(\Delta,\frac{1}{128})}$, and $\epsilon^{-1}>(25d\deg(\mathcal{G})|\mathcal{E}|)^{2}\max\left(\eta\|C\|_{\infty},68\right)$. If $\mathbb{L}_{2}$ is tight and $|\mathcal{V}^{*}_{2}|=1$, the EMP algorithm returns a MAP assignment after $\lceil\frac{4\mathcal{S}_{0}(\deg(\mathcal{G})+1)}{\epsilon^{2}}\rceil$ iterations for EMP-cyclic and after $\lceil\frac{4\mathcal{S}_{0}}{\epsilon^{2}}\rceil$ iterations for EMP-greedy.

When $C$ is integral, $\Delta\geq\frac{1}{2}$, yielding a bound of all known parameters. The main technical challenge in producing this result is to relate the termination condition of EMP to the $l_{1}$ distance between $\bm{\mu}^{(k)}$ and $\bm{\mu}^{*}$ (the MAP assignment), as this may lie outside the polytope $\mathbb{L}_{2}$. It does not suffice to provide convergence guarantees in function value as the goal of MAP inference is to produce integral assignments. The proof proceeds in two steps. First we show that $\mathbf{\mu}^{(k)}$ is the entropy-regularized solution to objective $C$ over a “slack” polytope $\mathbb{L}_{2}^{\nu^{(k)}}$. Where the slack vector $\nu^{(k)}$ corresponds to the constraint violations of $\mathbf{\mu}^{(k)}$. We use this characterization to “project” $\mathbf{\mu}^{(k)}$ onto a nearby $\mathbb{L}_{2}$ feasible point $\mathbf{\mu}^{(k)}(2)$. Second, we can use the properties of the primal objective to bound $\mathbf{\mu}^{(k)}(2)$ and $\mathbf{\mu}_{\eta}^{*}$. The proof is in the appendix.

7 NUMERICAL EXPERIMENTS

We illustrate our theoretical results in a practical application of the EMP algorithms. Ravikumar et al. (2010) already gave empirical evidence that the basic EMP-cyclic is competitive with standard solvers. Therefore, the objective of these experiments is to understand how graph and algorithm properties affect approximation (Theorem 1) and convergence (Theorem 2). We consider the family of multi-label Potts models (Wainwright et al., 2005) with $d=3$ labels on $\mathbb{L}_{2}$. For each trial, the cost vector is $C_{i}(x_{i})=\alpha_{i}(x_{i})$, $\forall i,x_{i}$ and

[TABLE]

where the parameters are random $\alpha_{i}(x_{i})\sim\text{Unif}(-0.5,0.5)$ and $\beta_{ij}\sim\text{Unif}\{-0.1,0.1\}$. The graphs considered are structured as $\sqrt{n}\times\sqrt{n}$ grids (Globerson and Jaakkola, 2008; Ravikumar et al., 2010; Erdogdu et al., 2017) and as Erdős-Rényi random graphs with edge probability $p=\frac{1.1\log n}{n}$. To evaluate recovery of the optimal MAP assignment, we first solved each graph with the ECOS LP solver (Domahidi et al., 2013) and selected graphs that were tight. Solving the LP to find the ground-truth was the main computational bottleneck. Further details can be found in Appendix E.

Approximation

In Figure 4, we evaluate the effect of regularization and graph size on the quality of the nearly converged solution from EMP for over $80$ iterations on grids. The box-plots indicate that large choices of $\eta$ often yield the exact MAP solution (cyan and purple). Moderate choices still yield competitive solutions but not optimal for larger graphs (orange and green). Low choices generally give poor solutions with high spread for all graph sizes (red and blue).

Convergence

We then investigate the effects of regularization on convergence for both variants. Figure 5 illustrates the distance of the rounded solution to the optimal MAP solution over projection steps on grids of size $n=2500$. EMP-greedy converges sharply and varying regularization has less of an effect on its convergence rate. Finally, in Figure 6, we look at Erdős-Rényi random graphs to observe the effect of the graph structure for both variants. We considered degree-limited random graphs with $\deg(\mathcal{G})=5$ and $\deg(\mathcal{G})=10$. The figure shows convergence over projection steps for graphs of size $n=400$. For both variants, the convergence rate deteriorates for higher degrees.

8 CONCLUSION

In this paper, we investigated the approximation effects of entropy regularization on MAP inference objectives. We combined these approximation guarantees with a convergence analysis of an edge-based message passing algorithm that solves the regularized objective to derive guarantees on the number of iterations sufficient to recover the true MAP assignment. We also showed empirically the effect of regularization and graph propertise on both the approximation and convergence. In future work, we wish to extend the analyses and proof techniques to higher order polytopes and general block-coordinate minimization algorithms.

Acknowledgements

We thank the anonymous reviewers and Marco Pavone for their invaluable feedback.

Appendix A Bregman Projection Derivation

The objective (Reg) can be equivalently interpreted as a Bregman projection. This interpretation has been explored by Ravikumar et al. (2010) as a basis for proximal updates and also Benamou et al. (2015) for the optimal transport problem. Here, we review the transformation because it is central to the algorithm of Ravikumar et al. (2010), upon which our main theoretical results are based.

By definition of the Bregman projection with respect to the negative entropy, $\Phi=-H$, we have

[TABLE]

where $\mathbbm{1}$ is a vector of ones of the same size as the marginal vector and $=_{+}$ denotes the two sides are equal up to a constant. Substituting this into (Reg) and multiplying through by $\eta$ yields the objective:

[TABLE]

Note the similarity to a projected mirror descent update over $\mathbb{L}_{m}$ starting from $\mathbbm{1}$ (Nemirovsky and Yudin, 1983; Bubeck, 2015). Using this insight and performing a single gradient update in the dual, we can transform the problem into a single Bregman projection of the vector. The unprojected marginal vector $\bm{\mu}^{\prime}$ satisfies

[TABLE]

where $\nabla\Phi(\bm{\mu})=-\nabla H(\bm{\mu})=\log\bm{\mu}$ is the dual map and $(\nabla\Phi)^{-1}(\bm{\mu})=\nabla\Phi^{*}(\bm{\mu})=\exp(\bm{\mu})$ is the inverse dual map. We have $\bm{\mu}^{\prime}=\exp(-\eta C)$ and the solution to the mirror descent update is $\mathcal{P}_{\mathbb{L}_{m}}(\exp(-\eta C))$. Therefore it is sufficient to solve the following Bregman projection problem:

[TABLE]

The projection, however, cannot be computed in closed form due to the complex geometry of $\mathbb{L}_{m}$. Sinkhorn-like algorithms such as those used in Cuturi (2013) are unavailable because the transportation polytopes $\mathcal{U}_{d}(\bm{\mu}_{i},\bm{\mu}_{j})$ are dependent on variables $\bm{\mu}_{i}$ and $\bm{\mu}_{j}$ which are also involved in the projection operation.

Appendix B Derivation of EMP Update Rules

We present the derivations of the update rules similar to Ravikumar et al. (2010) for a given edge $ij\in\mathcal{E}$ based on the Bregman projections onto the individual constraint sets $\mathcal{X}_{ij\rightarrow i}$, $\mathcal{X}_{ij,i}$, $\mathcal{X}_{ij\rightarrow j}$, $\mathcal{X}_{ij,j}$. We refer the reader to Ravikumar et al. (2010) for the original algorithm and derivation. We derive only the first two projections; the last two can be found by exchanging the indices.

• (a)

For the projection $\bm{\mu}^{\prime}=\mathcal{P}_{\mathcal{X}_{ij\rightarrow i}}(\bm{\mu})$, where

[TABLE]

there are no constraints on any edges or vertices other than $ij$ and $i$. Therefore, $\forall k\neq i$, $\bm{\mu}_{k}^{\prime}=\bm{\mu}_{k}$. Similarly, $\forall k\ell\neq ij$, $\bm{\mu}_{k\ell}^{\prime}=\bm{\mu}_{k\ell}$.

The Lagrangian of the projection is given in terms of primal variables $\bm{\mu}$ and dual variables $\alpha$:

[TABLE]

By the first-order optimality condition, the primal solution in terms of the dual variables is

[TABLE]

Substituting this solution back in to the Lagrangian, we have

[TABLE]

Again, by the first-order optimality condition, the dual solution is

[TABLE]

Substituting this value for $\alpha^{*}$ into the primal solution yields the desired result.

• (b)

Again, for the projection onto

[TABLE]

only $\bm{\mu}_{i}$ and $\bm{\mu}_{ij}$ are affected. $\mathcal{X}_{ij,i}$ enforces that the variables $\bm{\mu}_{ij}$ and $\bm{\mu}_{i}$ each sum to one. It is well known and easy to show that the Bregman projection with respect to the negative entropy is simply the $\bm{\mu}_{ij}$ and $\bm{\mu}_{i}$ normalized by their sums. This normalization can also be written as a multiplicative update of the same form by observing that

[TABLE]

where $\xi_{ij}^{*}=\log\sum_{x_{i},x_{j}}\bm{\mu}_{ij}(x_{i},x_{j})$ and $\xi_{i}^{*}=\log\sum_{x_{i}}\bm{\mu}_{i}(x_{i})$. Again, these can be derived via the Lagrangian.

Appendix C Extensions of EMP

C.1 Dual EMP

We may also equivalently interpret the multiplicative updates in Algorithm 1 and Algorithm 2 as additive updates of the dual variables. The dual interpretation is consistent with past work in dual MAP algorithms (Sontag et al., 2011) and may be more practical to avoid numerical issues in implementation. Instead of tracking the primal variables $\bm{\mu}$, we track a sum of the dual variables with $\zeta$ for each vertex and edge. Enforcing consistency between a given joint distribution and its marginals in (a) yields updated dual variable sums

[TABLE]

where again $\alpha^{*}(x_{i})=\frac{1}{2}\log\frac{\sum_{x_{j}}\bm{\mu}_{ij}(x_{i},x_{j})}{\bm{\mu}_{i}(x_{i})}$. The same is done for the vertex $j$ in (c) with indices exchanged. The normalization step in (b) yields

[TABLE]

where $\xi_{ij}^{*}=\log\sum_{x_{i},x_{j}}\bm{\mu}_{ij}(x_{i},x_{j})$ and $\xi_{i}^{*}=\log\sum_{x_{i}}\bm{\mu}_{i}(x_{i})$. Again, the same is done for (d). The primal marginal vector is recovered with

[TABLE]

We will later make explicit the dual formulation as it will aid in the theoretical analysis.

C.2 Clique Constraints

The version of EMP presented in the paper is for the $\mathbb{L}_{2}$ local polytope, which enforces only pairwise consistency among the variables with edges, but this can be fairly easily extended. In this section, we discuss higher order pseudo-marginals and their constraints. Consider the polytope that enforces consistency on all subsets of $\mathcal{V}$ of size $k$ and below, denoted by $\mathcal{C}$. We use the notation of Meshi et al. (2012). The constraint set is written as

[TABLE]

where $x_{c\setminus i}$ denotes a marginalization over all variables except $i$. For convenience, we may also now account for higher-order interactions in the model itself:

[TABLE]

The projection operation in (Proj) is the same for $C=-\theta$. Analogous update rules to Proposition 1 can derived with exactly the same procedure. For a given subset $c\in\mathcal{C}$ and vertex $i$, we have that $\bm{\mu}^{\prime}=\mathcal{P}_{ij\rightarrow i}(\bm{\mu})$ constitutes the update

[TABLE]

The normalization updates are identical as well. As in the presented EMP algorithm, we can design greedy and cyclic algorithms around these update equations. The theoretical analysis in Section 6 will focus on the case with edges only. We leave the general analysis of $\mathbb{L}_{\mathcal{C}}$ for future work.

Appendix D Omitted Proofs and Derivations from Section 6

D.1 Derivation of the Lyapunov function (8)

For convenience, $L$ is restated here:

[TABLE]

The Lagrangian of (Proj) with primal variables $\bm{\mu}$ and dual variables $(\lambda,\xi)$ can be written as

[TABLE]

where

[TABLE]

The partial derivatives with respect to $\bm{\mu}_{ij}(x_{i},x_{j})$ and $\bm{\mu}_{i}(x)$ are given by

[TABLE]

Setting the derivatives to zero gives the solution $\bm{\mu}$ in terms of the dual variables:

[TABLE]

By substituting $\bm{\mu}$ in $\mathcal{L}$, we obtain the Lyapunov function $L$.

D.2 Proof of Lemma 1

In this section we prove Lemma 1. We restate the result for the reader’s convenience.

Lemma 3**.**

For a given edge $ij\in\mathcal{E}$, let $\bm{\mu}^{\prime}$ and $(\lambda^{\prime},\xi^{\prime})$ denote the updated primal and dual variables after a projection from one of (a)–(d) in Proposition 1. We have the following improvements on $L$. If $\bm{\mu}^{\prime}$ is equal to:

1. (a)

$\mathcal{P}_{\mathcal{X}_{ij\rightarrow i}}(\bm{\mu})$, then $L(\lambda^{\prime},\xi^{\prime})-L(\lambda,\xi)=2h^{2}(\bm{\mu}_{ij}\mathbbm{1},\bm{\mu}_{i})$ 2. (b)

$\mathcal{P}_{\mathcal{X}_{ij,i}}(\bm{\mu})$, then $L(\lambda^{\prime},\xi^{\prime})-L(\lambda,\xi)\geq 0$ 3. (c)

$\mathcal{P}_{\mathcal{X}_{ij\rightarrow j}}(\bm{\mu})$, then $L(\lambda^{\prime},\xi^{\prime})-L(\lambda,\xi)=2h^{2}(\bm{\mu}_{ij}^{\top}\mathbbm{1},\bm{\mu}_{j})$ 4. (d)

$\mathcal{P}_{\mathcal{X}_{ij,j}}(\bm{\mu})$, then $L(\lambda^{\prime},\xi^{\prime})-L(\lambda,\xi)\geq 0$.

Proof.

Let $L$ and $L^{\prime}$ denote the values of the Lyapunov function before and after the projection in each case.

(a)  Due to the projection $\bm{\mu}^{\prime}=\mathcal{P}_{\mathcal{X}_{ij\rightarrow i}}(\bm{\mu})$, only $\bm{\mu}_{ij}$ and $\bm{\mu}_{i}$ change values.

[TABLE]

(b)  Due to the projection $\bm{\mu}^{\prime}=\mathcal{P}_{\mathcal{X}_{ij\rightarrow i}}(\bm{\mu})$ change, again only $\bm{\mu}_{ij}$ and $\bm{\mu}_{i}$, but they are simply normalized. From the derivation of the updates, we can see that only dual variables $\xi_{i}$ and $\xi_{ij}$ are updated in order for the normalization to occur. We have, from the update rule in Proposition 1

[TABLE]

The improvement on the Lyapunov function can then be written as

[TABLE]

where the second equality uses the fact that $\bm{\mu}_{ij}^{\prime}$ and $\bm{\mu}_{i}$ both sum to one. This last expression can be shown to be non-negative by recognizing the classical inequality $x-\log x-1\geq 0$ for all $x>0$.

(c)  The proof of improvement is identical to (a); however, we replace vertex $i$ with $j$ and all row sums $\bm{\mu}_{ij}\mathbbm{1}$ with column sum $\bm{\mu}_{ij}^{\top}\mathbbm{1}$.

(d)  The proof of improvement is identical to (b), but we replace $i$ with $j$ for the vertex marginal normalization.

∎

D.3 Fixed points of EMP

We start this section by noting that all fixed points of EMP correspond to valid (constraint satisfying) primal solutions and therefore must equal global optima of the dual function.

First note that any fixed point of EMP corresponds to a candidate solution all whose constraints are satisfied. Indeed, at optimality $\lambda^{*},\xi^{*}$ satisfy:

[TABLE]

with $\bm{\mu}^{*}_{\eta}\in\mathbb{L}_{2}$. Since all constraints are satisfied, for all projection types $\mathcal{P}$ in Lemma 1, $\mathcal{P}(\bm{\mu}_{\eta}^{*})=\bm{\mu}_{\eta}^{*}$.

For the converse, we proceed by contradiction. Let $\bm{\mu}$ be a fixed point of EMP. As such, all the normalization constraints (ensuring the edge and node distributions each sum to one) must be satisfied. Assume then that a constraint of type (a) or (c) is not satisfied. Without loss of generality let $ij\rightarrow i$ be the unsatisfied constraint. As a consequence of 1, the Lyapunov objective can be strictly increased by performing the corresponding Bregman projection, and therefore EMP couldn’t have possibly be at a fixed point. We summarize these observations in the following proposition:

Proposition 2**.**

All maxima of $L(\lambda,\xi)$ are fixed points of EMP and all fixed points of EMP are maxima of $L(\lambda,\xi)$.

D.4 Proof of Lemma 2

In this section we prove Lemma 2, we restate it here for readability:

Lemma 4**.**

Let $\lambda^{*}$, $\xi^{*}$ denote the maximizers of $L$. The difference in function value between the optimal value of $L$ and the value at the first iteration is upper bounded as

[TABLE]

Proof.

We start by showing the upper bound:

[TABLE]

We have that $(\lambda,\xi)=(0,0)$ when $\bm{\mu}=e^{-\eta C}$ before any updates to the primal variables. By Lemma 1, $L(0,0,0)\leq L(\lambda^{(1)},\xi^{(1)})$. Then we have

[TABLE]

We may establish an upper bound on $L(\lambda^{*},\xi^{*})$ by finding a feasible point in the primal objective (Proj). It is easy to verify that $\bm{\mu}$ is in $\mathbb{L}_{2}$ if $\forall ij\in\mathcal{E}$ and $\forall i\in\mathcal{V}$, $\bm{\mu}_{ij}(x_{i},x_{j})=\frac{1}{d^{2}}$ and $\bm{\mu}_{i}(x_{i})=\frac{1}{d}$. With this choice of $\bm{\mu}$, the value of (Proj) is

[TABLE]

where $\mathbb{E}_{U}$ denotes the uniform distribution. where the last inequality follows from the fact that $C_{ij}(0,0)=C_{i}(0)=0$. Therefore,

[TABLE]

We now proceed to show the following (direct) bound on $L(\lambda^{*},\xi^{*})-L(\lambda^{(1)},\xi^{(1)})$:

[TABLE]

We work under the assumption that at any time $k$, all the component distributions of $\bm{\mu}^{(k)}$ are normalized so its entries sum to $1$. Notice that in this case

[TABLE]

If we initialize our algorithm to $\lambda^{(1)}=0$, and $\xi^{(1)}$ be the normalization factors corresponding to this choice of $\lambda$, then

[TABLE]

Notice that at optimality $\lambda^{*},\xi^{*}$, for all $ij\in\mathcal{E}$ and, for all $x_{i},x_{j}$,

[TABLE]

And for all $i\in\mathcal{V}$ and for all $x$,

[TABLE]

Therefore, for all $ij\in\mathcal{E}$ and for all $x_{i},x_{j}$:

[TABLE]

For all $i\in\mathcal{V}$ and for all $x$:

[TABLE]

Summing Equations (14) and (15) over all $ij\in\mathcal{E}$, $i\in\mathcal{V}$ and $x_{i},x_{j},x\in\chi$ yields:

[TABLE]

And, therefore,

[TABLE]

Notice that the RHS of the equation above is positive since: $\sum_{i=1}^{\ell}\exp(a_{i})\geq\frac{1}{\ell}\sum_{i=1}^{\ell}\exp(a_{i})\geq\exp\left(\frac{\sum_{i=1}^{\ell}a_{i}}{\ell}\right)$ for all $\ell\in\mathbb{N}$ and all $a_{1},\cdots,a_{\ell}\in\mathbb{R}$. Combining Equations (13) and (17) and the observation that $L(0,0)\leq L(\lambda^{(1)},\xi^{1})$ (by virtue of Lemma 1) we obtain the final result. ∎

In the case when all entries of $C$ are positive it may be the case that $S\gg\|\exp\left(-\eta C\right)\|_{1}$.

D.5 Complete Proof of Theorem 2

In this section, we will complete the proof of Theorem 2 by handling the case of EMP-cyclic. We require two additional technical lemmas on the $l_{1}$ distance between updated variables. We will use $r(\cdot)$ and $c(\cdot)$ to denote row and column sums respectively of joint distribution matrices.

Lemma 5**.**

Let $a,b\in\Sigma_{d}$ be two points in the simplex and let $p\in\mathbb{R}^{d}_{+}$ s.t. $\min(a_{i},b_{i})\leq p_{i}\leq\max(a_{i},b_{i})$ for all $1\leq i\leq d$. Let $c\in\Sigma_{d}$ defined as $c=\frac{p}{\sum_{i}p_{i}}$. Then:

[TABLE]

Proof.

We only need to prove that $\|a-c\|\leq\|a-b\|_{1}$. From $\min(a_{i},b_{i})\leq p_{i}\leq max(a_{i},b_{i})$ we obtain:

[TABLE]

Let $t=\frac{1}{\sum_{i}p_{i}}$. The following relationships hold:

[TABLE]

Note that

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

The result follows.

∎

Let $A\in\Sigma_{d\times d}$ with elements $a_{ij}$ be a matrix representing joint distribution probabilities. For $p=\begin{bmatrix}p_{1}&\ldots&p_{d}\end{bmatrix}^{\top}\in\Sigma_{d}$, define

[TABLE]

where $z$ is a normalization term, such that the new probabilities matrix sums to one. The notation $r(A)_{i}$ denotes the $i$th element of row sum vector $r(A)$.

Lemma 6**.**

The following inequality holds on the difference between $A$ and $\widetilde{A}$:

[TABLE]

Proof.

[TABLE]

∎

This proof of Theorem 2 relies heavily on the primal and dual variables at given times throughout the algorithm. As such, it is necessary to define precise notation for these temporal events. We note that there are two loops in the algorithm: an outer loop that controls the iterations and an inner one that loops over all edges in $\mathcal{E}$. The outer loop’s current iteration is given by $k\geq 0$, as defined and updated in Algorithm 1. We denote the current step of the inner loop by $t$ where $1\leq t\leq 4|\mathcal{E}|$. This is due to the fact that there are four projections for each edge ($\mathcal{X}_{ij\rightarrow i}$, $\mathcal{X}_{ij,i}$, $\mathcal{X}_{ij\rightarrow j}$, and $\mathcal{X}_{ij,j}$) in one full iteration for $\mathbb{L}_{2}$. Thus the algorithm alternates between enforcing consistency between an edge and vertex and normalizing the local distributions.

The value of $\bm{\mu}$ at iteration $k$ and step $t$ within iteration $k$ is denoted by $\bm{\mu}^{(k,t)}$. For example, at the very start of the algorithm, we are at iteration $k=1$ and step $t=1$ with initial value $\bm{\mu}^{(1,1)}$, which is equal to $\exp(-\eta C)$ with normalized vertex marginal and edge joint distributions. The constraint set onto which a projection is made at $t$ in any iteration is denoted by $\mathcal{X}^{(t)}$. Note that we drop $k$ in the constraint set notation because the order in which the projections occur is always the same.

Proof of Theorem 2.

Let $k^{*}$ be the first iteration such that the termination condition in Algorithm 1 with respect to $\epsilon$ is met. For $k$ such that $1\leq k\leq k^{*}$, there exists $ij\in\mathcal{E}$ such that $\|r(\bm{\mu}_{ij}^{(k,1)})-\bm{\mu}_{i}^{(k,1)}\|_{1}\geq\epsilon$ or $\|c(\bm{\mu}_{ij}^{(k,1)})-\bm{\mu}_{j}^{(k,1)}\|_{1}\geq\epsilon$.

First consider the case where $\|c(\bm{\mu}_{ij}^{(k,1)})-\bm{\mu}_{j}^{(k,1)}\|_{1}\geq\epsilon$. Let $t$ be chosen such that $\mathcal{X}^{(t)}=\mathcal{X}_{ij\rightarrow j}$. Note that $\bm{\mu}_{j}^{(k,t)}$ can move within the $\epsilon$-ball of $c(\bm{\mu}_{ij}^{(k,t)})$ between times $1$ and $t$ of the $k$th iteration due to earlier projections involving vertex $j$. However, $\bm{\mu}_{ij}^{(k,t^{\prime})}=\bm{\mu}_{ij}^{(k,1)}$ for all $t^{\prime}\leq t-2$ because it is only updated at step $t-2$ where $\mathcal{X}^{(t-2)}=\mathcal{X}_{ij\rightarrow i}$. Then, by repeatedly applying the triangle inequality, we have

[TABLE]

where $\mathcal{T}^{(t)}_{j,r}$ and $\mathcal{T}^{(t)}_{j,c}$ are sets of times before $t$ where a projection (for row and column consistency, respectively) caused $\bm{\mu}_{j}$ to be updated:

[TABLE]

Therefore, $\bm{\mu}_{j}^{(k,t^{\prime}+2)}$ is the result of enforcing consistency with another edge of $i$ and then normalizing $\bm{\mu}_{j}$. Let $e_{t^{\prime}}$ denote the edge (incident on $j$) onto which projections are occurring at step $t^{\prime}\in\mathcal{T}_{j,r}^{(t)}\cup\mathcal{T}_{j,c}^{(t)}$. From Lemma 5, if $t^{\prime}\in\mathcal{T}_{j,r}{(t)}$, then

[TABLE]

If $t^{\prime}\in\mathcal{T}_{j,c}^{(t)}$, then

[TABLE]

Similarly, by combining Lemma 5 and Lemma 6, we have

[TABLE]

Note that since the variables are normalized at every even step, they are individually valid probability distributions, and so the Hellinger inequality can be applied. For distributions, $p$ and $q$, the inequality states

[TABLE]

Therefore,

[TABLE]

The last inequality follows from telescoping over all steps in iteration $k$ due to Lemma 1. This proof was for the case when $\|c(\bm{\mu}_{ij}^{(k,1)})-\bm{\mu}_{j}^{(k,1)}\|_{1}\geq\epsilon$. For the case when $\|r(\bm{\mu}_{ij}^{(k,1)})-\bm{\mu}_{i}^{(k,1)}\|_{1}\geq\epsilon$, the procedure is identical except we may ignore the term $\|c(\bm{\mu}_{ij}^{(k,t)})-c(\bm{\mu}_{ij}^{(k,t-2)})\|_{1}$ since $\bm{\mu}_{ij}$ is constant within iteration $k$ until the projection onto $X_{ij\rightarrow i}$. Thus, the improvement lower bound still holds.

Putting these results together with Lemma 2, we see that as long as a single constraint is violated above the $\epsilon$ threshold at the start of an iteration, it is possible to show that the value of $L$ increases by at least $\epsilon^{2}/4(\deg(\mathcal{G})+1)$ during the iteration. This implies that EMP-cyclic terminates in at most $\lceil\frac{4\mathcal{S}_{0}(\deg(\mathcal{G})+1)}{\epsilon^{2}}\rceil$ iterations.

∎

D.6 Proof of Theorem 3

We start by defining a version of $\mathbb{L}_{2}$ with slack vectors. Let $\nu$ be a vector indexed in a similar way as $\bm{\mu}$, where $\{\nu_{ij},\nu_{ji}\}_{ij\in\mathcal{E}}$. We define the slack $\nu$ as $\nu_{ij}=\bm{\mu}_{ij}\mathbbm{1}-\bm{\mu}_{i}$ and $\nu_{ji}=\bm{\mu}_{ij}^{\top}\mathbbm{1}-\bm{\mu}_{j}$. Then we define the slack polytope $\mathbb{L}_{2}^{\nu}$ as

[TABLE]

Notice that by definition the slack vectors $\nu$ satisfy that, for all $ij\in\mathcal{E}$, $\nu_{ij}^{\top}\mathbbm{1}=\nu_{ji}^{\top}\mathbbm{1}=0$. The main difference between $\mathbb{L}_{2}$ and $\mathbb{L}_{2}^{\nu}$ lies in that the joints do not marginalize exactly to the vertex probabilities but do so up to a slack. Consider the entropy-regularized linear program corresponding to $\mathbb{L}_{2}^{\nu}$:

[TABLE]

Introducing the exact same ensemble of dual variables $\lambda,\xi$ as in the Lyapunov function derivation, its dual function equals

[TABLE]

Furthermore, if $\lambda^{*},\xi^{*}$ were a set of optimal dual variables, the optimal primal $\bm{\mu}^{*}$ can be computed via

[TABLE]

They satisfy the same formulae as the problem without slack variables. Since dual optimality is equivalent to primal feasibility, whenever an iterate of EMP satisfies slack of $\nu$, its corresponding primal solution is optimal for (Reg-slack).

We start with a useful manipulation lemma:

Lemma 7**.**

Let $\nu,\nu^{\prime}$ be two slack vectors and let $\bm{\mu}\in\mathbb{L}_{2}^{\nu}$. Assume $\|v^{\prime}\|_{\infty}\leq\frac{1}{2d}$.

If for all $ij\in\mathcal{E}$ and $i\in\mathcal{V}$, $\bm{\mu}_{i}+\nu^{\prime}_{ij}\in\Sigma_{d}$, then there exists a vector $\bm{\mu}^{\prime}\in\mathbb{L}_{2}^{\nu^{\prime}}$ such that

[TABLE] 2. 2.

If $\nu=0$444We do not require that $\bm{\mu}_{i}+\nu^{\prime}_{ij}\in\Sigma_{d}$, then there exists a vector $\bm{\mu}^{\prime}\in\mathbb{L}_{2}^{\nu^{\prime}}$ such that

[TABLE]

Proof.

First we consider the case when for all $ij\in\mathcal{E}$, $\bm{\mu}_{i}+\nu_{ij}^{\prime}$ is a valid distribution (in other words, all its entries are in $[0,1]$ and its values sum to $1$). In this case, we can argue for the existence of $\bm{\mu}^{\prime}$ via the following:

Let $\bm{\mu}^{\prime}_{i}=\bm{\mu}_{i}$ for all $i\in\mathcal{V}$. Let $ij\in\mathcal{E}$ and observe that $\bm{\mu}_{ij}\mathbbm{1}=\bm{\mu}_{i}+\nu_{ij}$ and $\bm{\mu}_{ij}^{\top}\mathbbm{1}=\bm{\mu}_{j}+\nu_{ji}$. We invoke Lemma 7 in Altschuler et al. (2017) to claim the existence of $\bm{\mu}_{ij}^{\prime}$ such that $\bm{\mu}_{ij}^{\prime}\mathbbm{1}=\bm{\mu}_{i}+\nu_{ij}^{\prime}$ and $(\bm{\mu}_{ij}^{\prime})^{\top}\mathbbm{1}=\bm{\mu}_{j}+\nu_{ji}^{\prime}$ and

[TABLE]

Setting $\bm{\mu}^{\prime}$ to be the ensemble with values $\{\bm{\mu}^{\prime}_{i}\}_{i\in\mathcal{V}}$ and $\{\bm{\mu}_{ij}^{\prime}\}_{ij\in\mathcal{E}}$ the result follows.

Now we consider the case when there exist $ij\in\mathcal{E}$ such that $\bm{\mu}_{i}+v_{ij}^{\prime}$ does not lie in the probability simplex. In this case we will have to define $\bm{\mu}^{\prime}_{i}$ different from $\bm{\mu}_{i}$. Consider some $i\in\mathcal{V}$. Let $N(i)$ be the set of neighbouring vertices to $i$ and we abuse notation slightly and use $\nu_{ij}$ for $j\in N(i)$ to denote the slack on $i$ as of the edge marginal shared by $i$ and $j$. We define $\bm{\mu}_{i}^{\prime}$ in the following way:

If $\bm{\mu}_{i}+\nu^{\prime}_{ij}\in\Sigma_{d}$ for all $j\in N(i)$ then let $\bm{\mu}_{i}^{\prime}=\bm{\mu}_{i}$. 2. 2.

Otherwise, let $\{x_{1},\cdots,x_{r}\}\subseteq[d]$ be the entries of $\bm{\mu}_{i}$ such that for all $x_{\tau}\in\{x_{1},\cdots,x_{r}\}$ there exists at least one $j\in N(i)$ for which $\left[\bm{\mu}_{i}+\nu_{ij}^{\prime}\right](x_{\tau})\not\in[0,1]$. Therefore, we must define $\bm{\mu}_{i}^{\prime}$ such that

[TABLE]

which can be done by taking the convex combination of $\mu_{i}$ with the uniform distribution:

[TABLE]

Setting $\theta=d\max_{j}\|\nu_{ij}^{\prime}\|_{\infty}$ guarantees this outcome because we are given that $\|\nu\|_{\infty}\leq\frac{1}{2d}$. Furthermore, we have

[TABLE]

This, in turn, implies $\sum_{i\in V}\|\bm{\mu}_{i}-\bm{\mu}_{i}^{\prime}\|_{1}\leq 2d\|\nu^{\prime}\|_{1}$. Then, we apply the result of Altschuler et al. (2017) again to achieve existence of $\{\bm{\mu}_{ij}^{\prime}\}_{ij\in\mathcal{E}}$ such that

[TABLE]

Summing over these yields $\sum_{ij\in\mathcal{E}}\|\bm{\mu}_{ij}-\bm{\mu}_{ij}^{\prime}\|_{1}\leq 2\|\nu^{\prime}\|_{1}+2d\deg(\mathcal{G})\|\nu^{\prime}\|_{1}$. Therefore $\|\bm{\mu}-\bm{\mu}^{\prime}\|_{1}\leq 6d\deg(\mathcal{G})\|\nu^{\prime}\|_{1}$

∎

We additionally require a similar lemma which allows us to project from one polytope to another while bounding the probabilities away from zero.

Lemma 8**.**

Fix $\tau$ such that $0<\tau\leq\frac{1}{8d^{2}}$ and a slack vector $\nu$ such that $\|\nu\|_{\infty}\leq\frac{1}{4d}$. If $\bm{\mu}\in\mathbb{L}^{\nu}_{2}$, then there exists a vector $\bm{\mu}^{\prime}\in\mathbb{L}_{2}$ such that

[TABLE]

If $\bm{\mu}\in\mathbb{L}_{2}$, then there exists a vector $\bm{\mu}^{\prime}\in\mathbb{L}_{2}^{\nu}$ such that

[TABLE]

Proof.

We address each case individually.

We use the first result from Lemma 7, which yields $\widehat{\bm{\mu}}\in\mathbb{L}_{2}$ such that $\|\bm{\mu}-\widehat{\bm{\mu}}\|_{1}\leq 2\|\nu\|_{1}$. If the probabilities are already bounded below $\tau$, then we are done; however, we must handle the worst case. As in the proof of Lemma 7, we compute a convex combination of $\widehat{\bm{\mu}}$ with the uniform distribution to draw the distribution away from zero values. Define

[TABLE]

where we set $\theta=\tau d^{2}$ which ensures that $\theta\in[0,1]$ and $\bm{\mu}\geq\tau$. Then, note that

[TABLE]

By the triangle inequality, we have

[TABLE] 2. 2.

In the second case, we start by constructing a distribution $\delta$, which is nearly uniform but lives in the slack polytope is and bounded away from zero by at least $\tau$.

For each $i\in\mathcal{V}$, we take $\delta_{i}\in\Sigma_{d}$ to be the uniform distribution where $\delta_{i}(x)=\frac{1}{d}\geq\tau$. Since $\|\nu\|_{\infty}\leq\frac{1}{4d}$, we perturb the uniform distribution with $\nu$ for each $j\in N(i)$, generating $\delta_{i}^{j}:=\delta_{i}+\nu_{ij}$. Again, we are abusing notation slightly by using $\nu_{ij}$ to denote marginalization of edge $ij$ to vertex $i$. Note that $\delta_{i}^{j}\in\Sigma_{d}$, so we can define the product distribution $\delta_{ij}=\delta_{i}^{j}(\delta_{j}^{i})\top\in\mathcal{U}_{d}(\delta_{i}^{j},\delta_{j}^{i})$, which, by construction, marginalizes such that the full vector $\delta$ given by the ensemble $\{\delta_{i}\}_{i\in\mathcal{V}}$ and $\{\delta_{ij}\}_{ij\in\mathcal{E}}$ is in $\mathbb{L}_{2}^{\nu}$. Furthermore, each component can be bounded below as

[TABLE]

Now, as before, we know there exists $\widehat{\bm{\mu}}\in\mathbb{L}_{2}^{2}$ such that $\|\bm{\mu}-\widehat{\bm{\mu}}\|_{1}\leq 6d\deg(\mathcal{G})\|\nu\|_{1}$ from Lemma 7. Therefore, we can take the convex combination of $\bm{\mu}^{\prime}=(1-\theta)\bm{\mu}+\theta\delta$ to get $\bm{\mu}^{\prime}\in\mathbb{L}_{2}^{\nu}$ such that $\bm{\mu}^{\prime}\geq\frac{\theta}{4d^{2}}$ in all entries.

Taking $\theta=4d^{2}\tau\in[0,1]$ ensures that $\bm{\mu}^{\prime}\geq\tau$. Furthermore, the difference can be computed as

[TABLE]

Therefore, we have $\|\widehat{\bm{\mu}}-\bm{\mu}^{\prime}\|_{1}\leq 8(|\mathcal{E}|+n)d^{2}\tau$, which by triangle inequality implies

[TABLE]

∎

We have now the necessary ingredients to prove the first theorem of this section, which provides a bound on the $l_{1}$ distance between the final iterate $\bm{\mu}^{k}$ of Algorithms 1 and 2 and the solution $\bm{\mu}_{\eta}^{*}$ of (Reg). Crucially we analyze these iterates under the assumption all their component distributions $\bm{\mu}^{(k)}_{i}$ for $i\in\mathcal{V}$ and $\bm{\mu}^{(k)}_{ij}$ for ${ij}\in\mathcal{E}$ are normalized.

Theorem 4**.**

Let $\bm{\mu}^{(k)}$ is the $k$th iterate of EMP and let $\nu^{(k)}$ be the slack vector corresponding to $\bm{\mu}^{(k)}$ such that $\|\nu^{(k)}\|_{\infty}\leq\frac{1}{4d}$. In other words,

[TABLE]

Fix $\tau>0$ such that $\tau\leq\frac{1}{8d^{2}}$. Let $\bm{\mu}^{(k)}(2)$ be the pseudo-marginal vector in $\mathbb{L}_{2}$ produced by the first case of Lemma 8 when fed with $\bm{\mu}^{(k)}$ and $\tau$. Then,

[TABLE]

Proof.

By definition $\bm{\mu}^{(k)}\in\mathbb{L}_{2}^{\nu^{(k)}}$. In fact, $\bm{\mu}^{(k)}$ is the optimizer of the following regularized linear program:

[TABLE]

This observation follows because $\bm{\mu}^{(k)}$ is in $\mathbb{L}_{2}^{\nu^{(k)}}$ and its elements can be written as in (24) and (25), thus satisfying dual feasibility.

Recall that after every iteration all the component distributions are normalized. Recall that

[TABLE]

where $\Phi=-H$ is the negative entropy. The point $\bm{\mu}_{\eta}^{*}$ is the optimal point of the information projection $\exp\left(-\eta C\right)$ for points in $\mathbb{L}_{2}$. By the properties of information projections,

[TABLE]

Since for $\bm{\mu}^{(k)}(2)$ and $\bm{\mu}_{\eta}^{*}$, the sum of their entries is the same, by Pinsker’s inequality (applied to each of the component vertex and edge distributions) this in turn implies that

[TABLE]

Let $\bm{\mu}_{\eta}^{*}(\nu^{(k)})$ in $\mathbb{L}_{2}^{\nu^{(k)}}$ be the vector produced by Lemma 8 applied to $\bm{\mu}_{\eta}^{*}\in\mathbb{L}_{2}$. Note that we utilize the existence of $\bm{\mu}_{\eta}^{*}(\nu^{(k)})$ and $\bm{\mu}^{(k)}(2)$ for analysis but we need not actually compute them. Expanding $I$ yields

[TABLE]

Term $A_{2}$ is negative since $\bm{\mu}^{(k)}$ is the optimal point in the slack polytope. Because $\bm{\mu}_{\eta}^{*}(\nu^{(k)})$ and $\bm{\mu}^{(k)}(2)$ were constructed such that all their probabilities are lower bounded by $\tau$, it holds that the entropies are $\log\frac{1}{\tau}$-Lipschitz in $\|\cdot\|_{1}$ Terms $A_{1}$ and $A_{3}$ can be then bounded:

[TABLE]

The result then follow as

[TABLE]

∎

Theorem 4, combined with the EMP algorithm’s optimality condition can provide convergence guarantees for the case when $\mathbb{L}_{2}$ is tight and the solution is unique. We restate the main result, Theorem 3, for readability.

Theorem 5**.**

Let $\eta\geq\frac{2\log(16n^{2}d^{2})+16|\mathcal{E}|d^{2}}{\min(\Delta,\frac{1}{128})}$, and $\epsilon^{-1}>(25d\deg(\mathcal{G})|\mathcal{E}|)^{2}\max\left(\eta\|C\|_{\infty},68\right)$. If $\mathbb{L}_{2}$ is tight and $|\mathcal{V}^{*}_{2}|=1$, the EMP algorithm returns a MAP assignment after $\lceil\frac{4\mathcal{S}_{0}(\deg(\mathcal{G})+1)}{\epsilon^{2}}\rceil$ iterations for EMP-cyclic and after $\lceil\frac{4\mathcal{S}_{0}}{\epsilon^{2}}\rceil$ iterations for EMP-greedy.

Proof.

Let $\bm{\mu}^{(k)}$ be the last internal iterate of the EMP algorithm before rounding. Since the stopping condition has been met, the slack vector $\nu^{(k)}$ corresponding to $\bm{\mu}^{(k)}$ must satisfy $\|\left(\nu^{(k)}\right)_{ij}\|_{1}\leq\epsilon$ for all $ij\in\mathcal{E}$ so that $\|\nu^{(k)}\|_{1}\leq 2|\mathcal{E}|\epsilon$.

Let $\bm{\mu}^{(k)}(2)$ be defined as in Theorem 4 and choose $\tau=\frac{\epsilon}{10(|\mathcal{E}|+n)d^{2}}$555As long as $\epsilon\leq\frac{1}{4d}$ at least, this guarantees $\tau\leq\frac{1}{8d^{2}}$, so we are free to use Theorem 4. Then, the bound from Theorem 4 becomes

[TABLE]

where the last inequality used the fact that $\log(x)\leq n(x^{1/n}-1)$ for $n>0$. Choosing $\epsilon^{-1}>425d^{2}\deg(\mathcal{G})^{2}|\mathcal{E}|^{2}\max\left\{\eta\|C\|_{\infty},68\right\}$ ensures that

[TABLE]

Consequently for all $i\in\mathcal{V}$

[TABLE]

and for all $ij\in\mathcal{E}$

[TABLE]

We also have

[TABLE]

which implies $\|\bm{\mu}^{(k)}(2)-\bm{\mu}^{(k)}\|_{1}\leq\frac{1}{24}$ and $\|\bm{\mu}_{\eta}^{*}-\bm{\mu}^{*}\|_{1}\leq\frac{1}{32}$ (by the condition on $\eta$, see Theorem 1). Putting these inequalities together by triangle inequality,

[TABLE]

For all $i\in\mathcal{V}$. A similar statement holds for all $ij\in\mathcal{E}$:

[TABLE]

Therefore, assuming $\bm{\mu}^{*}$ (the solution of $\mathbb{L}_{2}$) is integral,

[TABLE]

and

[TABLE]

∎

Appendix E Experiment Details

In this section, we provide some additional details for the experiments in Section 7. As mentioned, empirical comparisons between state-of-the-art solvers and EMP-like algorithms have been studied extensively (Meshi et al., 2012; Ravikumar et al., 2010; Werner, 2007; Kappes et al., 2013). For instance, Meshi et al. (2012) found that the regularized star-based message passing algorithms greatly outperform standard optimization techniques such as FISTA and gradient descent, which do not exploit the coordinate structure of the problem.

The primary purpose of these experiments is to understand how the theoretical results in Section 6 manifest in a practical setting. In particular, we would like to understand how the convergence rates, in terms of the ability to round to the solution, behave as a function of the parameters of the problem such as graph size, choice of regularization $\eta$, and connectivity of the graph. In all experiments, we ran an LP solver on the graph in order to obtain the ground-truth MAP assignment. We only considered problems that were tight. The solver specifically is the ECOS solver through a CVXPY wrapper.

E.1 Grid Experiments

As mentioned, our first set of experiments considered solving the MAP problem on $\sqrt{n}\times\sqrt{n}$ grids, totallying $n$ vertices. The vertices were connected by edges to their vertical and horizontal neighbors in the grid. This setting is fairly standard in the literature (Erdogdu et al., 2017; Globerson and Jaakkola, 2008; Ravikumar et al., 2010).

We considered the MAP problem with $d=3$ labels and choose a cost vector $C$ in the family of multi-label Potts models, another well-studied application (Wainwright and Jordan, 2008). Potts models typically have diagonal potentials between edges. That is, we only penalize/reward when the labels on two connected vertices agree. We randomly generated the actual values of the vector. For vertex costs, we chose $C_{i}(x_{i})\sim\text{Unif}(-0.5,0.5)$ and for the edge costs we chose

[TABLE]

where $\beta_{ij}\sim\text{Unif}\{-0.1,0.1\}$. In the approximation results, we ran the algorithms until they had effectively converged after 80 iterations, where each iteration consisted of a full pass over the edges. For EMP-cyclic, this means simply going through all the edges once. For EMP-greedy, one iteration means the opportunity to update each edge exactly once, (although the algorithm will greedily select them in reality). Thus both algorithms update the same number of edges, though their choices will be different. Regardless, we found that 80 iterations was reasonably sufficient to observe the approximation properties. We measured the results in terms of the average Hamming distance between the LP’s solution, which is integral, and the rounded solution returned by the algorithms.

E.2 Random Graph Experiments

While the grid topology offers a consistent platform to evaluate the algorithms, we also considered randomly generated graphs, specifically Erdős-Rényi random graphs. These graphs are constructed by iterating through every pair of the $n$ vertices. Then, an edge is drawn between vertex $i$ and $j$ with probability $p$. Specifically, we chose $p=\frac{1.1\log n}{n}$, which is just large enough that the graph is almost surely connected. We found these to be useful hyperparameter because any lower and the graph would largely be disconnected. Any higher and typically we found the LP was not tight. We chose the same multi-label Potts model for generating the cost vector $C$.

With these experiments, we intended to understand how diverse graph topologies would affect convergence due to randomness. In particular, we restricted the degrees of the graph to $\deg(\mathcal{G})=5,10$ to observe how the algorithms behave on denser graphs.