Barrier Frank-Wolfe for Marginal Inference

Introduction

Markov random fields (MRFs) are used in many areas of computer science such as vision and speech. Inference in these undirected graphical models is generally intractable. Our work focuses on performing approximate marginal inference by optimizing the Tree Re-Weighted (TRW) objective (Wainwright et al., 2005). The TRW objective is concave, is exact for tree-structured MRFs, and provides an upper bound on the log-partition function.

Fast combinatorial solvers for the TRW objective exist, including Tree-Reweighted Belief Propagation (TRBP) (Wainwright et al., 2005), convergent message-passing based on geometric programming (Globerson and Jaakkola, 2007), and dual decomposition (Jancsary and Matz, 2011). These methods optimize over the set of pairwise consistency constraints, also called the local polytope. Sontag and Jaakkola (2007) showed that significantly better results could be obtained by optimizing over tighter relaxations of the marginal polytope. However, deriving a message-passing algorithm for the TRW objective over tighter relaxations of the marginal polytope is challenging. Instead, Sontag and Jaakkola (2007) use the conditional gradient method (also called Frank-Wolfe) and off-the-shelf linear programming solvers to optimize TRW over the cycle consistency relaxation. Rather than optimizing over the cycle relaxation, Belanger et al. (2013) optimize the TRW objective over the exact marginal polytope. Then, using Frank-Wolfe, the linear minimization performed in the inner loop can be shown to correspond to MAP inference.

The Frank-Wolfe optimization algorithm has seen increasing use in machine learning, thanks in part to its efficient handling of complex constraint sets appearing with structured data (Jaggi, 2013; Lacoste-Julien and Jaggi, 2015). However, applying Frank-Wolfe to variational inference presents challenges that were never resolved in previous work. First, the linear minimization performed in the inner loop is computationally expensive, either requiring repeatedly solving a large linear program, as in Sontag and Jaakkola (2007), or performing MAP inference, as in Belanger et al. (2013). Second, the TRW objective involves entropy terms whose gradients go to infinity near the boundary of the feasible set, therefore existing convergence guarantees for Frank-Wolfe do not apply. Third, variational inference using TRW involves both an outer and inner loop of Frank-Wolfe, where the outer loop optimizes the edge appearance probabilities in the TRW entropy bound to tighten it. Neither Sontag and Jaakkola (2007) nor Belanger et al. (2013) explore the effect of optimizing over the edge appearance probabilities.

Although MAP inference is in general NP hard (Shimony, 1994), it is often possible to find exact solutions to large real-world instances within reasonable running times (Sontag et al., 2008; Allouche et al., 2010; Kappes et al., 2013). Moreover, as we show in our experiments, even approximate MAP solvers can be successfully used within our variational inference algorithm. As MAP solvers improve in their runtime and performance, their iterative use could become feasible and as a byproduct enable more efficient and accurate marginal inference. Our work provides a fast deterministic alternative to recently proposed Perturb-and-MAP algorithms (Papandreou and Yuille, 2011; Hazan and Jaakkola, 2012; Ermon et al., 2013).

Contributions. This paper makes several theoretical and practical innovations. We propose a modification to the Frank-Wolfe algorithm that optimizes over adaptively chosen contractions of the domain and prove its rate of convergence for functions whose gradients can be unbounded at the boundary. Our algorithm does not require a different oracle than standard Frank-Wolfe and could be useful for other convex optimization problems where the gradient is ill-behaved at the boundary.

We instantiate the algorithm for approximate marginal inference over the marginal polytope with the TRW objective. With an exact MAP oracle, we obtain the first provably convergent algorithm for the optimization of the TRW objective over the marginal polytope, which had remained an open problem to the best of our knowledge. Traditional proof techniques of convergence for first order methods fail as the gradient of the TRW objective is not Lipschitz continuous.

We develop several heuristics to make the algorithm practical: a fully-corrective variant of Frank-Wolfe that reuses previously found integer assignments thereby reducing the need for new (approximate) MAP calls, the use of local search between MAP calls, and significant re-use of computations between subsequent steps of optimizing over the spanning tree polytope. We perform an extensive experimental evaluation on both synthetic and real-world inference tasks.

Background

Markov Random Fields: MRFs are undirected probabilistic graphical models where the probability distribution factorizes over cliques in the graph. We consider marginal inference on pairwise MRFs with $N$ random variables $X_{1},X_{2},\ldots,X_{N}$ where each variable takes discrete states $x_{i}\in\textrm{VAL}_{i}$. Let $G=(V,E)$ be the Markov network with an undirected edge $\{i,j\}\in E$ for every two variables $X_{i}$ and $X_{j}$ that are connected together. Let $\mathcal{N}(i)$ refer to the set of neighbors of variable $X_{i}$. We organize the edge log-potentials $\theta_{ij}(x_{i},x_{j})$ for all possible values of $x_{i}\in\textrm{VAL}_{i}$, $x_{j}\in\textrm{VAL}_{j}$ in the vector $\bm{\theta}_{ij}$, and similarly for the node log-potential vector $\bm{\theta}_{i}$. We regroup these in the overall vector $\vec{\bm{\theta}}$. We introduce a similar grouping for the marginal vector $\vec{\bm{\mu}}$: for example, $\mu_{i}(x_{i})$ gives the coordinate of the marginal vector corresponding to the assignment $x_{i}$ to variable $X_{i}$.

Tree Re-weighted Objective (Wainwright et al., 2005): Let $Z(\vec{\bm{\theta}})$ be the partition function for the MRF and $\mathcal{M}$ be the set of all valid marginal vectors (the marginal polytope). The maximization of the TRW objective gives the following upper bound on the log partition function:

Frank-Wolfe (FW) Algorithm: In recent years, the Frank-Wolfe (aka conditional gradient) algorithm has gained popularity in machine learning (Jaggi, 2013) for the optimization of convex functions over compact domains (denoted $\mathcal{D}$). The algorithm is used to solve $\min_{\bm{x}\in\mathcal{D}}f(\bm{x})$ by iteratively finding a good descent vertex by solving the linear subproblem:

and then taking a convex step towards this vertex: $\bm{x}^{(k+1)}=(1-\gamma)\bm{x}^{(k)}+\gamma\bm{s}^{(k)}$ for a suitably chosen step-size $\gamma\in$. The algorithm remains within the feasible set (is projection free), is invariant to affine transformations of the domain, and can be implemented in a memory efficient manner. Moreover, the FW gap $g(\bm{x}^{(k)}):=\langle-\nabla f(\bm{x}^{(k)}),\bm{s}^{(k)}-\bm{x}^{(k)}\rangle$ provides an upper bound on the suboptimality of the iterate $\bm{x}^{(k)}$. The primal convergence of the Frank-Wolfe algorithm is given by Thm. $1$ in Jaggi (2013), restated here for convenience: for $k\geq 1$, the iterates $\bm{x}^{(k)}$ satisfy:

where $C_{f}$ is called the “curvature constant”. Under the assumption that $\nabla f$ is $L$-Lipschitz continuousI.e. $\|\nabla{f}(\bm{x})-\nabla{f}(\bm{x}^{\prime})\|_{*}\leq L\|\bm{x}-\bm{x}^{\prime}\|$ for $\bm{x},\bm{x}^{\prime}\in\mathcal{D}$. Notice that the dual norm $\|\cdot\|_{*}$ is needed here. on $\mathcal{D}$, we can bound it as $C_{f}\leq L\operatorname{diam}_{||.||}(\mathcal{D})^{2}$.

Optimizing over Contractions of the Marginal Polytope

