Self-Guided Belief Propagation -- A Homotopy Continuation Method

TL;DR

Self-guided belief propagation (SBP) introduces a homotopy continuation approach that gradually incorporates pairwise potentials, leading to more accurate and stable probabilistic inference solutions in graphical models, especially when traditional BP fails to converge.

Contribution

This work presents SBP, a novel homotopy continuation method for belief propagation that guarantees convergence to a global optimum in certain models and improves accuracy in general.

Findings

SBP outperforms BP in accuracy when BP converges.

SBP finds unique, stable solutions when BP does not converge.

SBP converges to the global optimum for attractive models.

Abstract

Belief propagation (BP) is a popular method for performing probabilistic inference on graphical models. In this work, we enhance BP and propose self-guided belief propagation (SBP) that incorporates the pairwise potentials only gradually. This homotopy continuation method converges to a unique solution and increases the accuracy without increasing the computational burden. We provide a formal analysis to demonstrate that SBP finds the global optimum of the Bethe approximation for attractive models where all variables favor the same state. Moreover, we apply SBP to various graphs with random potentials and empirically show that: (i) SBP is superior in terms of accuracy whenever BP converges, and (ii) SBP obtains a unique, stable, and accurate solution whenever BP does not converge.

Figures

Tables

Table 1. TABLE I: Results for general models with J i ​ j ∈ { − 1 , 1 } subscript 𝐽 𝑖 𝑗 1 1 J_{ij}\in\{-1,1\} on grid graphs ( N = 25 𝑁 25 N\!=\!25 and N = 100 𝑁 100 N\!=\!100 ), complete graphs ( N = 10 𝑁 10 N\!=\!10 ), and random graphs ( N = 10 𝑁 10 N\!=\!10 ). We report the MSE to the exact marginals and the MSE B subscript MSE 𝐵 \text{MSE}_{{B}} to the Bethe approximation, convergence ratio, and the overall number of BP iterations. Only converged runs are considered for BP ∘ superscript BP \text{BP}^{\circ} and BP D ∘ superscript subscript BP D \text{BP}_{\text{D}}^{\circ} but all runs are considered for SBP a ​ l ​ l subscript SBP 𝑎 𝑙 𝑙 \text{SBP}_{all} , Gibbs a ​ l ​ l subscript Gibbs 𝑎 𝑙 𝑙 \text{Gibbs}_{all} , and ℱ B , a ​ l ​ l ∗ superscript subscript ℱ 𝐵 𝑎 𝑙 𝑙 \mathcal{F}_{B,all}^{*} .

		Grid Graph( $5\times 5$)			Grid Graph ($10\times 10$)			Complete Graph			Random Graph
	$\theta$	$0$	$0.1$	$0.4$	$0$	$0.1$	$0.4$	$0$	$0.1$	$0.4$	$0$	$0.1$	$0.4$
MSE	${\text{BP}}^{\circ }$	0.338	0.251	0.102	-	-	0.184	0.463	0.466	0.356	0.252	0.202	0.101
	${\text{BP}}_{\text{D}}^{\circ }$	0.226	0.198	0.066	0.186	0.240	0.154	0.463	0.473	0.422	0.128	0.116	0.083
	${\text{SBP}}_{all}$	0.000	0.029	0.047	0.000	0.026	0.077	0.000	0.055	0.074	0.000	0.048	0.049
	${ℱ}_{B,all}^{\ast }$	0.036	0.042	0.069	-	-	-	-	-	-	-	-	-
	${\text{Gibbs}}_{all}$	0.001	0.016	0.064	0.001	0.037	0.120	0.096	0.096	0.077	0.001	0.011	0.048
Convergence ratio	${\text{BP}}^{\circ }$	0.05	0.11	0.26	0.00	0.00	0.02	0.41	0.42	0.50	0.30	0.33	0.49
	${\text{BP}}_{\text{D}}^{\circ }$	0.11	0.16	0.69	0.01	0.02	0.12	0.41	0.41	0.50	0.62	0.64	0.80
Number of iterations	${\text{BP}}^{\circ }$	40	52	84	-	-	102	17	17	18	42	53	50
	${\text{BP}}_{\text{D}}^{\circ }$	1370	1449	1735	2711	2313	2599	211	207	234	1077	1057	873
	${\text{SBP}}_{all}$	5	182	146	5	149	209	5	51	110	5	149	131
	${\text{Gibbs}}_{all}$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$
${\text{MSE}}_B$	${\text{SBP}}_{all}$	0.036	0.037	0.022	-	-	-	-	-	-	-	-	-
${ℱ}_B^{\circ }\left({\zeta }_M\right)$ equals	SBP	100	10	23	-	-	-	-	-	-	-	-	-

Equations

$$P_{\mathbf{X}}(\mathbf{x})$$

$$\displaystyle=\frac{1}{\mathcal{Z}}\exp\big{(}-E(\mathbf{x})\big{)},$$

$$P_{\mathbf{Y}}(\mathbf{y})=\sum _{x_i:X_i\in \{\mathbf{X}\backslash \mathbf{Y}\}}{P}_{\mathbf{X}}(\mathbf{x}),$$

$$m_i=E(X_i)=P_{X_i}(+1)-P_{X_i}(-1),$$

$${\chi }_{ij}=E(X_i,X_j)=\sum _{x_i,x_j}{x}_i{x}_{j}P(x_i,x_j).$$

$$P_{\mathbf{X}}(\mathbf{x})$$

$${\mu }_{ij}^{n+1}(x_j)\propto \sum _{x_i\in \mathcal{X}}\Phi (x_i,x_j)\Phi (x_i)\prod _{X_k\in N(i)\backslash j}{\mu }_{ki}^n(x_i)$$

$${\mathbf{\mu }}^{n+1}$$

$$=(1-ϵ)\mathcal{B}\mathcal{P}({\mathbf{\mu }}^n)+ϵ{\mathbf{\mu }}^n.$$

$$\tilde{P}(x_i)=$$

$$\tilde{P}(x_i,x_j)=$$

$$\prod _{X_k\in N(i)\backslash j}{\mu }_{ki}^{\circ }(x_i)\cdot \prod _{X_l\in N(j)\backslash i}{\mu }_{lj}^{\circ }(x_j),$$

$${\tilde{P}}_B=\{\tilde{P}(x_i),\tilde{P}(x_i,x_j):X_i\in \mathbf{X},(i,j)\in \mathbf{E}\}.$$