Motivation: We wish to (1) use the fewest possible MAP calls, and (2) avoid regions near the boundary where the unbounded curvature of the function slows down convergence. A viable option to address (1) is through the use of correction steps, where after a Frank-Wolfe step, one optimizes over the polytope defined by previously visited vertices of $\mathcal{M}$ (called the fully-corrective Frank-Wolfe (FCFW) algorithm and proven to be linearly convergence for strongly convex objectives (Lacoste-Julien and Jaggi, 2015)). This does not require additional MAP calls. However, we found (see Sec. 5) that when optimizing the TRW objective over $\mathcal{M}$, performing correction steps can surprisingly hurt performance. This leaves us in a dilemma: correction steps enable decreasing the objective without additional MAP calls, but they can also slow global progress since iterates after correction sometimes lie close to the boundary of the polytope (where the FW directions become less informative). In a manner akin to barrier methods and to Garber and Hazan (2013)’s local linear oracle, our proposed solution maintains the iterates within a contraction of the polytope. This gives us most of the mileage obtained from performing the correction steps without suffering the consequences of venturing too close to the boundary of the polytope. We prove a global convergence rate for the iterates with respect to the true solution over the full polytope.

We describe convergent algorithms to optimize $\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ for $\vec{\bm{\mu}}\in\mathcal{M}$. The approach we adopt to deal with the issue of unbounded gradients at the boundary is to perform Frank-Wolfe within a contraction of the marginal polytope given by $\mathcal{M}_{\delta}$ for $\delta\in$, with either a fixed $\delta$ or an adaptive $\delta$.

$\mathcal{M}_{\delta}:=(1-\delta)\mathcal{M}+\delta\,\bm{u}_{0}$, where $\bm{u}_{0}\in\mathcal{M}$ is the vector representing the uniform distribution.

Marginal vectors that lie within $\mathcal{M}_{\delta}$ are bounded away from zero as all the components of $\bm{u}_{0}$ are strictly positive. Denoting $\mathcal{V}^{(\delta)}$ as the set of vertices of $\mathcal{M}_{\delta}$, $\mathcal{V}$ as the set of vertices of $\mathcal{M}$ and $f(\vec{\bm{\mu}}):=-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$, the key insight that enables our novel approach is that:

Therefore, to solve the FW subproblem (3) over $\mathcal{M}_{\delta}$, we can run as usual a MAP solver and simply shift the resulting vertex of $\mathcal{M}$ towards $\bm{u}_{0}$ to obtain a vertex of $\mathcal{M}_{\delta}$. Our solution to optimize over restrictions of the polytope is more broadly applicable to the optimization problem defined below, with $f$ satisfying Prop. 3.3 (satisfied by the TRW objective) in order to get convergence rates.

Solve $\min_{\bm{x}\in\mathcal{D}}f(\bm{x})$ where $\mathcal{D}$ is a compact convex set and $f$ is convex and continuously differentiable on the relative interior of $\mathcal{D}$.

(Controlled growth of Lipschitz constant over $\mathcal{D}_{\delta}$). We define $\mathcal{D}_{\delta}:=(1-\delta)\mathcal{D}+\delta\bm{u}_{0}$ for a fixed $\bm{u}_{0}$ in the relative interior of $\mathcal{D}$. We suppose that there exists a fixed $p\geq 0$ and $L$ such that for any $\delta>0$, $\nabla f(\bm{x})$ has a bounded Lipschitz constant $L_{\delta}\leq L\delta^{-p}\,\,\,\forall\bm{x}\in\mathcal{D}_{\delta}$.

Fixed $\delta$: The first algorithm fixes a value for $\delta$ a-priori and performs the optimization over $\mathcal{D}_{\delta}$. The following theorem bounds the sub-optimality of the iterates with respect to the optimum over $\mathcal{D}$.

Let $f$ satisfy the properties in Prob. 3.2 and Prop. 3.3, and suppose further that $f$ is finite on the boundary of $\mathcal{D}$. Then the use of Frank-Wolfe for $\min_{\bm{x}\in\mathcal{D}_{\delta}}f(\bm{x})$ realizes a sub-optimality over $\mathcal{D}$ bounded as:

where $\bm{x}^{*}$ is the optimal solution in $\mathcal{D}$, $C_{\delta}\leq L_{\delta}\operatorname{diam}_{||.||}(\mathcal{D}_{\delta})^{2}$, and $\omega$ is the modulus of continuity function of the (uniformly) continuous $f$ (in particular, $\omega(\delta)\downarrow 0$ as $\delta\downarrow 0$).

Adaptive $\delta$: The second variant to solve $\min_{\bm{x}\in\mathcal{D}}f(\bm{x})$ iteratively perform FW steps over $\mathcal{D}_{\delta}$, but also decreases $\delta$ adaptively. The update schedule for $\delta$ is given in Alg. 1 and is motivated by the convergence proof. The idea is to ensure that the FW gap over $\mathcal{D}_{\delta}$ is always at least half the FW gap over $\mathcal{D}$, relating the progress over $\mathcal{D}_{\delta}$ with the one over $\mathcal{D}$. It turns out that $\text{FW-gap-}\mathcal{D}_{\delta}=(1-\delta)\text{FW-gap-}\mathcal{D}+\delta\cdot g_{u}(\bm{x}^{(k)})$, where the “uniform gap” $g_{u}(\bm{x}^{(k)})$ quantifies the decrease of the function when contracting towards $\bm{u}_{0}$. When $g_{u}(\bm{x}^{(k)})$ is negative and large compared to the FW gap, we need to shrink $\delta$ (see step 5 in Alg. 1) to ensure that the $\delta$-modified direction is a sufficient descent direction.

We can show that the algorithm converges to the global solution as follows:

For a function $f$ satisfying the properties in Prob. 3.2 and Prop. 3.3, the sub-optimality of the iterates obtained by running the FW updates over $\mathcal{D}_{\delta}$ with $\delta$ updated according to Alg. 1 is bounded as:

A full proof with a precise rate and constants is given in App. D. The sub-optimality $h_{k}:=f(\bm{x}^{(k)})-f(\bm{x}^{*})$ traverses three stages with an overall rate as above. The updates to $\delta^{(k)}$ as in Alg. 1 enable us to (1) upper bound the duality gap over $\mathcal{D}$ as a function of the duality gap in $\mathcal{D}_{\delta}$ and (2) lower bound the value of $\delta^{(k)}$ as a function of $h_{k}$. Applying the standard Descent Lemma with the Lipschitz constant on the gradient of the form $L\delta^{-p}$ (Prop. 3.3), and replacing $\delta^{(k)}$ by its bound in $h_{k}$, we get the recurrence: $h_{k+1}\leq h_{k}-Ch_{k}^{p+2}$. Solving this gives us the desired bound.

Application to the TRW Objective: $\min_{\vec{\bm{\mu}}\in\mathcal{M}}-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ is akin to $\min_{\bm{x}\in\mathcal{D}}f(\bm{x})$ and the (strong) convexity of $-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ has been previously shown (Wainwright et al., 2005; London et al., 2015). The gradient of the TRW objective is Lipschitz continuous over $\mathcal{M}_{\delta}$ since all marginals are strictly positive. Its growth for Prop. 3.3 can be bounded with $p=1$ as we show in App. E.1. This gives a rate of convergence of $O(k^{-1/2})$ for the adaptive-$\delta$ variant, which interestingly is a typical rate for non-smooth convex optimization. The hidden constant is of the order $O(\|\theta\|\cdot|V|)$. The modulus of continuity $\omega$ for the TRW objective is close to linear (it is almost a Lipschitz function), and its constant is instead of the order $O(\|\theta\|+|V|)$.

Algorithm

We detail a few heuristics that aid practicality.

Fast Local Search: Fast methods for MAP inference such as Iterated Conditional Modes (Besag, 1986) offer a cheap, low cost alternative to a more expensive combinatorial MAP solver. We warm start the ICM solver with the last found vertex $\bm{s}^{(k)}$ of the marginal polytope. The subroutine LOCALSEARCH (Alg. 6 in Appendix) performs a fixed number of FW updates to the pseudomarginals using ICM as the (approximate) MAP solver.

Re-optimizing over the Vertices of $\mathcal{M}$ (FCFW algorithm): As the iterations of FW progress, we keep track of the vertices of the marginal polytope found by Alg. 2 in the set $V$. We make use of these vertices in the CORRECTION subroutine (Alg. 5 in Appendix) which re-optimizes the objective function over (a contraction of) the convex hull of the elements of $V$ (called the correction polytope). $\bm{x}^{(0)}$ in Alg. 2 is initialized to the uniform distribution which is guaranteed to be in $\mathcal{M}$ (and $\mathcal{M}_{\delta}$). After updating $\bm{\rho}$, we set $\bm{x}^{(0)}$ to the approximate minimizer in the correction polytope. The intuition is that changing $\bm{\rho}$ by a small amount may not substantially modify the optimal $\bm{x}^{*}$ (for the new $\bm{\rho}$) and that the new optimum might be in the convex hull of the vertices found thus far. If so, CORRECTION will be able to find it without resorting to any additional MAP calls. This encourages the MAP solver to search for new, unique vertices instead of rediscovering old ones.

Approximate MAP Solvers: We can swap out the exact MAP solver with an approximate MAP solver. The primal objective plus the (approximate) duality gap may no longer be an upper bound on the log-partition function (black-box MAP solvers could be considered to optimize over an inner bound to the marginal polytope). Furthermore, the gap over $\mathcal{D}$ may be negative if the approximate MAP solver fails to find a direction of descent. Since adaptive-$\delta$ requires that the gap be positive in Alg. 1, we take the max over the last gap obtained over the correction polytope (which is always non-negative) and the computed gap over $\mathcal{D}$ as a heuristic.

Theoretically, one could get similar convergence rates as in Thm. 3.4 and 3.5 using an approximate MAP solver that has a multiplicative guarantee on the gap (line 8 of Alg. 2), as was done previously for FW-like algorithms (see, e.g., Thm. C.1 in Lacoste-Julien et al. (2013)). With an $\epsilon$-additive error guarantee on the MAP solution, one can prove similar rates up to a suboptimality error of $\epsilon$. Even if the approximate MAP solver does not provide an approximation guarantee, if it returns an upper bound on the value of the MAP assignment (as do branch-and-cut solvers for integer linear programs, or Sontag et al. (2008)), one can use this to obtain an upper bound on $\log Z$ (see App. J).

Experimental Results

Setup: The L1 error in marginals is computed as: $\zeta_{\mu}:=\frac{1}{N}\sum_{i=1}^{N}|\mu_{i}(1)-\mu_{i}^{*}(1)|$. When using exact MAP inference, the error in $\log Z$ (denoted $\zeta_{\log Z}$) is computed by adding the duality gap to the primal (since this guarantees us an upper bound). For approximate MAP inference, we plot the primal objective. We use a non-uniform initialization of $\bm{\rho}$ computed with the Matrix Tree Theorem (Sontag and Jaakkola, 2007; Koo et al., 2007). We perform 10 updates to $\bm{\rho}$, optimize $\vec{\bm{\mu}}$ to a duality gap of $0.5$ on $\mathcal{M}$, and always perform correction steps. We use LOCALSEARCH only for the real-world instances. We use the implementation of TRBP and the Junction Tree Algorithm (to compute exact marginals) in libDAI (Mooij, 2010). Unless specified, we compute marginals by optimizing the TRW objective using the adaptive-$\delta$ variant of the algorithm (denoted in the figures as $M_{\delta})$.

MAP Solvers: For approximate MAP, we run three solvers in parallel: QPBO (Kolmogorov and Rother, 2007; Boykov and Kolmogorov, 2004), TRW-S (Kolmogorov, 2006) and ICM (Besag, 1986) using OpenGM (Andres et al., 2012) and use the result that realizes the highest energy. For exact inference, we use Gurobi Optimization (2015) or toulbar2 (Allouche et al., 2010).

Test Cases: All of our test cases are on binary pairwise MRFs. (1) Synthetic 10 nodes cliques: Same setup as Sontag and Jaakkola (2007, Fig. 2), with $9$ sets of $100$ instances each with coupling strength drawn from $\mathcal{U}[-\theta,\theta]$ for $\theta\in\{0.5,1,2,\ldots,8\}$. (2) Synthetic Grids: $15$ trials with $5\times 5$ grids. We sample $\theta_{i}\sim\mathcal{U}$ and $\theta_{ij}\in$ for nodes and edges. The potentials were $(-\theta_{i},\theta_{i})$ for nodes and $(\theta_{ij},-\theta_{ij};-\theta_{ij},\theta_{ij})$ for edges. (3) Restricted Boltzmann Machines (RBMs): From the Probabilistic Inference Challenge 2011.http://www.cs.huji.ac.il/project/PASCAL/index.php (4) Horses: Large ($N\approx 12000$) MRFs representing images from the Weizmann Horse Data (Borenstein and Ullman, 2002) with potentials learned by Domke (2013). (5) Chinese Characters: An image completion task from the KAIST Hanja2 database, compiled in OpenGM by Andres et al. (2012). The potentials were learned using Decision Tree Fields (Nowozin et al., 2011). The MRF is not a grid due to skip edges that tie nodes at various offsets. The potentials are a combination of submodular and supermodular and therefore a harder task for inference algorithms.

On the Optimization of $\mathcal{M}$ versus $\mathcal{M}_{\delta}$

On the Applicability of Approximate MAP Solvers

Horses: See Fig. 2 (right). The models are close to submodular and the local relaxation is a good approximation to the marginal polytope. Our marginals are visually similar to those obtained by TRBP and our algorithm is able to scale to large instances by using approximate MAP solvers.

Related Work for Marginal Inference with MAP Calls

Hazan and Jaakkola (2012) estimate $\log Z$ by averaging MAP estimates obtained on randomly perturbed inflated graphs. Our implementation of the method performed well in approximating $\log Z$ but the marginals (estimated by fixing the value of each random variable and estimating $\log Z$ for the resulting graph) were less accurate than our method (Fig. 1(e), 1(f)).

Discussion

Our work creates a flexible, modular framework for optimizing a broad class of variational objectives, not simply TRW, with guarantees of convergence. We hope that this will encourage more research on building better entropy approximations. The framework we adopt is more generally applicable to optimizing functions whose gradients tend to infinity at the boundary of the domain.

Our method to deal with gradients that diverge at the boundary bears resemblance to barrier functions used in interior point methods insofar as they bound the solution away from the constraints. Iteratively decreasing $\delta$ in our framework can be compared to decreasing the strength of the barrier, enabling the iterates to get closer to the facets of the polytope, although its worthwhile to note that we have an adaptive method of doing so.

Acknowledgements

RK and DS gratefully acknowledge the support of the Defense Advanced Research Projects Agency (DARPA) Probabilistic Programming for Advancing Machine Learning (PPAML) Program under Air Force Research Laboratory (AFRL) prime contract no. FA8750-14-C-0005. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the view of DARPA, AFRL, or the US government.

References

Appendix A Preliminaries

The supplementary material is divided into two parts:

(1) The first part is dedicated to the exposition of the theoretical results presented in the main paper. Section B details the variants of the Frank-Wolfe algorithm that we used and analyzed. Section C gives the proof to Theorem 3.4 (fixed $\delta$) while Section D gives the proof to Theorem 3.5 (adaptive $\delta$). Finally, Section E applies the convergence theorem to the TRW objective and investigates the relevant constants.

(2) The remainder of the supplementary material provides more information about the experimental setup as well as additional experimental results.

A.2 Descent Lemma

The following descent lemma is proved in \citesupbertsekas1999nonlinear (Prop. A24) and is standard for any convergence proof of first order methods. We provide a proof here for completeness. It also highlights the origin of the requirement that we use dual norm pairings between $\bm{x}$ and the gradient of $f(\bm{x})$ (because of the generalized Cauchy-Schwartz inequality).

Rearranging terms, we get the desired bound. ∎

Appendix B Frank-Wolfe Algorithms