$$E_B($$

$$-\sum _{(i,j)\in \mathbf{E}}\sum _{x_i,x_j}\tilde{P}(x_i,x_j)\ln \Phi (x_i,x_j).$$

$$S_B($$

$$+\sum _{X_i}(|N(i)|-1)\sum _{x_i}\tilde{P}(x_i)\ln \tilde{P}(x_i).$$

$${\mathcal{F}}_B({\tilde{P}}_B)=$$

$$\mathcal{L}=\{{\tilde{P}}_B:\sum _{x_i}\tilde{P}(x_i)=1,\sum _{x_j}\tilde{P}(x_i,x_j)=\tilde{P}(x_i)\}.$$

$${\mathcal{F}}_B^{\circ }={\mathcal{F}}_B({\tilde{P}}_B^{\circ }).$$

$$P(x_i)=\frac{1}{\sum _{m=1}^M{{\mathcal{Z}}_{\mathcal{B}}}^m}\sum _{m=1}^M{\tilde{P}}^m(x_i){{\mathcal{Z}}_{\mathcal{B}}}^m.$$

$${\Phi }_m(x_i,x_j)$$

$$=\Phi (x_i,x_j{)}^{{\zeta }_m}.$$

$${\mathbf{\mu }}_M^{\circ }={\mathcal{B}\mathcal{P}}_M^{\circ }\left({\mathcal{B}\mathcal{P}}_{M-1}^{\circ }\left(\cdots {\mathcal{B}\mathcal{P}}_1^{\circ }\left({\mathbf{\mu }}_1^{\text{init}}\right)\right)\right).$$

$$H(\mathbf{\mu },\zeta )=\mathbf{\mu }-\mathcal{B}\mathcal{P}(\mathbf{\mu })\text{where}\Psi =\Psi (\zeta ).$$

$$c(\zeta ):H(\mathbf{\mu },\zeta )=0$$

$${\mu }_{ij}^{n+1}(x_j)\propto \sum _{x_i\in \mathcal{X}}\Phi (x_i)\prod _{X_k\in N(i)\backslash j}{\mu }_{ki}^n(x_i).$$

$$\tilde{P}(x_i)$$

$$=\frac{e^{x_i{\theta }_i}}{e^{{\theta }_i}+e^{-{\theta }_i}}.$$

$$P(x_i)=\frac{1}{\mathcal{Z}}{e}^{x_i{\theta }_i}\sum _{x_j:X_j\in \{\mathbf{X}\backslash X_i\}}\prod _{X_j\in \{\mathbf{X}\backslash X_i\}}{e}^{x_j{\theta }_j}=\frac{e^{x_i{\theta }_i}}{e^{{\theta }_i}+e^{-{\theta }_i}}.$$

$$\frac{\partial {\mathcal{F}}_B(\zeta )}{\partial \zeta }$$

$$=-\frac{\partial }{\partial \zeta }\sum _{(i,j)\in \mathbf{E}}\sum _{x_i,x_j}\tilde{P}(x_i,x_j)\cdot \zeta J_{ij}\cdot x_i{x}_j$$

$$=-\sum _{(i,j)\in \mathbf{E}}{J}_{ij}\cdot {\chi }_{ij}.$$

$$\underset{\zeta \to {\zeta }_c^{+}}{\lim }{m}_i^s=\underset{\zeta \to {\zeta }_c^{+}}{\lim }{m}_i^t.$$

$$\underset{\zeta \to {\zeta }_c^{+}}{\lim }{m}_i^s\ne \underset{\zeta \to {\zeta }_c^{+}}{\lim }{m}_i^t.$$

$${\mu }_{ij}(X_j=+1)\ge {\mu }_{ij}(X_j=-1).$$

Full text

1 Introduction

Computing the marginal distributions and evaluating the partition function are two fundamental problems in the context of probabilistic graphical models. Both problems can be solved efficiently for tree-structured models but are NP-hard if the graphical model contains loops [7]. Since loops are present in many problems of relevance there is a need for efficient approximation methods.

Belief propagation (BP) provides a way to approximate the marginal distribution and may be viewed as an approach to try to minimize the Bethe free energy $\mathcal{F}_{B}$ [61]. BP has a long success story in many applications, including computer vision, speech processing, social network analysis, and error-correcting codes [25, 19, 35]. But, despite its success, BP does not always approximate the marginals well and may even fail to converge. Reasons for BP’s failure include the existence of multiple fixed points of varying accuracy (where it depends on implementation details to which one BP converges) and of unstable fixed points that cause BP to oscillate far away from any fixed point [51, 33, 16]. These limitations motivate the search for modifications of BP that have better convergence properties and find more accurate marginals.

One alternative is to perform a different optimization problem that minimizes the Bethe free energy $\mathcal{F}_{B}$, instead of applying BP [61]. While, for certain problem classes, polynomial-time algorithms exist, it is even problematic to approximate the global minimum for models with arbitrary potentials [5, 42, 55]. Hence, the pursuit for methods that approximate the marginals well – ideally with both run-time and convergence guarantees – is still ongoing.

In this work, we present self-guided belief propagation (SBP) to help address this gap. The observation that strong pairwise potentials often reduce accuracy and worsen the convergence properties [21] inspired us to construct a homotopy; i.e., we first consider only local potentials (where BP is exact) and subsequently modify the model by gradually increasing the pairwise potentials to their given values. SBP thus solves a deterministic sequence of models that iteratively refines the Bethe approximation towards a solution that is uniquely defined by the initial model.

We empirically evaluate SBP for grid-graphs, complete graphs, and random graphs with random potentials. Compared to BP, we observe superior performance in terms of accuracy; in fact, SBP achieves more accurate results than Gibbs sampling in a fraction of run-time. Additionally, SBP exhibits favorable convergence properties and excels for hard problems, for which – despite BP failing to converge – SBP provides accurate results. We also analyze SBP theoretically and demonstrate optimality of the selected fixed point for attractive models with unidirectional local potentials.111An attractive model is a probabilistic graphical model where all pairwise interactions are specified by positive couplings; a general model also contains pairwise interactions with negative couplings. Unidirectional local potentials are potentials that all favor the same state (see Section 2.2).

It is worth emphasizing that SBP is a relatively straightforward modification to standard BP that can be easily added to existing BP implementations. We thus expect that the ease of use lowers the hurdle for practical applications.

The paper is structured as follows: Section 2 provides background on probabilistic graphical models, belief propagation, and methods that minimize the free energy. Moreover, we present a unified taxonomy of model-classes according to their complexity. Our proposed algorithm is presented in Section 3. We evaluate SBP and discuss empirical observations in Section 4 and provide a more formal analysis in Section 5. Finally, we conclude the paper in Section 6.

2 Background

In this section, we briefly introduce probabilistic graphical models and specify the models considered in this work. We further introduce the BP algorithm and its connection to the Bethe approximation. Finally, we summarize important properties regarding the solution space of BP as a foundation for our theoretical analysis of SBP.

2.1 Probabilistic Graphical Models

Let us consider an undirected graph $\mathcal{G}=(\mathbf{X},\mathbf{E})$, where $\mathbf{X}=\{X_{1},\dots,X_{N}\}$ is the set of nodes, and $\mathbf{E}$ is the set of undirected edges. Then, two nodes $X_{i}$ and $X_{j}$ are joined by an edge if $(i,j)\in\mathbf{E}$; note that we consider each edge only once, i.e., $(i,j)=(j,i)$. We further denote the set of neighbors of $X_{i}$ by $N({i})=\{X_{j}\in\mathbf{X}:(i,j)\in\mathbf{E}\}$.

Let us define a probabilistic graphical model $\mathcal{U}=(\mathcal{G},\Psi)$, where $\Psi=\{\Phi({\mathbf{x}_{1}}),\ldots,\Phi({\mathbf{x}_{K}})\}$ is the set of all $K$ potentials and $\mathbf{X}$ is the set of random variables. In this work we focus on pairwise models where all potentials consist of two variables at most, so that the joint distribution $P_{\mathbf{X}}(\mathbf{x})$ factorizes according to

[TABLE]

where the potentials specify the energy $E(\mathbf{x})$.

We consider the following two problems: First, evaluating the partition function $\mathcal{Z}$ (i.e., the normalization function of the joint distribution). Consider any distribution $Q_{\mathbf{X}}(\mathbf{x})$ and the (Gibbs) free energy \mathcal{F}=\mathcal{F}(Q_{\mathbf{X}}(\mathbf{x}))=\sum_{\mathbf{x}}Q_{\mathbf{X}}(\mathbf{x})\big{(}E(\mathbf{x})+\ln Q_{\mathbf{X}}(\mathbf{x})\big{)}; then, evaluating the partition function and minimizing the free energy are equivalent since the minimum of $\mathcal{F}$ corresponds to $\mathcal{Z}$ according to $\min\mathcal{F}=\mathcal{F}^{*}=-\ln\mathcal{Z}$ (cf. [61]).

Second, we consider computing the marginal distribution

[TABLE]

where $\mathbf{Y}\subset\mathbf{X}$ may be any set of RVs. The identities of the random variables are often obvious from the notation of the values; therefore, we often use the shorthand notation $P(\mathbf{y})$ for the marginal distribution. Note that the problems of evaluating the partition function and computing the marginal distributions are inherently linked, as $\mathcal{F}$ obtains its minimum precisely for the marginal distributions $P_{\mathbf{X}}(\mathbf{x})$.

In general, both problems are intractable [7]. Consequently one often relaxes the problems and only approximates the marginal distributions and the partition function. This allows for an elegant iterative algorithm that was discovered multiple times in different fields. It is known as belief propagation (BP) in computer science, the sum-product algorithm in information theory, and the cavity- or the Bethe-method in physics; we refer the reader to [26, 30] for a good overview of the underlying principles.

2.2 MODEL SPECIFICATION

We consider binary pairwise models where every random variable $X_{i}$ takes values from $\mathcal{X}=\{-1,+1\}$ (i.e., Ising models). For these models, instead of considering the singleton marginals $\tilde{P}(x_{i})$ and the pairwise marginals $\tilde{P}(x_{i},x_{j})$ explicitly, it is often more convenient to work with the means $m_{i}$ and the correlations $\chi_{ij}$ defined according to

[TABLE]

[TABLE]

Let us define couplings $J_{ij}\in\mathbb{R}$ that are assigned to each edge $(i,j)\in\mathbf{E}$ and local fields $\theta_{i}\in\mathbb{R}$ that act on each variable $X_{i}\in\mathbf{X}$. These parameters then define the pairwise potentials $\Phi(x_{i},x_{j})=\exp(J_{ij}x_{i}x_{j})$ and the local potentials $\Phi(x_{i})=\exp(\theta_{i}x_{i})$ so that the corresponding joint distribution from (1) can be written as

[TABLE]

There are two different types of interactions between random variables: if $J_{ij}$ is negative then the edge $(i,j)$ is repulsive; if $J_{ij}$ is positive then the edge $(i,j)$ is attractive. Accordingly, one distinguishes between general models that contain both attractive and repulsive edges and attractive models222 Note that attractive models are also known as ferromagnetic models [30] or log-supermodular models [40]. that contain only attractive edges.

Attractive models are often more well-behaved than general models but, rather than only distinguishing between these two model-types, we find it beneficial to consider an even finer-grained distinction into the following three model-classes. In doing so, we additionally unify different naming conventions in the literature and provide a consistent taxonomy of model-classes.

2.2.1 Frustrated, Balanced, and Unidirectional Models

First, frustrated models are models that contain cycles for which the product over all couplings equates to a negative number, i.e., cycles with an odd number of repulsive edges (see Fig. 1a).333 Frustrated models are also referred to as spin glasses in the physics literature. These are some of the most complex models considered of the form (5). Note that if a model is frustrated, it must be a general model per definition.

Second, balanced models are models that are not frustrated. For a balanced model, it is possible to flip a certain subset of variables without affecting the represented distribution such that the resulting model is attractive (cf. [52]). To flip a variable $X_{i}$, one must switch its two states, and: reverse the sign of the corresponding local field $\theta_{i}$; and reverse the signs of the couplings $J_{ij}$ for all incident edges $(i,j):X_{j}\in N({i})$. For an illustrative example, consider the model in Fig. 1b. This model is a general one in its current form; reversing the sign of the rightmost variable further implies changing both attached edges from repulsive to attractive, which results in an equivalent attractive model (see Fig. 1c). Such an equivalent model exists if and only if the original model contains no cycles with an odd number of repulsive edges, i.e., if the original model was not frustrated [14].

Finally, unidirectional models are attractive models with unidirectional local potentials, i.e., either all $\theta_{i}\leq 0$ or $\theta_{i}\geq 0$ (see Fig. 1d).444 Unidirectional models are also referred to as attractive models in an external field, true ferromagnetic models, or consistent models in the literature. In a similar manner as for balanced models, certain general models have equivalent unidirectional models, i.e., by flipping variables to render the model attractive all local potentials are mirrored so that the attractive model is unidirectional; we call these models quasi-unidirectional. All theoretical findings for unidirectional models thus readily extend to quasi-unidirectional models.

The major advantage of considering unidirectional models is that they are relatively well understood. Thus unidirectional models are of particular relevance for our work – as well as for studying inference algorithms in general.

Intuitively, in unidirectional models, we expect that all variables will energetically favor the same state as indicated by the local potentials. Additionally, since all edges are attractive, no contradictions between neighboring variables will occur, thus further reinforcing this behavior. Moreover, various research directions indicate the relevance of this model class. We will list a couple of these indications below:

On the one hand, it is generally often easier to study homogeneous models, where all pairwise potentials are the same (i.e., $J_{ij}=J$) and all local potentials are the same (i.e., $\theta_{i}=\theta$). Provided that the model is unidirectional, the findings from homogeneous models carry over readily [21, 12]. Alternatively, one often considers models of the form (5) in an external field $\theta$, i.e, where all local potentials are the same (or have the same sign) [17]. Keeping the physical interpretation in mind, we might hope that models that are placed in an external field are easier to understand. Let us use the concept of flipping variables for transforming balanced models into attractive models again. Interestingly, if we consider models in an external field only, then flipping variables is not a viable option as this would imply reversing the sign of some $\theta_{i}$. Thus for models in a constant external field, the class of attractive models reduces precisely to the class of unidirectional models.

On the other hand, the complexity of approximating the partition function also depends on the model class and is in line with our distinction. That is for unidirectional models the partition function can be approximated efficiently in polynomial time [17]; for frustrated models, the partition function cannot be approximated efficiently; balanced models are of intermediate complexity  [12].

Finally, it is sometimes beneficial to study models without local potentials; in fact, one can find an auxiliary model without local potentials for every model of the form (5) by adding one variable to the graph (and connecting it to all existing variables), where the auxiliary model has the same marginals and the same partition function (up to a constant factor of two) [53, 18, 41]. Provided that the original model was attractive, then the auxiliary model remains attractive (or can be rendered attractive by flipping the additional variable) if and only if the original model was unidirectional.555 On the one hand, we need to allow for a flip of the additional variable for the case that all $\theta_{i}<0$; on the other hand, this makes the statement agnostic to the definition (i.e., it holds irrespective if we have $\mathcal{X}=\{-1,+1\}$ or $\mathcal{X}=\{0,1\}$).

2.3 BELIEF PROPAGATION (BP)

BP approximates the marginals by recursively exchanging messages between random variables. The messages $\mu^{n+1}_{ij}(x_{j})$ from $X_{i}$ to $X_{j}$ at iteration $n+1$ are updated according to (6) and are normalized so that $\sum_{x_{j}\in\mathcal{X}}\mu^{n}_{ij}(x_{j})=1$.

[TABLE]

Let $\bm{\mu}^{n}=\{\mu^{n}_{ij}(x_{j}):e_{ij}\in\mathbf{E},x_{j}\in\mathcal{X}\}$ be the set of all messages at iteration $n$ and let the mapping induced by (6) be denoted as $\bm{\mu}^{n+1}=\mathcal{BP}(\bm{\mu}^{n})$. If all successive messages remain unchanged, i.e., if $\bm{\mu}^{n+1}=\bm{\mu}^{n}$, then BP is converged to a fixed point $\bm{\mu}^{\circ}$. We further write $\bm{\mu}^{\circ}=\mathcal{BP}^{\circ}(\bm{\mu}^{\text{init}})$, where $\mathcal{BP}^{\circ}$ updates the messages until convergence. If BP fails to converge and the messages oscillate, one can try to achieve convergence by either changing the update-rule [10, 47, 23], or by replacing the messages with a convex combination of the last messages [34]. The latter method is known as damping ($\text{BP}_{\text{D}}$) where a damping parameter $\epsilon\in[0,1)$ specifies the new update rule

[TABLE]

After convergence, the singleton marginals $\tilde{P}(x_{i})$ and pairwise marginals $\tilde{P}(x_{i},x_{j})$ are approximated by

[TABLE]

where $Z_{i},Z_{ij}\in\mathbb{R}_{+}^{*}$ guarantee that all probabilities sum to one. We refer to the set of singleton and pairwise marginals as pseudomarginals ${\tilde{P}_{B}}$, where

[TABLE]

We further denote the pseudomarginals obtained at a fixed point of BP by ${\tilde{P}_{B}^{\circ}}$.666 With slight abuse of notation, we will refer to both the pseudomarginals ${\tilde{P}_{B}^{\circ}}$ and the messages $\bm{\mu}^{\circ}$ as fixed points throughout this work. Strictly speaking, however, the pseudomarginals are only evaluated at a fixed point rather than being a fixed point themselves.

2.4 THE BETHE APPROXIMATION & RELATED WORK

BP is closely connected to variational methods and concepts from statistical physics. In particular, a direct relationship between performing BP and minimizing the Bethe free energy exists (cf. [61]).

The Bethe free energy $\mathcal{F}_{B}({\tilde{P}_{B}})$ is evaluated over the pseudomarginals (i.e., the singleton- and pairwise marginals) and is a function of the average energy $E_{{B}}({\tilde{P}_{B}})$ and the Bethe entropy $S_{{B}}({\tilde{P}_{B}})$ that are defined by

[TABLE]

This subsequently defines the Bethe free energy according to

[TABLE]

Let us define the local polytope $\mathcal{L}$, i.e, the set of consistent pseudomarginals, where

[TABLE]

Minimizing $\mathcal{F}_{B}$ over the local polytope gives us the desired global minimum $\mathcal{F}_{B}^{*}=\min_{\mathcal{L}}\mathcal{F}_{B}({\tilde{P}_{B}})$. This, however, is not straightforward as the constrained $\mathcal{F}_{B}$ is often non-convex.

The main reason for considering the Bethe free energy is its immediate connection with the fixed points of BP $\bm{\mu}^{\circ}$ (and the associated pseudomarginals ${\tilde{P}_{B}^{\circ}}$): all stationary points $\mathcal{F}_{B}^{\circ}$ are in a one-to-one correspondence with the fixed points of BP. More precisely, we have

[TABLE]

In practice, we are particularly interested in stable fixed points (for which BP converges if initialized close enough); we index all stable fixed points ${\tilde{P}_{B}^{s}}$ by $s=1,\ldots,S$ and denote the set of all stable fixed points by $\bm{S}$. Likewise, we consider the set of all fixed points that constitute minima of the Bethe free energy $\bm{M}$ and index them by $m=1,\ldots,M$. Note that every stable fixed point ${\tilde{P}_{B}^{s}}$ must be a minimum of $\mathcal{F}_{B}$; the converse, however, need not be the case, i.e., not every local minimum corresponds to a stable fixed point, so that $\bm{S}\subseteq\bm{M}$ [50].

This correspondence between BP and $\mathcal{F}_{B}$ led to a better understanding of BP and inspired plenty of methods that minimize $\mathcal{F}_{B}$ directly [58, 62]. The minimization, however, is still highly non-trivial and requires good approximation methods in practice. Since strong pairwise potentials often reduce the accuracy of the Bethe approximation and are responsible for its non-convexity [20], one can correct the entropy term (12) by accounting for the strong potentials; this admits convex relaxations that provide provable convergent message passing algorithms [49, 28, 11, 15, 29]. There is, however, a trade-off between convergence-properties and accuracy in general. That is, if it can be minimized, the Bethe approximation often provides accurate results and outperforms convex free energy approximations [29, 55]. Thus, it is a relevant problem to directly approximate $\mathcal{F}_{B}$ in a way that allows for efficient minimization. Polynomial run-time algorithms exist that approximate $\mathcal{F}_{B}$ for restricted models: these include sparsity constraints [42] or require attractive models [55]. If both properties are fulfilled, i.e., for locally tree-like attractive models the Bethe approximation is exact and can be optimized efficiently [9]. Note that $\mathcal{F}_{B}$ provides an upper bound on $\mathcal{F}$ for attractive models [40, 60, 56].

We aim to estimate $\mathcal{F}_{B}$ and ${P_{B}}$; while the approximation of $\mathcal{F}_{B}$ in [55] is $\varepsilon$-accurate, no run-time guarantees exist for general models. Our work, on the contrary, provides an estimate in constant run-time (see Theorem 6 in Section 5); the approximation error, however, cannot be made arbitrarily small in general. Both methods overcome their respective disadvantages when restricting the models; i.e., both methods do efficiently minimize the Bethe approximation for unidirectional models.777 The restriction to unidirectional models is actually too conservative for [55], where the approximation of $\mathcal{F}_{B}$ takes polynomial run-time for balanced models as well.

2.5 Solution Space of Binary Pairwise Models

Over the years, a substantial amount of work has been devoted to getting a good understanding of how the solution space is affected by the specification of the graphical model. We are particularly interested in how the performance of BP and the energy landscape of $\mathcal{F}_{B}$ vary, depending on the specification of the potentials. Here, we briefly summarize the most important results with a particular focus on how the solution space changes when considering increasingly stronger couplings.

This discussion serves two purposes: on the one hand, it provides a good intuition and motivates the working principle of SBP discussed in Section 3; on the other hand, it prepares the required background for our theoretical analysis in Section 5. Moreover, the difference in the behavior of BP is often only considered between general and attractive models. Instead, we consider the difference in BP’s behavior between frustrated, balanced, and unidirectional models. We are convinced that this distinction is not only relevant for analyzing SBP in this work but that this is a crucial distinction one must make in the pursuit of understanding what makes some problems inherently hard to solve.

For as long as the couplings $J_{ij}$ are sufficiently small (relative to $\theta_{i}$), the Bethe free energy $\mathcal{F}_{B}$ is convex, BP has a unique fixed point, and BP converges [20, 33]. As the couplings increase, the performance tends to suffer as multiple fixed points emerge and as BP may fail to converge (even if the fixed point is still unique) [33, 21].888 Note that the differences in behavior coincide with different phases in statistical mechanics (cf. [64, 30, 48]). The precise effect of the strong couplings on the behavior of BP and the shape of $\mathcal{F}_{B}$ mainly depends on the respective model-class. Thus, we will split the following discussion accordingly.

For frustrated models, the number of stationary points scales exponentially with the number of variables. Often, multiple stationary points coexist with the same value of $\mathcal{F}_{B}$ but possibly completely different marginals; it may even happen that there is not a single unambiguous global minimum. Moreover, for frustrated models, the Bethe free energy is not lower bounded by the exact free energy, so that the global minimum need not approximate the free energy best. Such complex energy landscapes render any attempt of minimizing the Bethe free energy futile. Moreover, in particular for models with strong couplings, BP often fails to converge for frustrated models.

For balanced models, the number of stationary points still scales exponentially with the number of variables. For example, attractive models with arbitrary local potentials (i.e., random field Ising models) are some of the simplest forms of disordered systems [2]. Yet, in comparison with frustrated models, they are much easier to solve. Since every balanced model can be transformed into an equivalent attractive model, the Bethe free energy is lower bounded by the exact free energy. Moreover, many methods that minimize $\mathcal{F}_{B}$ are particularly efficient for attractive models and thus for balanced models as well. One of the properties that we will use throughout is that for attractive models all fixed points corresponding to local minima of $\mathcal{F}_{B}$ are stable [50, Theorem 6], [36, 46]. Although BP may still fail to converge on balanced models, it is thus at least guaranteed that BP will converge if the messages are initialized close enough to a fixed point.

For unidirectional models two fixed points exist at most. This not only renders approximating the free energy and the marginals relatively straightforward but also allows for a thorough evaluation of various algorithms. Since we devote a great part of our theoretical analysis of SBP to its performance on unidirectional models, we will now discuss the solution space of unidirectional models in more detail and prepare some foundations for the later analysis.

Most properties result from generalizing the properties of homogeneous models. As discussed in Section 2.2.1 these models qualitatively behave the same.

First, we treat the special case where all local potentials are specified by $\theta_{i}=0$ and where both fixed points – if they exist – are symmetric.

Lemma 1**.**

Let us consider a unidirectional model with loops, with $\theta_{i}=0$, and with couplings strong enough to admit multiple (i.e., two) fixed points. Then, both fixed points $s$ and $t$ are symmetric, i.e., $m_{i}^{s}=-m_{i}^{t}$, and both fixed points have the same Bethe free energy, i.e., $\mathcal{F}_{B}^{s}=\mathcal{F}_{B}^{t}$.

Proof:

Given that the local potentials are the same for both states, the message update rule from (6) simplifies to $\mu^{n+1}_{ij}(x_{j})\propto\sum_{x_{i}}\Phi(x_{i},x_{j})\prod_{X_{k}}\mu^{n}_{ki}(x_{i})$. Consequently, if $\mu^{s}_{ij}(x_{j})=\mu^{\circ}_{ij}(x_{j})$ are fixed point messages, so are the symmetric messages $\mu^{t}_{ij}(x_{j})=1-\mu^{\circ}_{ij}(x_{j})$. By computing the pseudomarginals according to (8)-(9) it follows immediately that the means are symmetric around 0 and that both fixed points have the same Bethe free energy (see Example 1). ∎

We will now present an illustrative example that describes the behavior of models with $\theta=0$. These models have been well-studied in the literature [32, 57, 20]. We present Example 1 to summarize the most important properties that serve as the basis for our theoretical analysis of SBP in Section 5.

Example 1**.**

Consider a homogeneous model on a complete graph with $N=4$ variables, positive couplings $J_{ij}=J>0$, and local potentials specified by $\theta_{i}=0$. For computing the Bethe free energy $\mathcal{F}_{B}$ according to (13), we enforce the symmetry of the model by setting all singleton marginals to the same value $\tilde{P}(x_{i})={\tilde{P}}$. We compute the Bethe-optimal pairwise marginals according to [59, Lemma 1] (see also [57]). This allows us to compute the pairwise marginals given only the singleton marginals.999 Note that an alternative representation of (5) is used in [59, 57]; there is, however, a straightforward mapping from one representation to the other.

The effect of the coupling strength is illustrated in Fig. 2. Note that for $\theta=0$, ${\tilde{P}}=(0.5,0.5)$ is a fixed point for every value of $J$; we will often refer to this fixed point in terms of its means $m_{i}=0$.101010 This fixed point exists whenever all local potentials are specified by $\theta=0$ and can be computed analytically [32]. Consequently, there is also a corresponding stationary point of $\mathcal{F}_{B}$. For small values of $J$, $\mathcal{F}_{B}$ is convex and this fixed point is unique (see Fig. 3a); as $J$ increases, $\mathcal{F}_{B}$ becomes flat (see Fig. 3b), before turning into a local maximum with two symmetric minima branching off (see Fig. 3c).

Second, when all local potentials are specified by $\theta_{i}>0$ the fixed points are not symmetric anymore. For later reference, we summarize some important properties of the fixed points in the following Lemma.

Lemma 2**.**

Let us consider a unidirectional model with $\theta_{i}>0$ and couplings strong enough to admit multiple (i.e., two) fixed points. Then, one fixed point has its means aligned with the local potentials, i.e., $m_{i}^{s}>0$; we call this fixed point consistent or state-preserving (cf. [12, 22]). Similar, as in Lemma 1, the means of the second fixed point have opposing signs, i.e., $m_{i}^{t}<0$. Both fixed points are not symmetric, however, and the consistent fixed point has larger means $|m_{i}^{s}|>|m_{i}^{t}|$ and larger correlations $\chi_{ij}^{s}>\chi_{ij}^{t}$. Accordingly, the consistent fixed point constitutes the global minimum of the Bethe free energy, i.e., $\mathcal{F}_{B}^{*}=\mathcal{F}_{B}^{s}$.

Proof:

First, the existence of a consistent fixed point, i.e., $m_{i}^{s}>0$, is an immediate consequence of bounding the location of the fixed points of BP in [54, Theorem 3],[61]. Second, the existence of a second point with opposing means, i.e., $m_{i}^{t}<0$, follows from the fact that only a single fixed point exists with positive means [9, Lemma 2.3]; thus the second fixed point, if it exists, must have negative means (see Example 2). Third, it is an immediate consequence of the update rule in (6) that the consistent fixed point is energetically favored, i.e., $\mu^{s}_{ij}(x_{j})>\mu^{t}_{ij}(-x_{j})$. Then, computing the marginals according to (8) implies that the consistent fixed point has larger means, i.e., $|m_{i}^{s}|>|m_{i}^{t}|$. The same line of reasoning shows that the correlations are larger as well. Finally, the consistent fixed point has a lower average energy, i.e., $E_{{B}}^{s}({\tilde{P}_{B}})<E_{{B}}^{t}({\tilde{P}_{B}})$ (see (11)). Therefore (and by symmetry of the Bethe entropy), it follows that $\mathcal{F}_{B}^{*}=\mathcal{F}_{B}^{s}$. ∎

Example 2**.**

Consider a unidirectional homogeneous model on a complete graph with $N=4$, positive couplings $J_{ij}=J>0$, and local potentials specified by $\theta_{i}=\theta=0.02$. Again, we enforce symmetry by setting all singleton marginals to the same value and compute the pairwise marginals and the Bethe free energy accordingly.

Similar as in Example 1, we illustrate the effect of the coupling strength on the shape of $\mathcal{F}_{B}$. In some sense we see a similar picture emerging; i.e., $\mathcal{F}_{B}$ is convex for small values of $J$ (see Fig. 3a) and non-convex for larger values of $J$ (see Fig. 3c).111111 Contrary to Example 1, no closed-form solution exists for computing the fixed points because of the local potentials. There is one crucial difference however: because of $\theta_{i}>0$ the minimum does not turn into a maximum with two minima branching of. Instead, we see that the global minimum persists for all values of $J$, whereas the energy landscape becomes flat before two additional stationary points emerge somewhere else (see Fig.3b).

We will see in Section 5 that unidirectional (non-homogeneous) models behave qualitatively the same; consequently we will make repeated use of the behaviour illustrated in Example 2 for analysing the behaviour of SBP.

Finally, note that every fixed point that corresponds to a local minimum is a stable fixed point of BP for unidirectional models. This is a direct consequence of [50, Theorem 6] and the fact that unidirectional models are attractive per definition. Even more importantly, BP will always converge for unidirectional models: therefore note that BP always converges for unidirectional models with $\theta_{i}=0$ [33] and that the existence of local potentials with $\theta_{i}\neq 0$ can only enhance the convergence properties [21]. Recently, it has also been shown that for unidirectional models BP will converge reasonably fast, if initialized properly [24, 9].

2.5.1 Combination of Fixed Points

If multiple fixed points exist, the performance of BP often varies considerably between different fixed points. Here we briefly introduce the RSB (replica symmetry breaking) assumption that states how all those fixed points can be combined to compute the exact marginals. Therefore, let us consider all local minima of the Bethe free energy $\mathcal{F}_{B}^{m}$, the pseudomarginals ${\tilde{P}_{B}^{m}}$ at the corresponding fixed point, and the associated partition function approximation $\mathcal{{Z}_{B}}^{m}$.

Conjecture 1** (RSB Assumption).**

Let us consider all $M$ local minima of the Bethe free energy $\mathcal{F}_{B}^{m}$, the pseudomarginals at the corresponding fixed point ${\tilde{P}_{B}^{m}}$, and the associated partition function approximation $\mathcal{{Z}_{B}}^{m}$. Then, the exact singleton marginals $P(x_{i})$ are given by a convex combination of all fixed points at minima of $\mathcal{F}_{B}$ according to

[TABLE]

This representation belongs to a set of postulates in [31]. One underlying assumption is that the system does exhibit multiple fixed points (unique fixed points would falsely imply exact marginals otherwise). Despite its non-rigorous flavor, Conjecture 1 has been verified for a wide range of problems [30, Ch.19]. Example 1 and Fig. 2c in particular give an intuitive illustration of the RSB assumption: for this model the local maximum (i.e., the unstable fixed point) corresponds to the exact marginals; yet, both local minima do not approximate the marginals well. Combining both minima according to (15), however, gives the exact singleton marginals because of the model’s symmetry properties.

Obtaining all minima of the Bethe free energy, however, is a complex task that is only possible for models with certain structures or potentials [6, 64, 20, 36]. Besides the application to optimization problems [4] the RSB assumption is thus still of limited practical relevance for estimating the marginals; nonetheless, it provides a powerful concept for assessing and comparing the marginals given a selection of fixed points [22].

3 Self-guided belief propagation (SBP)

The main aim of SBP is to select an accurate fixed point, i.e., a fixed point that approximates the singleton marginals well, or – if none are stable – to approximate an accurate fixed point. In this section we present an intuitive justification of the proposed method and subsequently introduce SBP in detail. We further present practical considerations and pseudocode of SBP. A formal treatment of SBP is presented in Section 5.

The current understanding of BP is that strong pairwise potentials negatively influence BP in the sense of deteriorating the convergence behavior [32, 20] and the approximation quality of the marginals [20, 22, 57] and that incorporating the potentials slowly [3] may reduce the overall number of iterations. Inspired by the recent observation that strong local potentials increase accuracy and lead to better convergence properties [21], we aim to reduce the influence of the pairwise potentials that negatively influence BP. In doing so, we hope that an accurate fixed point emerges if we start from a simple model with independent random variables and slowly increase the potentials’ strength.

SBP starts from a simple model with independent random variables and slowly incorporates the edge potentials’ strength. More precisely, SBP is a homotopy method that solves the simple problem first and – by repetitive application of BP – keeps track of the fixed point as the interaction strength is increased by a scaling term.

Formally, SBP considers an increasing length-$M$ sequence $\{\zeta_{m}\}$ where $m=1,\ldots,M$ such that $\zeta_{m}<\zeta_{m+1}$ and $\zeta_{m}\in[0,1]$ with $\zeta_{1}=0$ and $\zeta_{M}=1$. This further indexes a sequence of probabilistic graphical models $\{\mathcal{U}_{m}\}$ that converges to the model of interest $\mathcal{U}_{M}=\mathcal{U}$. Every probabilistic graphical model has a set of potentials $\Psi_{m}=\{\Phi_{m}(x_{i},x_{j}),\Phi_{m}(x_{i})\}$ associated, where $\Phi_{m}(x_{i})=\Phi(x_{i})$ and the pairwise potentials at index $m$ are exponentially scaled by

[TABLE]

We further write $\bm{\mu}^{\circ}_{m}$ to clarify that we consider a BP fixed point for the model $\mathcal{U}_{m}$. If multiple fixed points exist, the initialization of BP plays a crucial role in determining to which fixed point BP converges and how it performs. SBP always provides a favorable initialization by the preceding fixed point and aims to perform the composite function121212 The underlying assumption is the existence of a continuous path that connects the initial messages $\bm{\mu}^{\text{init}}_{{1}}$ to the final fixed point $\bm{\mu}^{\circ}_{M}$, where all intermediate initializations $\bm{\mu}^{\text{init}}_{m}$ lie along this path. For now, we assume a fixed number of models $M$ for practical reasons but we will consider SBP’s behavior in the limit of $M\rightarrow\infty$ in our formal analysis in Section 5.

[TABLE]

This may lead to problems if the fixed point becomes unstable for some value $m<M$, in which case we cannot rely on BP to keep track of the fixed point anymore. Instead, SBP provides the last stable fixed point in that case, i.e., $\bm{\mu}^{\circ}_{m}$ as the final estimate.

Example 3**.**

We illustrate how SBP approximates the marginals for a frustrated model, for which BP fails to converge, in Fig. 4. Initially, SBP obtains the pseudomarginals for $\zeta=0$ by running BP on the simple model. As $\zeta$ increases, SBP estimates the pseudomarginals along the solution path by running BP to keep track. For $\zeta>0.7$, however, BP does not converge anymore; SBP consequently terminates and provides the marginals of the last stable fixed point as the final approximate solution. Note that the approximated marginals are already close to the exact ones in this example; experiments show that this is often the case (see Section 4).

In other words, SBP relaxes the problem of minimizing $\mathcal{F}_{B}$ by making all variables independent (and the Bethe approximation exact). Then the problem is deformed into the original one by increasing $\zeta$ from zero to one. Thereby, a stationary point $\mathcal{F}_{B}^{\circ}$ emerges as a well-behaved path (see Property 1 in Section 5) and SBP keeps track of it with BP constantly correcting the stationary point.

3.1 Practical considerations

In practice the run-time of SBP is influenced by the difference between two successive fixed points $\bm{\mu}^{\circ}_{m}$ and $\bm{\mu}^{\circ}_{m-1}$ – the difference is primarily determined by the number of steps $M$. Ideally $M$ should be as large as possible. This, however, increases the run-time (see Theorem 6); in practice, we would choose $M$ as small as possible but as large as necessary. Moreover, one can adaptively increase the step size if two successive fixed points are close, i.e., if $\bm{\mu}^{\circ}_{m}\simeq\bm{\mu}^{\circ}_{m-1}$ (cf. [43, pp.23], [1]). Our experiments show that it is sufficient to use rather coarse steps; we used $M\leq 10$ for all reported experiments and using more steps would not have improved the accuracy.

Additionally, instead of initializing $\mathcal{BP}_{m}$ with its preceding fixed point messages, i.e.,

$\bm{\mu}^{\text{init}}_{m}=\bm{\mu}^{\circ}_{m-1}$ one can (e.g., by spline extrapolation) estimate $\bm{\mu}^{\text{init}}_{m}=f(\bm{\mu}^{\circ}_{m-1},\bm{\mu}^{\circ}_{m-2},\ldots,\bm{\mu}^{\circ}_{m-k})$ so that $\bm{\mu}^{\text{init}}_{m}\cong\bm{\mu}^{\circ}_{m}$ to reduce the overall number of iterations. We empirically observed that the benefit diminishes for $k>3$.

3.2 Pseudocode

Pseudocode of SBP is presented in Algorithm 1. Note that we initialize the messages $\bm{\mu}^{\text{init}}$ randomly. The sequence of messages over all models $\{\mathcal{U}_{m}\}$ is collected in $\{\bm{\mu}^{\circ}_{m}\}=\{\bm{\mu}^{\circ}_{1},\ldots,\bm{\mu}^{\circ}_{m}\}$. ExtrapolateMsg applies cubic spline extrapolation to estimate the initial messages of the subsequent model.

We further present the pseudocode for the adaptive step size controller in Algorithm 2. The step size is increased depending on $k$, where $k$ is the number of preceeding fixed point messages that are $\varepsilon$-close to the current fixed point messages $\bm{\mu}^{\circ}_{m}$

4 Experiments

We apply SBP to attractive (Section 4.2) and general (Section 4.3) models on $n\times n$ grid graphs of different sizes, on complete graphs with $N=10$ random variables, and on random graphs (i.e., Erdos Renyi graphs with an average degree $\hat{d}=3$) with $N=10$ random variables. These graphs are considered in order to render the computation of the exact marginals feasible and to make the results comparable to previous work [57, 44, 29, 45].

4.1 Experimental settings

We evaluate SBP and compare it to BP, $\text{BP}_{\text{D}}$ (BP with damping), and Gibbs sampling (run for $10^{5}$ iterations) using the mean squared error (MSE) between the approximate marginals ${\tilde{P}_{X_{i}}}$ and the exact marginals $P_{X_{i}}$, where, for binary random variables, $\text{MSE}=\frac{2}{N}\sum_{i=1}^{N}|P_{X_{i}}(+1)-{\tilde{P}_{X_{i}}}(+1)|^{2}$. We further compare the run-time by counting the overall number of BP iterations and the number of iterations for Gibbs sampling.131313 Computing the acceptance-probability requires similar run-time as one BP message update.

For BP and SBP we set the maximum number of iterations to $N_{BP}=10^{3}$ and use random scheduling. We randomly initialize $\bm{\mu}^{\text{init}}$ and use an adaptive step size with an initial step size of $step_{init}=0.1$ and a threshold of $\varepsilon=10^{-3}$. The overall number of steps is thus $M\leq 10$.

For $\text{BP}_{\text{D}}$ we choose a large damping factor $\epsilon=0.9$ in order to prioritize convergence over run-time and therefore increase the maximum number of iterations to $N_{BP}=10^{4}$. Carefully selecting a damping factor that depends on a given model may reduce the number of iterations until convergence but cannot increase the accuracy; moreover, if chosen too small, $\text{BP}_{\text{D}}$ may fail to converge at all.

The initial messages are randomly initialized 100 times for each model, before applying BP with and without damping. We consider BP (and $\text{BP}_{\text{D}}$) as converged for a model if at least a single message initialization (out of 100) exists for which BP converges. The convergence ratio is the number of experiments (or probabilistic graphical models) for which BP converged at least once divided by the overall number of experiments (i.e., $100$).

The reported error (MSE) and the number of iterations are averaged over all convergent runs of BP and $\text{BP}_{\text{D}}$ (i.e., BP$\circ$ and $\text{BP}_{\text{D}}^{\circ}$) while all runs that did not converge are discarded. On the contrary, we average the error and the number of iterations over all models for SBP ($\text{SBP}_{all}$), Gibbs sampling ($\text{Gibbs}_{all}$), and for minimization of the Bethe approximation ($\mathcal{F}_{B,all}^{*}$).

4.2 Attractive models

We consider grid graphs with $N=10\times 10$ random variables, random graphs with $N=10$ random variables, and complete graphs with $N=10$ random variables. We generate $L=100$ models for every value of $\beta\in\{0,0.5,\ldots,5\}$ and sample the potentials according to $\theta_{i}\sim\mathcal{U}(-0.5,0.5)$ and $J_{ij}\sim\mathcal{U}(0,\beta)$; i.e., overall we consider 1100 different models for each individual graph-structure. Note that BP is randomly initialized 100 times for every considered model. We compute the MSE for every value of $\beta$ and visualize the mean and the standard deviation of the MSE141414 Note that the MSE is not Gaussian distributed but we report the standard deviation for simplicity. as well as the number of iterations in Fig. 5a.

BP (orange) converges rapidly for all graphs considered; hence, there is no additional benefit in using $\text{BP}_{\text{D}}$ (green) that only increases the number of iterations. SBP (blue) slightly increases the number of iterations as compared to BP but converges in fewer iterations than $\text{BP}_{\text{D}}$. Note that SBP is guaranteed to capture the global optimum when applied to unidirectional models (see Theorem 9-10). But even if we do allow for random local potentials, we empirically observe that SBP consistently outperforms BP with respect to accuracy; specifically for models with strong couplings These models exhibit multiple stable fixed points [20] such that, depending on the initialization, BP often converges to inaccurate fixed points while SBP selects an accurate one.

4.3 General models

General models with random potentials are often frustrated and traditionally pose problems for BP and other methods that aim to minimize the Bethe approximation.

First, in order to evaluate the performance of SBP we consider uniform local potentials $\theta_{i}=\theta\in\{0,0.1,0.4\}$ and draw the couplings with equal probability from $J_{ij}\in\{-1,1\}$; the results are summarized in Tab. I. Although BP and $\text{BP}_{\text{D}}$ fail to converge for most models we observe that SBP stops after only a few iterations and significantly outperforms BP in terms of accuracy. In fact, SBP achieves accuracy competitive with Gibbs sampling but requires three orders of magnitude fewer iterations.

Second, we apply SBP to general graphs and evaluate whether SBP provides a good approximation of the pseudomarginals $\tilde{P}_{B}^{*}$ that correspond to the global minimum of the Bethe free energy $\mathcal{F}_{B}^{*}$. Therefore we consider grid graphs (of size $5\times 5$), which still allows us to approximate the global minimum $\mathcal{F}_{B}^{*}$ – and the related pseudomarginals $\tilde{P}_{B}^{*}$ – reasonably well by [55]. The results are summarized in Tab. I and show that SBP approximates $\tilde{P}_{B}^{*}$ within the accuracy of our reference method ($\text{MSE}_{{B}}$). We further report the number of times where SBP obtains the terminal fixed point, i.e., for $\mathcal{U}_{M}$, in Tab. I i.e.,\bigl{(}\mathcal{F}_{B}^{\circ}(\zeta_{M}) equals SBP \bigr{)}. It becomes obvious that SBP approximates the terminal fixed point reasonably well, despite frequently stopping for $\zeta_{m}<1$. Moreover, the MSE reveals that SBP does not only approximate the pseudomarginals $\tilde{P}_{B}^{*}$ well but concurrently provides an accurate approximation of the exact marginals ${P_{B}}$.

Third, we investigate how the approximation quality depends on the scaling term $\zeta_{m}$. Therefore, we depict the evolution of $\text{MSE}_{{B}}$ (to the approximate solution) and of the MSE (to the exact solution) in Fig. 6. We observe that $\text{MSE}_{{B}}$ (blue) decreases monotonically with every iteration, which empirically verifies that SBP proceeds along a well-behaved solution path (see Property 1). Note that $\text{MSE}_{{B}}$ decreases rapidly in the first iterations and SBP spends a major part of the overall run-time for slight improvements. The MSE to the exact solution, on the other hand, decreases first until it increases again as SBP incorporates stronger couplings. Stronger couplings tend to degrade the quality of the Bethe approximation in loopy graphs and lead to marginals that are increasingly biased towards one state [51, 20]. This explains why the MSE to the exact solution increases as SBP converges towards the terminal fixed point. One could exploit this behavior and restrict the run-time by stopping SBP after consumption of a fixed iteration budget; this may even increase the accuracy with respect to the exact solution.

Finally, we investigate the influence of the coupling strength: therefore we consider $\theta_{i}\sim\mathcal{U}(-0.5,0.5)$ and $J_{ij}\sim\mathcal{U}(-\beta,\beta)$. For every $\beta\in\{0,0.5,\ldots,5\}$ we consider $L=100$ models and present the averaged results in Fig. 5b. We restrict the results to $\beta\leq 2$ on the grid graph because BP did only converge sporadically for models with stronger couplings. Even for models where BP converged, SBP requires only slightly more iterations than BP and fewer than $\text{BP}_{\text{D}}$. The benefits of SBP become increasingly evident as the coupling strength increases. Again SBP (blue) significantly outperforms BP$\circ$ (orange) and $\text{BP}_{\text{D}}^{\circ}$ (green) on all graphs with respect to accuracy.

5 Theoretical Properties of SBP

Here, we present some more formal arguments and discuss the properties of SBP to understand under which conditions the algorithm (presented in Section 3) can be expected to perform well. We begin by formally defining the notion of a solution path in Section 5.1 before we discuss the properties of SBP’s solution path in Section 5.2. Finally, we restrict ourselves to unidirectional models, for which we analyze the accuracy of SBP in Section 5.3.

5.1 Solution Path: Definition

First, we fix our notation: For the model $\mathcal{U}_{m}$ with its potentials $\Psi_{m}$ we refer to the pseudomarginals by ${\tilde{P}_{B}^{\circ}}(\zeta_{m})$, and, with slight abuse of notation, we refer to the corresponding stationary point of the Bethe free energy by $\mathcal{F}_{B}^{\circ}(\zeta_{m})=\mathcal{F}_{B}({\tilde{P}_{B}^{\circ}}(\zeta_{m}))$. It is beneficial to study the behavior of SBP as $M$ tends towards infinity. Therefore we consider the unit interval $\zeta\in[0,1]$ to be the compact support of the functions $\mathcal{F}_{B}(\zeta)$ and ${\tilde{P}_{B}}(\zeta)$. SBP is inspired by the idea to proceed along a so-called solution path as $\zeta$ increases from zero to one in order to obtain the marginal distributions for the model of interest. Therefore, we consider a continuous homotopy function $H(\bm{\mu},\zeta):\mathbb{R}^{|\bm{\mu}|+1}\rightarrow\mathbb{R}^{|\bm{\mu}|}$. Note that the scaling factor modifies the graphical model $\mathcal{U}=(\mathcal{G},\Psi(\zeta)$, for which the pairwise potentials are specified according to (16); we define the homotopy accordingly by

[TABLE]

Consequently, for a given value $\zeta$, $H(\bm{\mu},\zeta)$ is zero precisely for the fixed points $\bm{\mu}^{\circ}$ of $\mathcal{U}=(\mathcal{G},\Psi(\zeta)$. It follows that a solution path

[TABLE]

exists that (i) has a start point $c(\zeta=0)=\bm{\mu}:H(\bm{\mu},\zeta=0)=0$, (ii) an endpoint $c(\zeta=1)=\bm{\mu}:H(\bm{\mu},\zeta=1)=0$ , and (iii) is continuous over $\zeta\in[0,1]$, i.e., the solution path connects the start- with the endpoint.

SBP then proceeds along a solution path $c(\zeta)$ that implicitly defines the pseudomarginals ${\tilde{P}_{B}^{\circ}}(\zeta)$ by (8)-(9). In particular, we refer to the start- and endpoint by ${\tilde{P}_{B}^{\circ}}(\zeta=0)$ and ${\tilde{P}_{B}^{\circ}}(\zeta=1)$ respectively.

The following example illustrates the solution path for a unidirectional model.

Example 4**.**

Consider a unidirectional grid graph of size $3\times 3$ with $J=2$ and $\theta=0.1$. For this model, we can compute the fixed points of BP; i.e., the solutions to $H(\bm{\mu},\zeta)=0$ with the numerical homotopy continuation method as in [20].

Fig. 7 shows how the solutions evolve as $\zeta$ goes from zero to one. Similar to Example 2, this illustrates how increasing the interactions between variables changes the energy landscape. That is, for small values of $\zeta$ we have a unique fixed point until – as $\zeta$ increases – two additional fixed points emerge; one of which is a minimum and one a maximum of $\mathcal{F}_{B}$. Although multiple solutions to $H(\bm{\mu},\zeta)=0$ exist for larger values of $\zeta$, all except for one lack a start point and are therefore of no relevance for any method that proceeds along a solution path defined by the homotopy in (18).

5.2 Solution Path: Properties

The following proposition summarizes the main properties of the solution path that is specified and followed by SBP.

Properties 1** (Convergence Properties).**

SBP proceeds along a well-defined solution path $c(\zeta)$. More precisely, we have:

1. (1)

BP has a unique fixed point $\bm{\mu}^{\circ}_{1}$ for $\zeta=0$, so that SBP has a unique start point ${\tilde{P}_{B}^{s}}(\zeta=0)$ (see Theorem 3). 2. (2)

A smooth (i.e., continuous) solution path originates from ${\tilde{P}_{B}^{s}}(\zeta=0)$ (see Theorem 4). 3. (3)

SBP proceeds along this (unique) solution path $c(\zeta)$ and terminates after a fixed number of iterations (see Theorem 6).

In the sequel, we validate all points of Property 1.

5.2.1 Uniqueness of the Start Point

First, we show that the initial model with $\zeta=0$ has a unique fixed point. This is an immediate (and known) consequence of the fact that all variables are independent; we still include the following Theorem for the sake of completeness.

Theorem 3** (Property 1.1).**

*For the initial model $\mathcal{U}_{1}$ with parameters $\Psi(\zeta=0)$, BP has a unique fixed point $\bm{\mu}^{\circ}_{1}$. That is, only a single start point exists for the solution path of SBP. Moreover, this start point is exact, i.e., ${\tilde{P}_{B}^{\circ}}(\zeta=0)={P_{B}}(\zeta=0)$, and BP is guaranteed to converge.151515 Note that this concurs with the sandwich-bound [54, Th.4] that reduces to $P_{X_{i}}(+1)={\theta_{i}}/(e^{\theta_{i}}+e^{-\theta_{i}})={\tilde{P}_{X_{i}}}(+1)$ for $J_{ij}=0$. *

Proof:

First, we show that all messages converge to the identical value whenever all $J_{ij}=0$. Therefore, note that we can omit the pairwise potential in the update equation of BP and rewrite (6) to

[TABLE]

It follows that (20) is independent of the state $x_{j}$, i.e., the computed messages are the same for all states of $X_{j}$. Thus, after normalization, $\mu_{ij}(x_{j})=\frac{1}{|\mathcal{X}|}=0.5$ for all edges $(i,j)\in\mathbf{E}$.

Second, the pseudomarginals are then obtained according to

[TABLE]

Finally, we must show that the exact marginals are the same, i.e, that $P(x_{i})=\tilde{P}(x_{i})$. Therefore, note that the joint in (5) factorizes only over the local potentials, which allows us to reorder the summation and compute the exact marginals according to

[TABLE]

∎

Theorem 3 thus reduces the problem of initializing SBP to computing the fixed point messages $\bm{\mu}^{\circ}_{1}$ of the initial model with $\zeta=0$; this can be done in linear time. Moreover, a unique start point implies that only a single solution path emerges from it.

5.2.2 Smooth Solution Path

Second, we show that – for any model – the Bethe free energy is a smooth function of the scaling parameter $\zeta$ in the sense that no discontinuities exist. This implies the existence of a smooth solution path; let us make this notion of a smooth solution path more precise:

Theorem 4** (Property 1.2).**

Consider the pseudomarginals ${\tilde{P}_{B}^{\circ}}(\zeta)$ along the (one) solution path that originates from the start point ${\tilde{P}_{B}^{\circ}}(\zeta=0)$. Then, both the pseudomarginals ${\tilde{P}_{B}^{\circ}}(\zeta)$ and the associated stationary point $\mathcal{F}_{B}^{\circ}(\zeta)$ are continuous functions of the scaling parameter $\zeta$ along the whole solution path, i.e., for $\zeta\in[0,1]$.

Proof:

First, we show that the Bethe free energy $\mathcal{F}_{B}({\tilde{P}_{B}},\zeta)$ is continuously differentiable. Therefore, the derivative with respect to $\zeta$ must exist and be continuous as well. Note that $\mathcal{F}_{B}({\tilde{P}_{B}},\zeta)$ is a high-dimensional function defined over ${\tilde{P}_{B}}\in\mathcal{L}$, i.e., over all locally consistent marginals. Now, let us consider (13) with the pairwise potentials defined by (16). Then, when taking the derivative with respect to $\zeta$, all terms depending only on singleton marginals vanish and the derivative is given according to

[TABLE]

Since (24) is a finite sum over finite terms, this proves that $\mathcal{F}_{B}(\zeta)$ is continuously differentiable.161616 Moreover – according to its definition in (13) – $\mathcal{F}_{B}({\tilde{P}_{B}},\zeta)$ is a polynomial, which guarantees that it is not just continuously differentiable but in fact an analytical function [39]. Note that this is in accordance with the fact that true phase transitions (singularities in the derivative of the free energy) can only occur in the thermodynamic limit, where (24) becomes a sum over infinitely many terms that equates to infinity. More precisely, (24) proves that $\mathcal{F}_{B}({\tilde{P}_{B}},\zeta)$ is continuously differentiable with respect to $\zeta$ for every point on the energy surface.

We specifically consider the local minimum $\mathcal{F}_{B}^{\circ}(\zeta)$ along the solution path, i.e., the minimum that emerges from the unique start point $\mathcal{F}_{B}^{\circ}(\zeta=0)$.

Given that the Bethe free energy varies in a continuous fashion along this solution path (for $\zeta\in[0,1]$), the local minimum $\mathcal{F}_{B}^{\circ}(\zeta)$ and its position vary in a continuous fashion as well [8].171717 This is actually not restricted to problems with a unique minimum; instead sets of minima are considered in [8] to prove that minima of constrained minimization problems vary continuously under continuous changes of the objective function (see also [38, Section 5]). For existence of this fixed point, note that minima of $\mathcal{F}_{B}$ persist when increasing $\zeta$ (similar to increasing the inverse temperature; see [63]). Consequently, because of the one-to-one correspondence between stationary points of the Bethe free energy $\mathcal{F}_{B}^{\circ}(\zeta)$ and fixed points of BP ${\tilde{P}_{B}^{\circ}}(\zeta)$, the pseudomarginals are also continuous. Note, however, despite each minimum varying continuously with increasing $\zeta$, the global minimum may not (as a former local minimum becomes the global minimum).

∎

In Theorem 4, we have shown that continuity of the Bethe free energy $\mathcal{F}_{B}(\zeta)$ and the corresponding pseudomarginals ${\tilde{P}_{B}^{\circ}}(\zeta)$ implies a smooth solution path. We implicitly assumed that this further implies a well-behaved solution path. In general this is the case but there are two scenarios that might be problematic as we will discuss now:

First, it is absolutely critical that the stationary points of $\mathcal{F}_{B}$ are zero-dimensional, i.e., that every minimum is a distinct point, and that there are finitely many stationary points. Fortunately, this is the case and the number of stationary points of $\mathcal{F}_{B}$ is always finite [50]; the same is true for the fixed points of BP (for as long as the messages are properly normalized) [27].

Second, a potential issue is how SBP copes with bifurcations of the solution path. A bifurcation is a point along the solution path at which a minimum $\mathcal{F}_{B}^{\circ}$ turns into a maximum with two new minima branching off; i.e., one minimum splitting up into two new minima (see Example 1 for the illustration of a typical bifurcation). In that case, SBP will follow the evolution of a single minimum and discard the other one. Which fixed point is selected mainly depends on the implementation details of BP (this is because we rely on BP to find the next fixed point after our prediction). Note that different minima are usually wide apart (see Example 1 or [22]). If, however, two minima would be close together, SBP might switch between branches because of numerical issues (e.g. because of quantization errors or a step-size too large). Since this is not a fundamental property of the solution path, we will assume sufficient numerical precision and rule out this phenomenon. Thus, the existence of bifurcations is no particular problem and SBP retains the property of a unique solution path.

Although bifurcations are not problematic from a practical perspective, it is still relevant for our theoretical analysis of SBP to understand if and where bifurcations of the solution path exist. Unfortunately, the current state of knowledge is still relatively limited in this respect. Note, however, that at least for unidirectional models it is actually possible to make precise statements regarding the existence of bifurcations.

Lemma 5**.**

Let $\mathcal{G}$ be a unidirectional model. Then, if the local potentials of all nodes are specified by $\theta_{i}=0$, a bifurcation occurs. Else, if some local potentials are nonzero, no bifurcation exists along the solution path.

Proof:

Since we consider unidirectional models there are two minima of $\mathcal{F}_{B}$ at most. We denote these minima by $\mathcal{F}_{B}^{s}$ and $\mathcal{F}_{B}^{t}$. We will explain now, why these two minima branch off the previously unique minimum as we increase the scaling parameter $\zeta$ if and only if $\theta_{i}=0$. Therefore, we denote the critical scaling parameter, for which the two minima begin to exist, by $\zeta_{c}$. Then, by definition, a bifurcation exists if and only if both minima coincide for $\zeta_{c}$, i.e, if

[TABLE]

First, remember that both fixed points are symmetric if all $\theta_{i}=0$, i.e., $m_{i}^{s}=-m_{i}^{t}$ (see Example 1 and Fig. 2). Because of this symmetry, both fixed points will bifurcate from a fixed point with zero mean $m_{i}=0$. Note, that this specific case is well understood and serves as the archetype for the existence of bifurcations on the models we consider in this work [33].

Second, it remains to show that this is not the case if $\theta_{i}\neq 0$. Therefore, we show that both minima do not coincide for $\zeta_{c}$, i.e.,

[TABLE]

This is an immediate consequence of the fact that both fixed points are not symmetric anymore (see Section 2.5). In fact, $m_{i}^{s}>0$ is strictly positive, whereas $m_{i}^{t}<0$ is strictly negative; as a consequence, both fixed points can never coincide, which rules out the existence of bifurcations along the solution path. ∎

The implications of Lemma 5 are illustrated in Example 1 and 2. Note how even small values for $\theta_{i}\neq 0$ break the symmetry between both fixed points and consequently imply the absence of bifurcations.

To summarize, Theorem 4 and Lemma 5 guarantee the existence of a smooth solution path that connects the start point ${\tilde{P}_{B}^{\circ}}(\zeta=0)$ to the endpoint ${\tilde{P}_{B}^{\circ}}(\zeta=1)$. This is a property of fundamental importance as, at least in principle, it allows path-tracking algorithms (as SBP) to proceed along this well-behaved solution path. The only remaining question is whether SBP is actually capable of proceeding along the solution path and if it obtains the endpoint; i.e., if SBP converges.

5.2.3 Convergence Properties of SBP

So far, we have substantiated the claim that a smooth solution path emerges from the trivial start system leading to the target system. The existence of this smooth solution path alone, however, is not sufficient to guarantee that SBP obtains the marginals of the target system. As already discussed in Section 3, SBP relies on BP to proceed along the solution path; ideally, keeping the overall amount of BP iterations low. If – at some point (for $\zeta<1$) – the fixed point ${\tilde{P}_{B}^{\circ}}(\zeta)$ becomes unstable, this approach will inevitably break down. Nonetheless, SBP will proceed along the solution path for as long as it has a stable fixed point and approximates the marginals precisely at the onset of instability. That is, even if the terminal fixed point ${\tilde{P}_{B}^{\circ}}(\zeta=1)$ is unstable, SBP will always converge and approximate the marginals in a limited number of iterations.181818 This is in stark contrast to plain BP that fails to converge and thus cannot be used to approximate the marginals in such a case.

Theorem 6** (Property 1.3).**

There exists some $\zeta\!\leq\!1$ so that SBP converges to $\tilde{P}_{B}^{\circ}(\zeta)\in\mathcal{L}$ in $\mathcal{O}(MN_{BP}\!)$.

Proof:

SBP increases $\zeta$ as long as BP converges in less than $N_{BP}$ iterations, and stops otherwise. Consequently, BP corrects the accuracy of the fixed point for each value $\zeta_{m}$ within a bounded number of iterations. The run-time of SBP is further determined by the choice of $M$, i.e., the step-size (see Section 3.1). That is, BP is run until convergence for at most $M$ times and thus, SBP converges in at most $M\cdot N_{BP}$ iterations. ∎

SBP is consequently capable of tracking the fixed point that emerges as $\zeta$ increases and requires $MN_{BP}$ iterations at most. SBP may, however, only converge to a surrogate model for $\zeta_{m}<1$ and is not guaranteed to obtain the pseudomarginals of the endpoint. One can characterize this error by computing a bound on $|\mathcal{F}_{B}^{\circ}(\zeta_{m})-\mathcal{F}_{B}^{\circ}(\zeta=1)|$ given the difference between $\Psi(\zeta_{m})$ and $\Psi(\zeta=1)$ (cf. [16, Theorem 16]).

SBP works particularly well for attractive models, for which every minimum of $\mathcal{F}_{B}$ is also a stable fixed point of BP [50]. More precisely, if for some scaling term $\zeta_{m}$ the fixed point of the solution path becomes unstable, i.e., if it turns into a local maximum, two stable fixed points will branch of [50, Theorem 6]. Thus for attractive models (and also for balanced models) a smooth solution path with a stable fixed point always exits so that SBP is consequently capable of tracking the solution path until its end.

Corollary 6.1**.**

Let $\mathcal{G}$ be an attractive or unidirectional model. Then, every point along the solution path going from $\zeta=0$ to $\zeta=1$ is stable. Consequently, some number of steps $M$ exists for which SBP converges and obtains the marginals ${\tilde{P}_{B}^{\circ}}(\zeta=1)$.

Proof.

According to Theorem 4, a smooth solution path $c(\zeta)$ exists that goes from the start point ${\tilde{P}_{B}^{\circ}}(\zeta=0)$ to the endpoint ${\tilde{P}_{B}^{\circ}}(\zeta=1)$, where every point along $c(\zeta)$ corresponds to a local minimum of the Bethe free energy $\mathcal{F}_{B}^{\circ}$. Every local minimum is also stable for attractive models [50]. Therefore, the correction by BP will converge in every step if initialized sufficiently close to the fixed point, i.e., if the step size is chosen small enough. ∎

5.3 Accuracy of SBP

The convergence properties of SBP (Property 1) are of fundamental importance for understanding SBP. Yet, so far we have only been concerned if SBP converges but neglected any discussion regarding the quality of the obtained solution. Thus, we will now theoretically analyze how well the fixed point obtained by SBP approximates the marginals. Unfortunately, however, our understanding of what characterizes a fixed point that approximates the marginals well is relatively poor. Indeed, it is not even obvious how the accuracy of the approximated marginals relates to the accuracy of the approximated free energy [22].

Therefore, we will mainly focus our discussion on the simplest model class for which the solution space is well understood (i.e., unidirectional models). Note that for unidirectional models SBP is actually guaranteed to converge (see Section 2.5). Regarding the accuracy of SBP – when applied to unidirectional models – we summarize the properties as follows:

Properties 2** (Accuracy of SBP for Unidirectional Models).**

Let $\mathcal{G}$ be a unidirectional model, i.e., all edges are attractive ($J_{ij}>0$) and all local potentials are specified by $\theta_{i}$ with the same sign. Without loss of generality we consider the case where all $\theta_{i}\geq 0$.191919 Note that identical results can be obtained in a straightforward manner for $\theta_{i}<0$ because of symmetry properties. Then, the solution path ends at a fixed point that approximates the marginals and the free energy well. More specifically, SBP either obtains the exact marginals (if all $\theta_{i}=0$) or finds the global minimum of the Bethe free energy $\mathcal{F}_{B}^{*}$ that – for unidirectional models – corresponds to the fixed point with the most accurate marginals. Let us now break down the argument into the following three points:

1. (1)

For $\theta_{i}=0$, SBP obtains the exact marginals, i.e., ${\tilde{P}_{B}^{\circ}}={P_{B}}$ (see Theorem 7). 2. (2)

For $\theta_{i}\neq 0$, SBP finds the global minimum of $\mathcal{F}_{B}(\zeta)$ and gives the best approximation to the free energy $\mathcal{F}$ (see Theorem 9). 3. (3)

For $\theta_{i}\neq 0$, the global minimum of $\mathcal{F}_{B}(\zeta)$ also gives the best approximation of the marginals, i.e., no other minimum of $\mathcal{F}_{B}(\zeta)$ has more accurate marginals (see Theorem 10).

First, we show that SBP obtains the exact marginals in the special case where all local potentials are specified by $\theta_{i}=0$

Theorem 7** (Property 2.1).**

Consider a unidirectional model with $\theta_{i}=0$. Then, the pseudomarginals obtained by SBP ${\tilde{P}_{B}^{s}}(\zeta=1)$ are identical to the exact ones ${P_{B}}(\zeta=1)$.

Proof:

For attractive models with $\theta_{i}=0$ it is straightforward to compute the exact marginals. The exact marginals of all variables are uniformly distributed and have zero mean $m_{i}(\zeta)=0$ for all values of $\zeta$; note that these marginals also correspond to a stationary point $\mathcal{F}_{B}^{\circ}(\zeta)$ [33]. For sufficiently small values of $J_{ij}$ this stationary point is a (unique) minimum but for large values of $J_{ij}$ it turns into a maximum (see [30, pp.385] and Example 1).

Even though the fixed point with $m_{i}(\zeta)=0$ is unstable, SBP will still converge to this fixed point. This is rather surprising, as – in general – BP can neither converge to an unstable fixed point nor would it remain there because of quantization errors. There are two reasons why SBP obtains the potentially unstable fixed point in the special case of models with $\theta_{i}=0$ nonetheless. First, note that the (unique) fixed point of the initial mode is identical to the exact solution; hence all variables have zero mean at the start point, i.e., $m_{i}^{s}(\zeta=0)=0$. Second, all messages $\mu_{ij}(x_{j})\in\bm{\mu}^{\circ}(\zeta=0)$ take the same value of $\mu^{\circ}_{ij}(x_{j})=1/2$, since $\theta_{i}=0$. The extrapolated messages will stay the same for every iteration of SBP because – as discussed above – they already constitute a valid (albeit possibly unstable) fixed point (see Example 1). Note, however, that these fixed point messages can be represented without quantization error in binary arithmetic. Thus SBP will remain on the fixed point corresponding to the exact marginals. ∎

The above statement arguably holds only for very restricted models and is thus of limited practical relevance. Indeed, not only does standard BP obtain the same fixed point if initialized with identical messages but the exact marginals are also known and can be computed analytically if all $\theta_{i}=0$, obviating the need for approximate inference methods.

Our next results will discuss the accuracy of SBP for the more general class of unidirectional models both with respect to approximating the marginals and the Free energy (or partition function respectively).

Therefore, we first describe how the scaling term $\zeta$ affects the marginals along the solution path $c(\zeta)$ by generalizing Griffiths’ inequality  [13] to the fixed points of BP. In particular, this analysis reveals that the means $m_{i}^{\circ}$ are monotonically increasing with $\zeta$.

Lemma 8** (Monotonicity of the Means).**

Consider two attractive models $\mathcal{U}_{m}$ and $\mathcal{U}_{n}$ on the same graph $\mathcal{G}$, associated scaling parameters $\zeta_{m}<\zeta_{n}$, and local potentials specified by $\theta_{i}>0$. Note that we are effectively considering two nearly equivalent models with $\mathcal{U}_{n}$ having strictly larger couplings than $\mathcal{U}_{m}$. Now, let us consider a BP fixed point of $\mathcal{U}_{m}$ with positive means $m_{i}^{\circ}(\zeta_{m})>0$. Then, the corresponding fixed point of $\mathcal{U}_{n}$ has larger means $m_{i}^{\circ}(\zeta_{n})>m_{i}^{\circ}(\zeta_{m})$ and correlations $\chi_{ij}^{\circ}(\zeta_{n})\geq\chi_{ij}^{\circ}(\zeta_{m})$.

Proof:

The main ingredient of this proof is to show how the ratio between the messages for their respective states $\bar{\mu}_{ij}=\frac{\mu_{ij}(X_{j}=+1)}{\mu_{ij}(X_{j}=-1)}$ increases monotonically with $\zeta$. Then, the fact that the means and the correlations increase monotonically follows directly from the definition of the pseudomarginals in terms of the messages (see (8) and (9)).

First, we show that the message ratio $\bar{\mu}_{ij}$ increases monotonically: therefore, let us denote the messages on both models by $\mu_{(m)ij}(x_{j})$ and $\mu_{(n)ij}(x_{j})$ respectively. Furthermore, without loss of generality, we assume that all couplings of $\mathcal{U}_{n}$ are $\varepsilon$-larger than those of $\mathcal{U}_{m}$. Then, for every $\varepsilon$ there exists $\delta>0$ so that whenever $0<\bar{\mu}_{(n)ij}-\bar{\mu}_{(m)ij}<\delta$, we have $0<J_{(n)ij}-J_{(m)ij}<\varepsilon$. Note that by assumption $m_{i}\in(0,1]$ so that

[TABLE]

First, we show that for all $(i,j)\in\mathbf{E}$

[TABLE]

Therefore, consider the update rules of (6) for both states

[TABLE]

[TABLE]

In (29) the larger product is multiplied by $e^{\varepsilon}$ and the smaller product is divided by $e^{\varepsilon}$. For (30) it is exactly the other way round so that the ratio between the messages increases which proofs (28). We shall denote the imposed difference $\delta^{\prime}\in\mathbb{R}_{+}^{*}$ on the messages by

[TABLE]

Second, we show that $m_{i}^{\circ}(\zeta_{n})>m_{i}^{\circ}(\zeta_{m})$, which is an immediate consequence of plugging (31) and (32) into (8) (see Example 1 that illustrates how all fixed points, i.e., minima of $\mathcal{F}_{B}$, are pulled towards $m_{i}^{\circ}$ extreme values of zero and one as $J$ increases).

Finally, it remains to show that $0\stackrel{{\scriptstyle(i)}}{{\leq}}\chi_{ij}^{\circ}(\zeta_{m})\stackrel{{\scriptstyle(ii)}}{{\leq}}\chi_{ij}^{\circ}(\zeta_{n})$. Without loss of generality we assume that all variables have equal degree $d+1$, constant coupling strength $J_{ij}=J$, and constant local fields $\theta_{i}=\theta$. First we show that (i) holds, i.e., $\chi=\chi_{ij,0}^{\circ}$ is positive. Let us express the marginals by (9) and denote the messages by $\mu=\mu_{(m)ij}(X_{j}=1)$. It follows that $\mu_{(m)ij}(X_{j}=-1)=(1-\mu)$ and that

[TABLE]

Let us further represent the messages by $\mu=1/2+\beta$ with $\beta\in[0,1/2]$. It follows that

[TABLE]

where $(a)$ follows from neglecting all exponential terms and thus upper bounding the positive term and lower bounding the negative term (with equality if and only if $J=0$ and $\theta=0$) and $(b)$ is a direct consequence of the square in (34). Now let us show that (ii) holds, i.e., $\chi$ increases monotonically, by taking the derivative of (33), so that

[TABLE]

where $(a)$ follows from neglecting the, strictly positive, first term in (36), and $(b)$ is a direct consequence from (27). ∎

Theorem 9** (Property 2).**

Consider a unidirectional model, i.e, an attractive model with $\theta_{i}>0$. Then the minimum of the Bethe free energy along the solution path $\mathcal{F}_{B}^{s}(\zeta)$ is always the global minimum, i.e., $\mathcal{F}_{B}^{s}(\zeta)=\min\mathcal{F}_{B}(\zeta)$. Moreover, the Bethe free energy of this stationary point decreases monotonically with $\zeta$.

Proof:

A unique start point ${\tilde{P}_{B}^{s}}(\zeta=0)$ exists by Theorem 3. For $\theta_{i}>0$ this fixed point has only positive means $m_{i}^{s}(\zeta=0)>0$ and correlations $\chi_{ij}^{s}(\zeta=0)>0$.

Consequently, Lemma 8 applies, which implies that the corresponding means $m_{i}^{s}(\zeta)$ and correlations $\chi_{ij}^{s}(\zeta)$ are monotonically increasing along the solution path. This further implies that the Bethe free energy $\mathcal{F}_{B}^{\circ}(\zeta)$ decreases along the solution path. Let $\mathcal{F}_{B}^{s}(\zeta=1)$ correspond to the endpoint of the solution path $c(\zeta)$ that emerges from the origin; then, it immediately follows that the error with respect to the endpoint $\mathcal{F}_{B}^{s}(\zeta=1)$ decreases along the solution path: i.e., consider two arbitrary values $m,n\in[0,1]$ such that $n>m$, then $|\mathcal{F}_{B}^{s}(\zeta_{m})-\mathcal{F}_{B}^{s}(\zeta=1)|\geq|\mathcal{F}_{B}^{s}(\zeta_{n})-\mathcal{F}_{B}^{s}(\zeta=1)|$.

It remains to show that SBP obtains the fixed point $\mathcal{F}_{B}^{s}(\zeta)$ that is the global minimum of the Bethe free energy for a given scaling parameter $\zeta$. Therefore, consider the fact that unidirectional models have two fixed points at most, only one of which has positive means (see Lemma 2). As discussed above, SBP proceeds along a solution path belonging to the consistent fixed point with $m_{i}^{s}(\zeta)>0$. The second fixed point $t$, if it exists, has negative means $m_{i}^{t}(\zeta)<0$. Given that SBP obtains the consistent fixed point, it follows from Lemma 2 that SBP obtains the global minimum, i.e., $\mathcal{F}_{B}^{t}\geq\mathcal{F}_{B}^{s}=\mathcal{F}_{B}^{*}$. ∎

Theorem 9 makes it explicit that – for unidirectional models – SBP favors the fixed point that defines the global minimum of $\mathcal{F}_{B}$. This is in accordance with the observation that the average energy is linear in the pseudomarginals (see (11)), which further implies that varying the local potentials essentially tilts the energy landscape (cf. [37]). Note that this behavior is also illustrated in Example 1 and 2 accordingly: i.e, as opposed to Fig. 2, local potentials with $\theta_{i}>0$ tilt the energy landscape in Fig. 3 so that the consistent fixed point constitutes the global minimum of $\mathcal{F}_{B}$.

Although finding the global minimum of $\mathcal{F}_{B}$ is desirable, one is ultimately interested in approximating the exact free energy as good as possible. Fortunately, both properties coincide for unidirectional models.

Corollary 9.1** (Accuracy of the Bethe Approximation).**

Consider a unidirectional model. Then there is no other fixed point of BP ${\tilde{P}_{B}^{k}}$ that approximates the free energy $\mathcal{F}$ better than the one obtained by SBP, i.e.,

[TABLE]

Proof:

This is a direct consequence of SBP obtaining the global minimum of the Bethe free energy (see Theorem 9) and the fact that the Bethe free energy upper bounds the free energy for attractive models (cf. [40]). ∎

This implies that the fixed point obtained by SBP minimizes the approximation error of the free energy $|\mathcal{F}_{B}^{k}(\zeta)\!-\!\mathcal{F}^{*}\!|$. Moreover – because of the direct relationship between the Bethe free energy and the Bethe approximation of the partition function – this implies that the fixed point obtained by SBP also approximates the partition function best.

Note, however, that statements about the approximation quality of the free energy do not necessarily extend to the accuracy of the pseudomarginals [22]. Nonetheless, at least for unidirectional models, this is the case and SBP obtains the fixed point with the most accurate marginals.

Theorem 10** (Accuracy of Marginals).**

Consider a unidirectional model and assume that the RSB assumption holds (see Conjecture 1). Then there is no other fixed point of BP ${\tilde{P}_{B}^{k}}$ that approximates the marginals better than the one obtained by SBP, i.e,

[TABLE]

Proof:

Note that we are only considering unidirectional models here. For these models, although not formally verified, the RSBP assumption is accepted to hold. Thus, if multiple fixed points do exist this allows us to express the exact marginals as a convex combination of all fixed points according to Conjecture 1.

Moreover, remember that the Bethe free energy of unidirectional models has two minima at most (see Section 2.5). This – and the fact that the convex combination of both fixed points yields the exact marginals by assumption – implies that the fixed point obtained by SBP (i.e., the one that minimizes the Bethe free energy) approximates the marginals best. ∎

To conclude, the fixed point along the solution path is always stable for unidirectional models. Thus, SBP follows the solution path to its endpoint. In doing so, SBP not only decreases the Bethe free energy and obtains the global minimum but it also provides the fixed point with the most accurate marginals.

6 Conclusion

In this paper, we introduced an iterative algorithm to perform approximate inference: self-guided belief propagation (SBP) is a simple and robust method that gradually accounts for the pairwise potentials and guides itself towards a unique, stable, and accurate solution.

We provide a comprehensive theoretical analysis in order to validate the underlying assumptions: in particular, we showed that: (i) a smooth solution path exists and originates from a unique start point that is obtained by neglecting the pairwise potentials; (ii) this solution path is well-behaved and can be tracked by SBP; and (iii) for unidirectional models, the solution of SBP approximates the exact solution well and corresponds to the global optimum of the Bethe approximation.

The theoretical analysis does not fully answer for which models SBP is expected to perform particularly well. Although SBP obtains the most accurate fixed point for unidirectional models, the marginals will not be approximated well for models with strong couplings, for the simple reason that accurate fixed points do not exist at all. We expect that SBP is particularly advantageous for all models that have a large amount of fixed points with varying accuracy. Thus SBP will perform well for dense models with random couplings (particularly if they are centered around zero) for which the marginals are largely affected by the local potentials. Not only do we expect that SBP approximates the marginals well in that case but it is actually essential for estimating marginals as BP will often fail to converge.

In order to verify SBP empirically, we applied SBP to a range of models, not only attractive but also general ones. In our experiments, SBP consistently improves the accuracy in comparison with BP (with and without damping). Besides, SBP approximates the exact marginals well on graphical models for which BP does not converge at all.

Acknowledgment

The authors greatly appreciate the valuable discussions with Bernhard C. Geiger, KNOW-Center, and acknowledge his help in making the formal definition of SBP much more concise. Moreover, we thank Alexander Ihler, University of California Irvine, for the encouraging discussion during UAI’17 that was pivotal for actually pursuing the – at this time – rather unpolished idea. We further appreciate the help of Florian Kulmer in implementing and empirically evaluating our algorithm and are particularly indebted to Harald Leisenberger for carefully reading the manuscript.

This work was supported by the Graz University of Technology LEAD project “Dependable Internet of Things in Adverse Environments”. Adrian Weller acknowledges support from the David MacKay Newton research fellowship at Darwin College, The Alan Turing Institute under EPSRC grant EP/N510129/1 & U/B/000074, and the Leverhulme Trust via the Leverhulme Centre for the Future of Intelligence (CFI).