In this section, we present the various algorithms that we use to do fully corrective Frank-Wolfe (FCFW) with adaptive contractions over the domain $\mathcal{D}$, as was done in our experiments.

To implement the approximate correction steps in the fully corrective Frank-Wolfe (FCFW) algorithm, we use the Frank-Wolfe algorithm with away steps \citepsupWolfe:1970wy, also known as the modified Frank-Wolfe (MFW) algorithm \citepsupGuelat:1986fq. We give pseudo-code for MFW in Algorithm 3 (taken from (Lacoste-Julien and Jaggi, 2015)). This variant of Frank-Wolfe adds the possibility to do an “away step” (see step 5 in Algorithm 3) in order to avoid the zig zagging phenomenon that slows down Frank-Wolfe when the solution is close to the boundary of the polytope. For a strongly convex objective (with Lipschitz continuous gradient), the MFW was known to have asymptotic linear convergence \citepsupGuelat:1986fq and its global linear convergence rate was shown recently (Lacoste-Julien and Jaggi, 2015), accelerating the slow general sublinear rate of Frank-Wolfe. When performing a correction over the convex hull over a (somewhat small) set of vertices of $\mathcal{D}_{\delta}$, this convergence difference was quite significant in our experiments (MFW converging in a small number of iterations to do an approximate correction vs. FW taking hundreds of iterations to reach a similar level of accuracy). We note that the TRW objective is strongly convex when all the edge probabilities are non-zero (Wainwright et al., 2005); and that it has Lipschitz gradient over $\mathcal{D}_{\delta}$ (but not $\mathcal{D}$).

The gap computed in step 6 of Algorithm 3 is non-standard; it is a sufficient condition to ensure the global linear convergence of the outer FCFW algorithm when using Algorithm 3 as a subroutine to implement the approximate correction step. See Lacoste-Julien and Jaggi (2015) for more details.

The MFW algorithm requires more bookkeeping than standard FW: in addition to the current iterate $\bm{x}^{(k)}$, it also maintains both the active set $\mathcal{S}^{(k)}$ (to search for the “away vertex”) as well as the barycentric coordinates $\bm{\alpha}^{(k)}$ (to know what are the away step-sizes that ensure feasibility – see step 13) i.e. $\bm{x}^{(k)}=\sum_{\bm{v}\in\mathcal{S}^{(k)}}\alpha^{(k)}_{\bm{v}}\bm{v}$.

B.2 Fully Corrective Frank-Wolfe (FCFW) with Adaptive-δ𝛿\delta

We give in Algorithm 4 the pseudo-code to perform fully corrective Frank-Wolfe optimization over $\mathcal{D}$ by iteratively optimizing over $\mathcal{D}_{\delta}$ with adaptive-$\delta$ updates. If $\delta$ is kept constant (skipping step 10), then Algorithm 4 implements the fixed $\delta$ variant over $\mathcal{D}_{\delta}$. We describe the algorithm as maintaining the correction set of atoms $V^{(k+1)}$ over $\mathcal{D}$ (rather than $\mathcal{D}_{\delta}$), as $\delta$ is constantly changing. One can easily move back and forth between $V^{(k+1)}$ and its contraction $V_{\delta}=(1-\delta^{(k)})V^{(k+1)}+\delta^{(k)}\bm{u}_{0}$, and so we note that an efficient implementation might work with either representation cheaply (for example, by storing only $V^{(k+1)}$ and $\delta$, not the perturbed version of the correction polytope). The approximate correction over $V_{\delta}$ is implemented using the MFW algorithm described in Algorithm 3, which requires a barycentric representation $\bm{\alpha}^{(k)}$ of the current iterate $\bm{x}^{(k)}$ over the correction polytope $V_{\delta}$. Our notation in Algorithm 4 uses the elements of $\mathcal{V}$ as indices, rather than their contracted version; that is, we maintain the property that $\bm{x}^{(k)}=\sum_{\bm{v}\in\mathcal{V}}\alpha_{\bm{v}}^{(k)}[(1-\delta^{(k)})\bm{v}+\delta^{(k)}\bm{u}_{0}]$. As $V_{\delta}$ changes when $\delta$ changes, we need to update the barycentric representation of $\bm{x}^{(k)}$ accordingly – this is done in step 11 with the following equation. Suppose that we decrease $\delta$ to $\delta^{\prime}$. Then the old coordinates $\bm{\alpha}$ can be updated to new coordinates $\bm{\alpha}^{\prime}$ for the new contraction polytope as follows:

This ensures that $\sum_{\bm{v}}\alpha_{\bm{v}}\bm{v}_{(\delta)}=\sum_{\bm{v}}\alpha^{\prime}_{\bm{v}}\bm{v}_{(\delta^{\prime})}$, where $\bm{v}_{(\delta)}:=(1-\delta)\bm{v}+\delta\bm{u}_{0}$, and that the coordinates form a valid convex combination (assuming that $\delta^{\prime}\leq\delta$), as can be readily verified.

Appendix C Bounding the Sub-optimality for Fixed δ𝛿\delta Variant

The pseudocode for optimizing over $\mathcal{D}_{\delta}$ for a fixed $\delta$ is given in Algorithm 4 (by ignoring the step 10 which updates $\delta$). It is stated with a stopping criterion $\epsilon$, but it can alternatively be run for a fixed number of $K$ iterations. The following theorem bounds the suboptimality of the iterates with respect to the true optimum $\bm{x}^{*}$ over $\mathcal{D}$. If one can compute the constants in the theorem, one can choose a target contraction amount $\delta$ to guarantee a specific suboptimality of $\epsilon^{\prime}$; otherwise, one can choose $\delta$ using heuristics. Note that unlike the adaptive-$\delta$ variant, this algorithm does not converge to the true solution as $K\rightarrow\infty$ unless $\bm{x}^{*}$ happens to belong to $\mathcal{D}_{\delta}$. But the error can be controlled by choosing $\delta$ small enough.

Let $f$ satisfy the properties in Problem 3.2 and suppose its gradient is Lipschitz continuous on the contractions $\mathcal{D}_{\delta}$ as in Property 3.3. Suppose further that $f$ is finite on the boundary of $\mathcal{D}$.

Then $f$ is uniformly continuous on $\mathcal{D}$ and has a modulus of continuity function $\omega$ quantifying its level of continuity, i.e. $|f(\bm{x})-f(\bm{x}^{\prime})|\leq\omega(\|\bm{x}-\bm{x}^{\prime}\|)\,\,\forall\bm{x},\bm{x}^{\prime}\in\mathcal{D}$, with $\omega(\sigma)\downarrow 0$ as $\sigma\downarrow 0$.

Let $\bm{x}^{*}$ be an optimal point of $f$ over $\mathcal{D}$. The iterates $\bm{x}^{(k)}\in\mathcal{D}_{\delta}$ of the FCFW algorithm as described in Algorithm 4 for a fixed $\delta>0$ has sub-optimality over $\mathcal{D}$ bounded as:

where $C_{\delta}\leq\operatorname{diam}(\mathcal{D}_{\delta})^{2}L_{\delta}$. Note that different norms can be used in the definition of $\omega(\cdot)$ and $C_{\delta}$.

Let $\bm{x}^{*}_{(\delta)}$ be an optimal point of $f$ over $\mathcal{D}_{\delta}$. As $f$ has a Lipschitz continuous gradient over $\mathcal{D}_{\delta}$, we can use any standard convergence result of the Frank-Wolfe algorithm to bound the suboptimality of the iterate $\bm{x}^{(k)}$ over $\mathcal{D}_{\delta}$. Algorithm 4 (with a fixed $\delta$) describes the FCFW algorithm which guarantees at least as much progress as the standard FW algorithm (by step 15 and 20a), and thus we can use the convergence result from Jaggi (2013) as already stated in (4): $f(\bm{x}^{(k)})-f(\bm{x}^{*}_{(\delta)})\leq\frac{2C_{\delta}}{(k+2)}$ with $C_{\delta}\leq\operatorname{diam}(\mathcal{D}_{\delta})^{2}L_{\delta}$, where $L_{\delta}$ comes from Property 3.3. This gives the first term in (7). Note that if the function $f$ is strongly convex, then the FCFW algorithm has also a linear convergence rate (Lacoste-Julien and Jaggi, 2015), though we do not cover this here.

Finally, we explain why $f$ is uniformly continuous. As $f$ is a (lower semi-continuous) convex function, it is continuous at every point where it is finite. As $f$ is said to be finite at its boundary (and it is obviously finite in the relative interior of $\mathcal{D}$ as it is continuously differentiable there), then $f$ is continuous over the whole of $\mathcal{D}$. As $\mathcal{D}$ is compact, this means that $f$ is also uniformly continuous over $\mathcal{D}$. ∎

We note that the modulus of continuity function $\omega$ quantifies the level of continuity of $f$. For a Lipschitz continuous function, we have $\omega(\sigma)\leq L\sigma$. If instead we have $\omega(\sigma)\leq C\sigma^{\alpha}$ for some $\alpha\in$, then $f$ is actually $\alpha$-Hölder continuous. We will see in Section E.2 that the TRW objective is not Lipschitz continuous, but it is $\alpha$-Hölder continuous for any $\alpha<1$, and so is “almost” Lipschitz continuous. From the theorem, we see that to get an accuracy of the order $\epsilon$, we would need $(\delta\operatorname{diam}(\mathcal{D}))^{\alpha}<\epsilon$, and thus a contraction of $\delta<\frac{\epsilon^{(1/\alpha)}}{\operatorname{diam}(\mathcal{D})}$.

Appendix D Convergence with Adaptive-δ𝛿\delta

In this section, we show the convergence of the adaptive-$\delta$ FW algorithm to optimize a function $f$ satisfying the properties in Problem 3.2 and Property 3.3 (Lipschitz gradient over $\mathcal{D}_{\delta}$ with bounded growth).

The adaptive update for $\delta$ (given in Algorithm 1) can be used with the standard Frank-Wolfe optimization algorithm or also the fully corrective Frank-Wolfe (FCFW) variant. In FCFW, we ensure that every update makes more progress than a standard FW step with line-search, and thus we will show the convergence result in this section for standard FW (which also applies to FCFW). We describe the FCFW variant with approximate correction steps in Algorithm 4, as this is what we used in our experiments.

We first list a few definitions and lemmas that will be used for the main convergence convergence result given in Theorem D.6. We begin with the definitions of duality gaps that we use throughout this section. The Frank-Wolfe gap is our primary criterion for halting and measuring the progress of the optimization over $\mathcal{D}$. The uniform gap is a measure of the decrease obtainable from moving towards the uniform distribution.

The Frank-Wolfe (FW) gap is defined as: $g(\bm{x}^{(k)}):=\langle-\nabla f(\bm{x}^{(k)}),\bm{s}^{(k)}-\bm{x}^{(k)}\rangle$.

The uniform gap is defined as: $g_{u}(\bm{x}^{(k)}):=\langle-\nabla f(\bm{x}^{(k)}),\bm{u}_{0}-\bm{x}^{(k)}\rangle$.

The FW gap over $\mathcal{D}_{\delta}$ is: $g_{(\delta^{(k)})}(\bm{x}^{(k)}):=\langle-\nabla f(\bm{x}^{(k)}),\bm{s}^{(k)}_{(\delta)}-\bm{x}^{(k)}\rangle$.

The name for the uniform gap comes from the fact that the FW gap over $\mathcal{D}_{\delta}$ can be expressed as a convex combination of the FW gap over $\mathcal{D}$ and the uniform gap:

For iterates progressing as in Algorithm 4 with adaptive update on $\delta$ as given in Algorithm 1, the gap over $\mathcal{D}_{\delta}$ and $\mathcal{D}$ are related as : $g_{(\delta^{(k)})}(\bm{x}^{(k)})\geq\frac{g(\bm{x}^{(k)})}{2}$.

The duality gaps $g(\bm{x}^{(k)})$ and $g_{(\delta^{(k)})}(\bm{x}^{(k)})$ computed as defined in (D.1) during Algorithm 4 are related by equation (8).

(2) When $g_{u}(\bm{x}^{(k)})<0$, from the update rule in lines 5 to 7 in Algorithm 1, we have $\delta^{(k)}\leq\frac{g(\bm{x}^{(k)})}{-4g_{u}(\bm{x}^{(k)})}\implies\delta^{(k)}g_{u}(\bm{x}^{(k)})\geq-\frac{g(\bm{x}^{(k)})}{4}$. Therefore, $g_{(\delta^{(k)})}(\bm{x}^{(k)})\geq\frac{3}{4}g(\bm{x}^{(k)})-\frac{g(\bm{x}^{(k)})}{4}=\frac{g(\bm{x}^{(k)})}{2}$.

Therefore, the gap over $\mathcal{D}_{\delta}$ and $\mathcal{D}$ are related as : $g_{(\delta^{(k)})}(\bm{x}^{(k)})\geq\frac{g(\bm{x}^{(k)})}{2}$. ∎

Another property that we will use in the convergence proof is that $-g_{u}$ is upper bounded for any convex function $f$:Note that on the other hand, $g_{u}(\bm{x})$ might go to infinity as $\bm{x}$ gets close to the boundary of $\mathcal{D}$ as the gradient of $f$ is allowed to be unbounded. Fortunately, we only need an upper bound on $-g_{u}$, not a lower bound.

Let $f$ be a continuously differentiable convex function on the relative interior of $\mathcal{D}$. Then for any fixed $\bm{u}_{0}$ in the relative interior of $\mathcal{D}$, $\exists B$ s.t.

In particular, we can take the finite value:

As $f$ is convex, its directional derivative is a monotone increasing function in any direction. Let $\bm{u}_{0}$ and $\bm{x}$ be points in the relative interior of $\mathcal{D}$; then their gradient exists and we have by the monotonicity property:

This inequality is valid for all $\bm{x}$ in the relative interior of $\mathcal{D}$, and can be extended to the boundary by taking limits (with potentially the RHS become minus infinity, but this is not a problem). Finally, by the definition of the dual norm (generalized Cauchy-Schwartz), we have $\left\langle\nabla f(\bm{u}_{0}),\bm{u}_{0}-\bm{x}\right\rangle\leq\|\nabla f(\bm{u}_{0})\|_{*}\|\bm{u}_{0}-\bm{x}\|\leq\|\nabla f(\bm{u}_{0})\|_{*}\operatorname{diam}_{\|\cdot\|}(\mathcal{D})$. ∎

Finally, we need a last property of Algorithm 4 that allows us to bound the amount of perturbation $\delta^{(k)}$ of the polytope at every iteration as a function of the sub-optimality over $\mathcal{D}$.

Let $B$ be a bound such that $-8g_{u}(\bm{x})\leq B$ for all $\bm{x}\in\mathcal{D}$ (given by Lemma D.3). Then at every stage of Algorithm 4, we have that:

For the last equality, we used the fact that $h_{k}$ is non-increasing since Algorithm 4 decreases the objective at every iteration (using the line-search in step 14). ∎

We now bound the generalization of a standard recurrence that will arise in the proof of convergence. This is a generalization of the technique used in \citetsupteo2007erm (also used in the context of Frank-Wolfe in the proof of Theorem C.4 in Lacoste-Julien et al. (2013)). The basic idea is that one can bound a recurrence inequality by the solution to a differential equation. We provide a detailed proof of the bound for completeness here.

Let $1<a\leq b$. Suppose that $h_{k}$ is any non-negative sequence that satisfies the recurrence inequality:

Then $h_{k}$ is strictly decreasing (unless it equals zero) and can be bounded for $k\geq 0$ as:

Taking the continuous time analog of the recurrence inequality, we consider the differential equation:

Finally, whenever $h_{k}>0$, we have that $h_{k+1}<h_{k}$ from the recurrence inequality, and so $h_{k}$ is strictly decreasing as claimed. ∎

Consider the optimization of $f$ satisfying the properties in Problem 3.2 and Property 3.3. Let $\widetilde{C}:=L\operatorname{diam}_{\|\cdot\|}(\mathcal{D})^{2}$, where $L$ is from Property 3.3. Let $B$ be the upper bound on the negative uniform gap: $-8g_{u}(\bm{x})\leq B$ for all $\bm{x}\in\mathcal{D}$, as used in Lemma D.4 (arising from Lemma D.3). Then the iterates $\bm{x}^{(k)}$ obtained by running the Frank-Wolfe updates over $\mathcal{D}_{\delta}$ with line-search with $\delta$ updated according to Algorithm 1 (or as summarized in a FCFW variant in Algorithm 4), have suboptimality $h_{k}$ upper bounded as:

$h_{k}\leq\left(\frac{1}{2}\right)^{k}h_{0}+\frac{\widetilde{C}}{\delta_{0}^{p}}$ for $k$ such that $h_{k}\geq\max\{B\delta_{0},\frac{2\widetilde{C}}{\delta_{0}^{p}}\}$,

$h_{k}\leq\frac{2\widetilde{C}}{\delta_{0}^{p}}\left[\frac{1}{\frac{1}{4}(k-k_{0})+1}\right]$ for $k$ such that $B\delta_{0}\leq h_{k}\leq\frac{2\widetilde{C}}{\delta_{0}^{p}}$,

$h_{k}\leq\left[\frac{\max(\widetilde{C},B\delta_{0}^{p+1})B^{p}}{\frac{p+1}{\max(8,p+2)}(k-k_{1})+1}\right]^{\frac{1}{p+1}}=O(k^{-\frac{1}{p+1}})$ for $k$ such that $h_{k}\leq B\delta_{0}$,

Let $\bm{x}_{\gamma}:=\bm{x}^{(k)}+\gamma\bm{d}_{k}^{\hskip 0.35002pt\textnormal{FW}}$ with $\bm{d}_{k}^{\hskip 0.35002pt\textnormal{FW}}$ defined in step 12 in Algorithm 4. Note that $\bm{x}_{\gamma}\in\mathcal{D}_{\delta}$ with $\delta=\delta^{(k)}$ for all $\gamma\in$. We apply the Descent Lemma A.1 on this update to get:

We have $L\|\bm{d}_{k}^{\hskip 0.35002pt\textnormal{FW}}\|^{2}\leq\widetilde{C}$ by assumption and $\langle\nabla f(\bm{x}^{(k)}),\bm{d}_{k}^{\hskip 0.35002pt\textnormal{FW}}\rangle=-g_{(\delta)}(\bm{x}^{(k)})$ by definition. Moreover, $\bm{x}^{(k+1)}$ is defined to make at least as much progress than the line-search result $\min_{\gamma\in}f(\bm{x}_{\gamma})$ (line 14 and 15), and so we have:

For the final inequality, we used Lemma D.2 which relates the gap over $\mathcal{D}_{\delta}$ to the gap over $\mathcal{D}$.

Subtracting $f(\bm{x}^{*})$ from both sides and using $g(\bm{x}^{(k)})\geq h_{k}$ by convexity, we get:

Stage 1: The $\min$ in the denominator is $\delta_{0}$ and $h_{k}$ is big: $h_{k}\geq\max\{B\delta_{0},\frac{2\widetilde{C}}{\delta_{0}^{p}}\}$.

Stage 2: The $\min$ in the denominator is $\delta_{0}$ and $h_{k}$ is small: $B\delta_{0}\leq h_{k}\leq\frac{2\widetilde{C}}{\delta_{0}^{p}}$.

Stage 3: The $\min$ in the denominator is $\frac{h_{k}}{B}$, i.e.: $h_{k}\leq B\delta_{0}$.

Since $h_{k}$ is decreasing, once we leave a stage, we no longer re-enter it. The overall strategy for each stage is as follows. For each recurrence that we get, we select a $\gamma^{*}$ that realizes the tightest upper bound on it.

Since we are restricted that $\gamma^{*}\in$, we have to consider when $\gamma^{*}>1$ and $\gamma^{*}\leq 1$. For the former, we bound the recurrence obtained by substituting $\gamma=1$ into (11). For the latter, we substitute the form of $\gamma^{*}$ into the recurrence and bound the result.

We consider the case where $h_{k}\geq B\delta_{0}$. This yields:

The bound is minimized by setting $\gamma^{*}=\frac{h_{k}\delta_{0}^{p}}{2\widetilde{C}}$. On the other hand, the bound is only valid for $\gamma\in$, and thus if $\gamma^{*}>1$, i.e. $h_{k}>\frac{2\widetilde{C}}{\delta_{0}^{p}}$ (stage 1), then $\gamma=1$ will yield the minimum feasible value for the bound. Unrolling the recursion (12) for $\gamma=1$ during this stage (where $h_{l}>\frac{2\widetilde{C}}{\delta_{0}^{p}}$ for $l<k$ as $h_{k}$ is decreasing), we get:

giving the bound for the iterates in the first stage.

We can compute an upper bound on the number of steps it takes to reach a suboptimality of $\frac{2\widetilde{C}}{\delta_{0}^{p}}$ by looking at the minimum $k$ which ensures that the bound in (13) becomes smaller than $\frac{2\widetilde{C}}{\delta_{0}^{p}}$, yielding $k_{\max}=\max(0,\lceil\log_{\frac{1}{2}}\frac{\widetilde{C}}{h_{0}\delta_{0}^{p}}\rceil)$. Therefore, let $k_{0}\leq k_{\max}$ be the first $k$ such that $h_{k}\leq\frac{2\widetilde{C}}{\delta_{0}^{p}}$.

Stage 2

In stage 2, we suppose that $B\delta_{0}\leq h_{k}\leq\frac{2\widetilde{C}}{\delta_{0}^{p}}$. This means that $\gamma^{*}=\frac{h_{k}\delta_{0}^{p}}{2\widetilde{C}}\leq 1$.

Substituting $\gamma=\gamma^{*}$ into (12) yields: $h_{k+1}\leq h_{k}-h_{k}^{2}\frac{\delta_{0}^{p}}{8\widetilde{C}}$.

Using the result of Lemma D.5 with $a=2$, $b=4$ and $C_{0}=\frac{2\widetilde{C}}{\delta_{0}^{p}}$, we get the bound:

It is worthwhile to point out at this juncture that the bound obtained for stage $2$ is the same as the one for regular Frank-Wolfe, but with a factor of $4$ worse due to the factor of $\frac{1}{2}$ in front of the FW gap which appeared due to Lemma D.2.

Stage 3

Here, we suppose $h_{k}\leq B\delta_{0}$. We can compute a bound on the number of steps $k_{1}$ needed get to stage 3 by looking at the number of steps it takes for the bound in stage 2 to becomes less than $B\delta_{0}$:

As before, moving forward, our notation on $k$ represents the number of steps taken after $k_{1}$ steps.

Then, the master inequality (11) becomes:

To simplify the rest of the analysis, we replace $\widetilde{C}B^{p}$ with $F:=\max(B\delta_{0}^{p+1},\widetilde{C})B^{p}$. We then get the bound:

which is minimized by setting $\gamma^{*}:=\frac{h_{k}^{p+1}}{2F}$. Since $F\geq B^{p+1}\delta_{0}^{p+1}$ (by construction) and $h_{k}^{p+1}\leq(B\delta_{0})^{p+1}$ (by the condition to be in stage 3), we necessarily have that $\gamma^{*}\leq 1$. We chose the value of $F$ to avoid having to consider the possibility $\gamma^{*}>1$ as we did in the distinction between stage 1 and stage 2.

Hence, substituting $\gamma=\gamma^{*}$ in (14), we get:

Using the result of Lemma D.5 with $a=p+2$, $b=\max(8,p+2)$ and $C_{0}=F$, we get the bound:

Interestingly, the obtained rate of $O(1/\sqrt{k})$ for $p=1$ (for the TRW objective e.g.) is the standard rate that one would get for the optimization of a general non-smooth convex function with the projected subgradient method (and it is even a lower bound for some class of first-order methods; see e.g. Section 3.2 in \citetsupnesterov2004lectures). The fact that our function $f$ does not have Lipschitz continuous gradient on the whole domain brings us back to the realm of non-smooth optimization. It is an open question whether Algorithm 4 has an optimal rate for the class of functions defined in the assumptions of Theorem D.6.

Appendix E Properties of the TRW Objective

In this section, we explicitly compute bounds for the constants appearing in the convergence statements for our fixed-$\delta$ and adaptive-$\delta$ algorithms for the optimization problem given by:

In particular, we compute the Lipschitz constant for its gradient over $\mathcal{M}_{\delta}$ (Property 3.3), we give a form for its modulus of continuity function $\omega(\cdot)$ (used in Theorem 3.4), and we compute $B$, the upper bound on the negative uniform gap (as used in Lemma D.3).

We first motivate our choice of norm over $\mathcal{M}$. Recall that $\vec{\bm{\mu}}$ can be decomposed into $|V|+|E|$ blocks, with one pseudo-marginal vector $\bm{\mu}_{i}\in\Delta_{\textrm{VAL}_{i}}$ for each node $i\in V$, and one vector $\bm{\mu}_{ij}\in\Delta_{\textrm{VAL}_{i}\textrm{VAL}_{j}}$ per edge $\{i,j\}\in E$, where $\Delta_{d}$ is the probability simplex over $d$ values. We let $c$ be the cliques in the graph (either nodes or edges). From its definition in (2), $f(\vec{\bm{\mu}}):=-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ decomposes as a separable sum of functions of each block only:

where $K_{c}$ is $(1-\sum_{j\in\mathcal{N}(i)}\rho_{ij})$ if $c=i$ and $\rho_{ij}$ if $c=\{i,j\}$. The function $g_{c}$ also decomposes as a separable sum:

We first consider one scalar component of the separable $g_{c}(\bm{\mu}_{c})$ function given in (16) (i.e. for one $\mu_{c}(x_{c})$ coordinate). Its derivative is $K_{c}(1+\log(\mu_{c}(x_{c}))-\theta_{c}(x_{c})$ with second derivative $\frac{K_{c}}{\mu_{c}(x_{c})}$. If $\vec{\bm{\mu}}\in\mathcal{M}_{\delta}$, then we have $\mu_{c}(x_{c})\geq\delta u_{0}(x_{c})=\frac{\delta}{n_{c}}$, where $n_{c}$ is the number of possible values that the assignment variable $x_{c}$ can take. Thus for $\vec{\bm{\mu}}\in\mathcal{M}_{\delta}$, we have that the $x_{c}$-component of $g_{c}$ is Lipschitz continuous with constant $|K_{c}|n_{c}/\delta$. We thus have:

The Lipschitz constant is thus indeed $\frac{L}{\delta}$ with $L:=\sum_{c}|K_{c}|n_{c}$. Let us first consider the sum for $c\in V$; we have $K_{i}=1-\sum_{j\in\mathcal{N}(i)}\rho_{ij}$. Thus:

Here we used the fact that $\rho_{ij}$ came from the marginal probability of edges of spanning trees (and so with $|V|-1$ edges). Similarly, we have $\sum_{ij\in E}|K_{ij}|\leq|V|$. Combining these we get:

E.2 Modulus of Continuity Function

We begin by computing a modulus of continuity function for $-x\log x$ with an additive linear term.

Let $g(x):=-Kx\log x+\theta x$. Consider $x,x^{\prime}\in$ such that $|x-x^{\prime}|\leq\sigma$, then:

Without loss of generality assume $x^{\prime}>x$, then we have two cases:

Case i. If $x>\sigma$, then we have that the Lipschitz constant of $g(x)$ is $L_{\sigma}=|\theta|+|K||(1+\log\sigma)|$ (obtained by taking the supremum of its derivative). Therefore, we have that $|g(x^{\prime})-g(x)|\leq L_{\sigma}\sigma$. Note that $L_{\sigma}\sigma\to 0$ when $\sigma\to 0$ even if $L_{\sigma}\to\infty$, since $L_{\sigma}$ grows logarithmically.

Case ii. If $x\leq\sigma$, then $x^{\prime}\leq x+\sigma\leq 2\sigma$. Therefore:

Now, we have that $-x\log x$ is non-negative for $x\in$. Furthermore, we have that $-x\log x$ is increasing when $x<\exp(-1)$ and decreasing afterwards. First suppose that $2\sigma\leq\exp(-1)$; then $-x^{\prime}\log x^{\prime}\geq-x\log x\geq 0$ which implies:

In the case $2\sigma>\exp(-1)$, then we have:

Combining these two possibilities, we get:

For small $\sigma$, the dominant term of the function $\omega_{g}(\sigma)$ in Lemma E.2 is of the form $C\cdot-\sigma\log\sigma$ for a constant $C$. If we require that this be smaller than some small $\xi>0$, then we can choose an approximate $\sigma$ by solving for $x$ in $-Ax\log x=\xi$ yielding $x=\exp(W_{-1}\frac{\xi}{A})$ where $W_{-1}$ is the negative branch of the Lambert W-function. This is almost linear and yields approximately $x=O(\xi)$ for small $\xi$. In fact, we have that $\omega_{g}(\sigma)\leq C^{\prime}\sigma^{\alpha}$ for any $\alpha<1$, and thus $g$ is “almost” Lipschitz continuous.

That is, for $\vec{\bm{\mu}},\vec{\bm{\mu}}^{\prime}\in\mathcal{M}$ with $\|\vec{\bm{\mu}}^{\prime}-\vec{\bm{\mu}}\|_{\infty}\leq\sigma$, we have:

$\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ can be decomposed into functions of the form $-Kx\log x+\theta x$ (see (15) and (16)) and so we apply the Lemma E.2 element-wise. Let $c$ index the clique component in the marginal vector.

where we recall $n_{c}$ is the number of values that $x_{c}$ can take. By re-using the bound on $\sum_{c}|K_{c}|n_{c}$ from (18), we get the result. ∎

E.3 Bounded Negative Uniform Gap

For the negative TRW objective $f(\vec{\bm{\mu}}):=-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$, the bound $B$ on the negative uniform gap as given in Lemma D.3 for $\bm{u}_{0}$ being the uniform distribution can be taken as:

E.4 Summary

We now give the details of suboptimality guarantees for our suggested algorithm to optimize $f(\vec{\bm{\mu}}):=-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ over $\mathcal{M}$. The (strong) convexity of the negative TRW objective is shown in (Wainwright et al., 2005; London et al., 2015). $\mathcal{M}$ is the convex hull of a finite number of vectors representing assignments to random variables and therefore a compact convex set. The entropy function is continously differentiable on the relative interior of the probability simplex, and thus the TRW objective has the same property on the relative interior of $\mathcal{M}$. Thus $-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ satisfies the properties laid out in Problem 3.2.

For the optimization of $-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ over $\mathcal{M}_{\delta}$ with $\delta\in(0,1]$, the suboptimality is bounded as:

Using $\operatorname{diam}_{\|\cdot\|_{\infty,1}}(\mathcal{M})\leq 2$, and $L_{\delta}$ from Lemma E.1, we can compute $C_{\delta}\leq\operatorname{diam}(\mathcal{M})^{2}L_{\delta}$. Lemma E.3 computes the modulus of continuity $\omega(\sigma)$. The rate then follows directly from Theorem C.1. ∎

Consider the optimization of $-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ over $\mathcal{M}$ with the optimum given by $\vec{\bm{\mu}}^{*}$. The iterates $\vec{\bm{\mu}}^{(k)}$ obtained by running the Frank-Wolfe updates over $\mathcal{M}_{\delta}$ using line-search with $\delta$ updated according to Algorithm 1 (or as summarized in a FCFW variant in Algorithm 4), have suboptimality $h_{k}=\text{TRW}(\vec{\bm{\mu}}^{*};\vec{\bm{\theta}},\bm{\rho})-\text{TRW}(\vec{\bm{\mu}}^{(k)};\vec{\bm{\theta}},\bm{\rho})$ upper bounded as:

$h_{k}\leq\left(\frac{1}{2}\right)^{k}h_{0}+\frac{\widetilde{C}}{\delta_{0}}$ for $k$ such that $h_{k}\geq\max\{B\delta_{0},\frac{2\widetilde{C}}{\delta_{0}}\}$,

$h_{k}\leq\frac{2\widetilde{C}}{\delta_{0}}\left[\frac{1}{\frac{1}{4}(k-k_{0})+1}\right]$ for $k$ such that $B\delta_{0}\leq h_{k}\leq\frac{2\widetilde{C}}{\delta_{0}}$,

$h_{k}\leq\left[\frac{\max(\widetilde{C},B\delta_{0}^{2})B}{\frac{1}{4}(k-k_{1})+1}\right]^{\frac{1}{2}}=O(k^{-\frac{1}{2}})$ for $k$ such that $h_{k}\leq B\delta_{0}$,

$\widetilde{C}:=16|V|\max_{(ij)\in E}(\textrm{VAL}_{i}\textrm{VAL}_{j})$

$k_{0}$ and $k_{1}$ are the number of steps to reach stage 2 and 3 respectively which are bounded as: $k_{0}\leq\max(0,\lceil\log_{\frac{1}{2}}\frac{\widetilde{C}}{h_{0}\delta_{0}}\rceil)$ $k_{0}\leq k_{1}\leq k_{0}+\max\left(0,\lceil\frac{8\widetilde{C}}{B\delta_{0}^{2}}\rceil-4\right)$

Using $\operatorname{diam}_{\|\cdot\|_{\infty,1}}(\mathcal{M})\leq 2$, we bound $\widetilde{C}\leq L\operatorname{diam}_{\|\cdot\|_{\infty,1}}(\mathcal{M})^{2}$ with $L$ (from Property 3.3) derived in Lemma E.1. We bound $-8g_{u}(\vec{\bm{\mu}}^{(k)})$ (the upper bound on the negative uniform gap) using the value derived in Lemma E.4. The rate then follows directly from Theorem D.6 using $p=1$ (see Lemma E.1 where $L_{\delta}\leq\frac{L}{\delta}$). ∎

The dominant term in Lemma E.6 is $\widetilde{C}B\,k^{-\frac{1}{2}}$, with $\widetilde{C}B=O(\|\vec{\bm{\theta}}\|_{1,\infty}|V|)$. We thus find that both bounds depend on norms of $\vec{\bm{\theta}}$. This is unsurprising since large potentials drive the solution of the marginal inference problem away from the centre of $\mathcal{M}$, corresponding to regions of high entropy, and towards the boundary of the polytope (lower entropy). Regions of low entropy correspond to smaller components of the marginal vector, which in turn result in larger and poorly behaved gradients of $-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$, which slows down the resulting optimization.

Appendix F Correction and Local Search Steps in Algorithm 2

Algorithm 5 details the CORRECTION procedure used in line 16 of Algorithm 2 to implement the correction step of the FCFW algorithm. It uses the modified Frank-Wolfe algorithm (FW with away steps), as detailed in Algorithm 3. Algorithm 6 depicts the LOCALSEARCH procedure used in line 17 of Algorithm 2. The local search is performing FW over $\mathcal{M}_{\delta}$ for a fixed $\delta$ using the iterated conditional mode algorithm as an approximate FW oracle. This enables the finding in a cheap of way of more vertices to augment the correction polytope $V$.

Appendix G Comparison to perturbAndMAP

Perturb & MAP. We compared the performance between our method and perturb & MAP for inference on $10$ node Synthetic cliques. We expand on the method we used to evaluate perturbAndMAP in Figure 1(e) and 1(f). We re-implemented the algorithm to estimate the partition function in Python (as described in Hazan and Jaakkola (2012), Section 4.1) and used toulbar2 (Allouche et al., 2010) to perform MAP inference over an inflated graph where every variable maps to five new variables. The log partition function is estimated as the mean energy of 10 exact MAP calls on the expanded graph where the single node potentials are perturbed by draws from the Gumbel distribution. To extract marginals, we fix the value of a variable to every assignment, estimate the log partition function of the conditioned graph and compute beliefs based on averaging the results of adding the unary potentials to the conditioned values of the log partition function.

Appendix H Correction Steps for Frank-Wolfe over ℳℳ\mathcal{M}

Recall that the correction step is done over the correction polytope, the set of all vertices of $\mathcal{M}$ encountered thus far in the algorithm. On experiments conducted over $\mathcal{M}$, we found that using a better correction algorithm often hurt performance. This potentially arises in other constrained optimization problems where the gradients are unbounded at the boundaries of the polytope. We found that better correction steps over the correction polytope (the convex hull of the vertices explored by the MAP solver, denoted $V$ in Algorithm 2), often resulted in a solution at or near a boundary of the marginal polytope (shared with the correction polytope). This resulted in the iterates becoming too small. We know that the Hessian of $\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ is ill conditioned near the boundaries of the marginal polytope. Therefore, we hypothesize that this is because the gradient directions obtained when the iterates became too small are simply less informative. Consequently, the optimization over $\mathcal{M}$ suffered. We found that the duality gap over $\mathcal{M}$ would often increase after a correction step when this phenomenon occurred. The variant of our algorithm based on $\mathcal{M}_{\delta}$ is less sensitive to this issue since the restriction of the polytope bounds the smallest marginal and therefore also controls the quality of the gradients obtained.

Appendix I Additional Experiments

Figure 5(a), 5(b) depicts the comparison of convergence of algorithm variants over $\mathcal{M}$ and $\mathcal{M}_{\delta}$ (same setup as Figure 1(a), 1(b). Here, we plot $\zeta_{\mu}$.

Appendix J Bounding log⁡Z𝑍\log Z with Approximate MAP Solvers

Suppose that we use an approximate MAP solver for line 7 of Algorithm 2. We show in this section that if the solver returns an upper bound on the value of the MAP assignment (as do branch-and-cut solvers for integer linear programs), we can use this to get an upper bound on $\log Z$. For notational consistency, we consider using Algorithm 2 for $\min_{\bm{x}\in\mathcal{D}}f(\bm{x})$, where $f(\bm{x})=-\text{TRW}(\vec{\bm{\mu}};\vec{\bm{\theta}},\bm{\rho})$ is convex, $\bm{x}=\vec{\bm{\mu}}$, and $\mathcal{D}=\mathcal{M}$.

The property that the duality gap may be used as a certificate of optimality (Jaggi, 2013) gives us:

Adding the gap onto the TRW objective yields an upper bound on the optimum (which from Equation 1 is an upper bound on $\log Z$), i.e. $\log Z\leq-f(\bm{x}^{*})$. From our definition of the duality gap $g(\bm{x}^{(k)})$ (line 8 in Algorithm 2) and (24), we have:

where $\bm{s}^{(k)}=\operatorname*{\arg\min}_{\bm{v}\in\mathcal{D}}\left\langle\nabla f(\bm{x}^{(k)}),\bm{v}\right\rangle=\operatorname*{\arg\max}_{\bm{v}\in\mathcal{D}}\left\langle-\nabla f(\bm{x}^{(k)}),\bm{v}\right\rangle$ (line 7 in Algorithm 2). Thus, if the approximate MAP solver returns an upper bound $\kappa$ such that $\max_{\bm{v}\in\mathcal{D}}\left\langle-\nabla f(\bm{x}^{(k)}),\bm{v}\right\rangle\leq\kappa$, then we get the following upper bound on the log-partition function: