# Circular Belief Propagation for Approximate Probabilistic Inference

Vincent Bouttier, Renaud Jardri, Sophie Deneve

► To cite this version:

Vincent Bouttier, Renaud Jardri, Sophie Deneve. Circular Belief Propagation for Approximate Probabilistic Inference. 2026. hal-05537654

HAL Id: hal-05537654

https://hal.science/hal-05537654v1

Preprint submitted on 5 Mar 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Circular Belief Propagation for Approximate
Probabilistic Inference

Vincent Bouttier¹,²*, Renaud Jardri¹,², Sophie Denève¹†

¹Group for Neural Theory, Ecole Normale Supérieure, Paris, France

²Université de Lille, Lille, France

renaud.jardri@univ-lille.fr

Abstract

Belief Propagation (BP) is a simple probabilistic inference algorithm, consisting of passing messages between nodes of a graph representing a probability distribution. Its analogy with a neural network suggests that it could have far-ranging applications for neuroscience and artificial intelligence. Unfortunately, it is only exact when applied to cycle-free graphs, which restricts the potential of the algorithm. In this paper, we propose Circular Belief Propagation (CBP), an extension of BP which limits the detrimental effects of message reverberation caused by cycles by learning to detect and cancel spurious correlations and belief amplifications. We show in numerical experiments involving binary probabilistic graphs that CBP far outperforms BP and reaches good performance compared to that of previously proposed algorithms.

# 1 Introduction

Probabilistic inference arises in many applications, from medical diagnosis to Wi-Fi protocols. In all methods performing inference, the Belief Propagation algorithm [Pearl, 1988] is among the most successful ones. This procedure consists of propagating local information by passing messages between nodes of a graph representing the probability distribution. Belief Propagation is known to achieve exact inferences when the graph is acyclic. The key reason for this lies in the fact that a message $m_{i \to j}$ is computed based on all messages $m_{k \to i}$ received by the sender node except the opposite message $m_{j \to i}$, therefore preventing this message from being reverberated and thus counted twice. One drawback of BP is that it often performs poorly in cyclic graphs, for which correcting for these loops of length 2 is no longer sufficient [Murphy et al., 1999, Weiss, 2000]. Messages pass in cycles from node to node and return to the original node, causing the same piece of information to be counted several times. This is known as the “double-counting” problem. Several variants or generalization of BP have been proposed to tackle this problem ever since the interest in BP started to grow [Minka, 2001, Sudderth et al., 2003, Ihler and McAllester, 2009].

We present in this work the Circular Belief Propagation algorithm. This algorithm proposes to counter the effect of cycles (which introduce spurious correlations between messages) by actively decorrelating the messages. We show that Circular BP significantly outperforms BP in cyclic probabilistic graphs, including in cases where BP does not converges, and performs approximate inference with a very impressive quality even for fully dense probabilistic graphs.

This paper is organized as follows. In Sections 2-5, we provide background information, define the Circular BP algorithm, relate it to existing algorithms, and provide sufficient conditions for its convergence. In sections 6 and 7, we show the performance of CBP on both synthetic and real-world problems involving binary variables, using supervised or unsupervised procedures. The Appendices

<footnote>*Now at Iktos SA.</footnote>

<footnote>†Now at Institut de Neurosciences de la Timone.</footnote>
provide, among others, formulations of the algorithm in the general case (non-pairwise graphs and continuous variables).

# 2 Background

## 2.1 Inference in probabilistic graphical models

The object of study is a probability distribution $p(\mathbf{x})$, where $\mathbf{x} = (x_1, x_2, \dots, x_n)$. The distribution can be decomposed into a product of conditionally independent factors that each describes interactions between distinct subsets of variables [Koller and Friedman, 2009, Wainwright and Jordan, 2008]:

$$p(\mathbf{x}) = \frac{1}{Z} \prod_{(i,j)} \psi_{ij}(x_i, x_j) \prod_i \psi_i(x_i) \quad (1)$$

where we consider only unitary and pairwise interactions (Appendices show an extension of this special case to any factorization of $p(\mathbf{x})$, that is, any Markov random field involving higher-order potentials, e.g., $\psi_{ijk}(x_i, x_j, x_k)$). $\psi_i$ represents prior knowledge about variable $x_i$, while $\psi_{ij}$ describes the interaction between $x_i$ and $x_j$. $Z$ is a normalization constant. As Figure 1A shows, the probability distribution can be represented graphically as a graphical model called *factor graph*, composed of variable nodes $x_i$ and factor nodes $\psi_{ij}$ and $\psi_i$.

We focus in this work on the particular inference problem of finding the marginals of $p(\mathbf{x})$, $p_i(x_i) \equiv \sum_{\mathbf{x}\setminus x_i} p(\mathbf{x})$, given partial and noisy information (incorporated into $\{\psi_i(x_i)\}$) about $\mathbf{x}$. Unfortunately, direct calculation of the marginals $p_i(x_i)$ takes exponential time in the number of variables $n$. That is why algorithms, like the *Belief Propagation algorithm* and sampling methods, have been developed to solve this inference problem more efficiently. However, their drawback is that they produce only approximate marginal probabilities $b_i(x_i) \approx p_i(x_i)$.

## 2.2 Belief Propagation

Belief Propagation [Pearl, 1988], or sum-product algorithm, is a message-passing algorithm which performs approximate inference on a probabilistic graph. It approximates the marginals of the distribution by making variable nodes $x_i$ share all the probabilistic information available with the rest of the network by sending messages to other variable nodes, eventually propagating information to the whole network. The algorithm consists of running iteratively the following update message equation on the graph (see Figure 1B):

$$m_{i \to j}^{\text{new}}(x_j) \propto \sum_{x_i} \psi_{ij}(x_i, x_j) \psi_i(x_i) \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}^{\text{old}}(x_i) \quad (2)$$

where $\mathcal{N}(i)$ is the set of neighbors of node $x_i$ in the graph. Messages are for instance initialized uniformly over the nodes ($m_{i \to j}(x_j) = 1/\mathcal{N}(j)$). Once messages have converged, approximate marginal probabilities (or *beliefs*) are computed as:

$$b_i(x_i) \propto \psi_i(x_i) \prod_{k \in \mathcal{N}(i)} m_{k \to i}(x_i) \quad (3)$$

A crucial feature of the BP algorithm is the message exclusion principle ($k \in \mathcal{N}(i) \setminus j$ in Equation (2)): to compute $m_{i \to j}$, all messages $m_{k \to i}$ coming to node $x_i$ are taken into account and combined, except the message in the opposite direction $m_{j \to i}$. In other words, node $x_i$ sends to $x_j$ everything it knows ($b_i$) except what $x_j$ communicated to $x_i$. This feature intuitively explains why BP produces exact marginals in the case of acyclic graphs: information propagates into the network without “coming back”. However, in the presence of cycles, the same evidence travels multiple times through loops of the graph, and is mistaken for new evidence [Pearl, 1988]; see legend of Figure 1. This intuitively explains why BP (also called *Loopy BP* in this case) is incorrect in the general case. More precisely, the algorithm does not necessarily converge on cyclic graphs, and when it converges it produces incorrect marginals. A number of empirical studies [Murphy et al., 1999, Weiss, 2000, Mooij and Kappen, 2004, Litvak et al., 2009] have led to a better understanding of the convergence and performance of BP depending on the probabilistic graph.
<figure><img src="image_2.png" /><figcaption>Figure 1: Belief Propagation and Circular Belief Propagation algorithms applied to a probabilistic graph. (A) The probability distribution $p(x)$ is represented by a factor graph with pairwise potentials $\psi_{ij}$ and unitary potentials $\psi_i$. (B) BP aims at estimating marginals $p_i(x_i)$ by exchanging messages in the graph. The message $m_{1\to2}$ depends on three components: the messages received by node $x_1$ from its neighbors except $x_2$, including the external message (see full black lines), and the interaction $\psi_{12}$. Estimated marginal $b_i(x_i)$ of a node $x_i$ is formed based on all messages received by the node. BP is not exact when applied to cyclic graphs, for two reasons. First, messages get counted multiple times: $m_{1\to2}$ naturally travels back to $x_1$ because of the cycle $x_1 - x_2 - x_3 - x_4$. Second, opposite messages are correlated: $m_{1\to2}$ depends on $m_{5\to1}$ which depends on $m_{4\to5}$ which depends on $m_{1\to4}$ which depends on $m_{2\to1}$. (C) Contrary to BP, Circular BP (partially) takes $m_{2\to1}$ into account to compute $m_{1\to2}$. Parameter $\kappa$ fights the belief amplifications caused by messages being reverberated. Parameters $\alpha$ decorrelates opposite messages.</figcaption></figure>

BP has initially been used because of its empirical performance, but lacked a deep theoretical comprehension as to why the algorithm worked so well in practice. However, in the early 2000s, a theoretical foundation of BP emerged: stable fixed points of BP are minima of the Bethe free energy from statistical physics [Yedidia et al., 2001, Heskes, 2002]. The Bethe free energy approximates the Gibbs free energy $G(b) := D_{KL}(b||p) - \log(Z) = \sum_x b(x)E(x) + \sum_x b(x)\log(b(x))$ where the energy $E(x) := -\log(p(x)Z)$, by considering that the variational distribution $b(x)$ can be represented by an acyclic probabilistic graph (Bethe approximation), as further explained in the Appendices. Minimizing the Gibbs free energy also minimizes the distance (KL divergence) between $b$ and $p$, and the idea (or hope) is that marginals of $b$ and $p$ are therefore similar: $b_i(x_i) \approx p_i(x_i)$. This tight link between BP and the Bethe free energy opened the way to investigate the goodness of BP through the estimation of the BP error [Wainwright et al., 2003, Taga and Mase, 2006, Ihler, 2007, Mooij and Kappen, 2009, Shi et al., 2010], and the properties of the Bethe free energy and its fixed points [Heskes et al., 2002, Heskes, 2004, Watanabe and Fukumizu, 2009, Watanabe, 2011, Weller et al., 2014, Knoll and Pernkopf, 2017].

# 3 Circular Belief Propagation

## 3.1 Definition

The Circular Belief Propagation (CBP) algorithm is an extension of the Belief Propagation (BP) algorithm and performs message-passing on probabilistic graphs to carry out approximate inference. It consists of initializing messages $m_{i\to j}$ to some distribution, and then running the following updates:

$$m_{i \to j}^{\text{new}}(x_j) \propto \sum_{x_i} \psi_{ij}(x_i, x_j)^{\beta_{ij}} (\psi_i(x_i)^{\gamma_i} \times \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}^{\text{old}}(x_i) \times m_{j \to i}^{\text{old}}(x_i)^{1-\alpha_{ij}/\kappa_i})^{\kappa_i} \quad (4)$$

After convergence (or a given number of iterations), approximate marginals are computed using:

$$b_i(x_i) \propto (\psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i)} m_{k \to i}(x_i))^{\kappa_i} \quad (5)$$

BP corresponds to the particular case $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$. The main difference between CBP and BP is the occurrence of the message in the opposite direction $m_{j\to i}$ to compute $m_{i\to j}$; see Figure 1B. Note that $\alpha_{ij}$ and $\beta_{ij}$ are assigned to the *undirected edge* $(i, j)$, while $\kappa_i$ and $\gamma_i$ are assigned to the variable node $x_i$.
## 3.2 Circular BP for binary distributions

For probability distributions $p(x)$ over binary variables ($x_i \in \{-1, +1\}$), CBP takes a very simple form in the log-domain. For simplicity, we consider here a particular class of such distributions called *Ising models* (also known as *Boltzmann machines*) from statistical physics [Ising, 1925, Baxter, 1982] in which pairwise factors take a specific form: $\psi_{ij}(x_i, x_j) \propto \exp(J_{ij}x_i x_j)$. CBP writes as follows:

$$ M_{i \to j}^{\text{new}} = f(B_i - \alpha_{ij} M_{j \to i}, \beta_{ij} J_{ij}) \quad (6) $$

$$ B_i = \kappa_i \left( \sum_{j \in \mathcal{N}(i)} M_{j \to i} + \gamma_i M_{\text{ext} \to i} \right) \quad (7) $$

(see Appendix E.6 for the demonstration). $M_{i \to j} \equiv \frac{1}{2} \log \left( \frac{m_{i \to j}(x_j = +1)}{m_{i \to j}(x_j = -1)} \right)$ represents the information about variable $x_j$ brought by variable $x_i$, and $B_i \equiv \frac{1}{2} \log \left( \frac{b_i(x_i = +1)}{b_i(x_i = -1)} \right)$ is by definition half of the log-likelihood ratio. The approximate marginal probabilities are given by $b_i(x_i = \pm 1) = \sigma(\pm 2B_i)$, i.e., $b_i(x_i) \propto \exp(B_i x_i)$. Lastly, $M_{\text{ext} \to i} \equiv \frac{1}{2} \log \left( \frac{\psi_i(x_i = +1)}{\psi_i(x_i = -1)} \right)$ represents alternatively prior information over a variable $x_i$, or here a noisy sensory observation providing information about the variable (ext stands for “external”).

Function $f$ takes a simple form (see also Mooij and Kappen [2007a]):

$$ f(x, J) = \operatorname{arctanh} \left( \tanh(J) \tanh(x) \right) \quad (8) $$

$f$ is a sigmoidal function of $x$ controlled by parameter $J$ representing the bounded “trust” $J$ between nodes; see Figure S1B. Note that $f$ resembles the usual non-linearity $f(x, w) = w \tanh(x)$.

## 3.3 Interpretation of Circular BP

Parameters $\alpha$, $\kappa$, $\beta$ and $\gamma$ can be interpreted easily. $\beta$ and $\gamma$ act as rescaling factors (respectively for the pairwise and unitary factors). $\kappa$ rescales the beliefs in order to counter the effect of message amplifications in loops (double-counting); see Figure 1C. But the specificity of BP is the parameter $\alpha$ which controls how the opposite message $m_{j \to i}$ is taken into account for the computation of $m_{i \to j}$: positively ($\alpha_{ij} < \kappa_i$), negatively ($\alpha_{ij} > \kappa_i$) or not at all ($\alpha_{ij} = \kappa_i$) as in BP. $\alpha$ therefore acts as a loop correction factor. Equation (6) means that node $i$, encoding for variable $x_i$, sends to node $j$ everything it knows ($B_i$) except a rescaled version of what $j$ communicated to $i$ ($M_{j \to i}$).

# 4 Related work

## 4.1 Circular BP from Jardri and Denève [2013]

A particular case of Circular BP was defined in Jardri and Denève [2013], which corresponds to $(\kappa, \beta, \gamma) = (1, 1, 1)$. The important conceptual difference is that the authors propose a way of impairing BP, not of improving it; see more in the discussion.

## 4.2 Link to reweighted Bethe free energy

Circular BP tightly relates to the family of reweighted BP algorithms for approximate inference, including Fractional BP [Wiegerinck and Heskes, 2002], Reweighted BP [Loh and Wibisono, 2014], Tree-reweighted BP [Wainwright. et al., 2002, 2003, Wainwright et al., 2005], Convex BP [Hazan and Shashua, 2008, 2010], Power Expectation Propagation [Minka and Lafferty, 2002, Minka, 2004], and $\alpha$-BP [Liu et al., 2019, 2020] (see the Appendices for more detail). All these reweighted BP algorithms are variational inference techniques which aim at minimizing a particular reweighted Bethe free energy approximating the Gibbs free energy. First, the variational entropy component $-\sum_x b(x) \log(b(x))$ of the Gibbs free energy is approximated as if the variational distribution $b(x)$ was:

$$ b(\mathbf{x}) \approx \prod_{i,j} \left( \frac{b_{ij}(x_i, x_j)}{b_i(x_i) b_j(x_j)} \right)^{1/\alpha_{ij}} \prod_i (b_i(x_i))^{1/\kappa_i} \quad (9) $$
(which corresponds to BP's Bethe approximation for $\alpha = 1$ and $\kappa = 1$). Second, the variational average energy component $\sum_x b(x)E(x)$ of the Gibbs free energy is approximated as if $p(x)$ was:

$$p(x) \approx \prod_{i,j} (\psi_{ij}(x_i, x_j))^{\beta_{ij}} \prod_i (\psi_i(x_i))^{\gamma_i} \quad (10)$$

Altogether, these hypotheses result in a parametric approximation of the Gibbs free energy, which generalizes the Bethe free energy from BP. The algorithm aiming at minimizing this resulting approximate Gibbs free energy writes for Ising models:

$$M_{i \to j}^{\text{new}} = \frac{1}{\alpha_{ij}} f(B_i - \alpha_{ij} M_{j \to i}, \alpha_{ij} \beta_{ij} J_{ij}) \quad (11)$$

$$B_i = \kappa_i \left( \sum_{j \in \mathcal{N}(i)} M_{j \to i} + \gamma_i M_{\text{ext} \to i} \right) \quad (12)$$

where function $f$ is the same as above. Message update equations are similar between Reweighted BP algorithms and Circular BP (see Equations (6) and (7)) and therefore so are the results of the algorithms (see Figure S11). The algorithms differ in the specific form of the message update for parameter $\alpha$. In Reweighted BP algorithms, $\alpha$ corresponds to a weight in the reweighted Bethe free energy, the quantity that these algorithms try to minimize. Instead, in Circular BP, $\alpha$ relates to the idea of cancelling the dependencies between opposite messages arising because of loops (see Figure 1C).

To find parameters for which these Reweighted BP algorithms outperform BP, several techniques are used. Tree-Reweighted BP [Wainwright et al., 2005] uses a gradient ascent method to optimize $\alpha$, but involves solving a max-weight spanning tree problem at each iteration. Fractional BP [Wiegerinck and Heskes, 2002] fits $\alpha$ by using an unsupervised learning rule based on linear response theory.

# 5 Convergence properties of Circular BP

Convergence of Circular BP is crucial to carry out approximate inference. In general, if message-passing does not converge, then beliefs have very little to do with the correct marginals [Murphy et al., 1999]. It is therefore important to control the convergence properties, and if possible, to ensure the convergence of our proposed algorithm. We state here sufficient conditions for the convergence of Circular BP in an Ising model, and provide ways to find such parameters. Notably, Theorem 5.3 shows that choosing $\alpha$ and $\kappa$ uniformly, equal and small enough guarantees the convergence of CBP, for any set of parameters $(\gamma, \beta)$. More precisely, in this case, the algorithm has only one fixed point and converges to it with at least a linear rate.

We start by defining matrix $A$ whose coefficients are:

$$A_{i \to j, k \to l} = |\kappa_i| \tanh|\beta_{ij} J_{ij}| \delta_{il} \mathbb{1}_{\mathcal{N}(i)}(k) \left|1 - \frac{\alpha_{ij}}{\kappa_i}\right|_{j=k} \quad (13)$$

where $\mathbb{1}_{\mathcal{N}(i)}(k) = 1$ if $k \in \mathcal{N}(i)$, otherwise $0$, $\delta_{il} = 1$ if $i = l$, otherwise $0$, and $|1 - \frac{\alpha_{ij}}{\kappa_i}|_{j=k} = 1 - \frac{\alpha_{ij}}{\kappa_i}$ if $j = k$, otherwise 1.

**Theorem 5.1.** If for any induced operator norm $\|\cdot\|$ (sometimes called natural matrix norm), $\|A\| < 1$, then CBP has a unique fixed point and CBP converges to it with at least a linear rate.

**Theorem 5.2.** If $\forall (i, j)$, $\alpha_{ij}/\kappa_i \le 1$ and the spectral radius of $A$, $\rho(A) < 1$, then Circular BP converges to the unique fixed point.

The proofs of Theorems 5.1 and 5.2 are provided in the Appendices and follow closely the ones provided in Mooij and Kappen [2007a] for BP, that is, the special case $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$.

A consequence of Theorem 5.2 is the following fundamental result, which distinguishes Circular BP from related approaches like Power EP, Fractional BP and $\alpha$-BP:

**Theorem 5.3.** For a given weighted graph (with weights $J_{ij}$), it is possible to find parameters $\alpha$ and $\kappa$ such that Circular BP converges for any external input $M_{\text{ext}}$ and any choice of parameters $(\gamma, \beta)$.
*Proof.* Let us take $\alpha_{ij} = \kappa_i \equiv v \in \mathbb{R}_+$. In this case, $A_{i \to j, k \to l} = v \tanh|\beta_{ij} J_{ij}| \delta_{il} \mathbb{1}_{\mathcal{N}(i) \setminus j}(k)$. When $v \to 0$, all coefficients of A go to zero. The spectral norm (induced by the $l_2$-norm) is a continuous application, and $\rho(0) = 0$, thus $\rho(A) \to 0$ when $v \to 0$. We conclude by using Theorem 5.2 as $\alpha_{ij}/\kappa_i = 1 \le 1$: there exists $v^*$ such that for all $p < p^*$, CBP converges. $\square$

# 6 Approach for optimizing the parameters

The goal is, given the interactions $J = \{J_{ij}\}$ between variables of the probability distribution, to perform approximate inference for any external input $M_{\text{ext}} = \{M_{\text{ext} \to i}\}$.<fnref n="3" /> We consider two ways for finding the optimal parameters of a given weighted graph. The supervised learning method minimizes the MSE loss between the approximate and true marginals, but can only be applied to small or sparsely connected graph where exact inference can be performed on the training examples. Additionally, we propose an unsupervised learning method where parameters are estimated using local learning rules, and that can generalize to arbitrarily large and complex graphs.

## 6.1 Supervised learning procedure

Parameters of the model $(\alpha, \kappa, \beta, \gamma)$ are learnt so that the approximate marginals or beliefs $\{b_i(x_i)\}$ are as close as possible to the true marginals $\{p_i(x_i)\}$. We fit these $2n_{\text{nodes}} + 2n_{\text{edges}}$ parameters by minimizing the MSE loss or squared L2 norm $L(b, p) = \frac{1}{n_{\text{nodes}}} \sum_{i=1}^{n_{\text{nodes}}} [b_i(x_i = +1) - p_i(x_i = +1)]^2$ between beliefs $\{b_i(x_i)\}$ obtained with CBP and the true marginals $\{p_i(x_i)\}$. True marginals $p_i(x_i)$ are computed using the Junction Tree algorithm, an exact inference algorithm implemented on the pgmpy library [Ankan and Panda, 2015]. Approximate marginals $b_i(x_i)$ are obtained by running CBP for $T = 100$ time steps. The model is trained on PyTorch [Paszke et al., 2019] using backpropagation through time with a gradient-descent based method, resilient backpropagation [Riedmiller and Braun, 1993]. PyTorch does automatic differentiation to update parameters. We use a learning rate of 0.001 and stop the learning once the validation loss saturates. Parameters are initialized such that CBP converges with this choice of parameters (see section 5). One initialization tested, among others, is to initialize $\beta$ and $\gamma$ at the value 1, while $\alpha$ and $\kappa$ are initialized at the same small enough value $v$: we decrease $v$ from value 1 until $\rho(A) < 1$, which ensures that CBP converges (see Theorem 5.3). The model is also trained using function least squares of SciPy [Virtanen et al., 2020]. We select the parameters of the method (PyTorch or SciPy) performing better on the validation set.

## 6.2 Unsupervised learning procedure

The inconvenient of the supervised learning procedure is the need to generate training examples, i.e., to compute the true marginals. This has an exponential complexity in the size of the graph and therefore does not scale to big dense graphs. We propose a way of learning parameters $\alpha$ and $\kappa$ with:

$$ \Delta\alpha_{ij} = \eta_1 [M_{j \to i}(B_i - \alpha_{ij} M_{j \to i}) + M_{i \to j}(B_j - \alpha_{ij} M_{i \to j})] \quad (14) $$

$$ \Delta\kappa_i = -\eta_2 M_{\text{ext} \to i}(B_i - M_{\text{ext} \to i}) \quad (15) $$

where we keep $(\beta, \gamma) = (1, 1)$ fixed. Messages and beliefs are taken after $T = 100$ iterations, and $\eta_1$ and $\eta_2$ are learning rates. The second part of Equation (14) ensures that $\alpha$ is symmetric.

The proposed learning rules try to minimize redundancies in the information sent to other nodes. We place ourselves in a scenario where there should not be correlations between messages sent in opposite directions: training examples are pure noise, that is, are uninformative. The learning rule on $\alpha$ makes the information that node $i$ sends to $j$ $(B_i - \alpha_{ij} M_{j \to i})$ orthogonal to the information that node $j$ received from $i$ $(M_{j \to i})$. The learning rule on $\kappa$ aims at countering the amplification of inputs messages, so that overall the information received by node $i$ (including the effects of cycles) is $M_{\text{ext} \to i}$. In the resulting message-passing algorithm, redundancies caused by cycles are minimized. See Appendix C for more details.

<footnote n="3">³Note that the training does not aim at generalizing to different graph structures or different graph weights, as the fitted parameters $(\alpha, \kappa, \beta, \gamma)$ depend on the interactions $J$.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: **Results of Circular BP on Erdős-Renyi graphs.** Estimated marginals on the test set for BP and Circular BP, for both unsupervised and supervised learning procedures. One point represents the belief of a node, on one of the test examples, and one of the 30 randomly generated graphs.</figcaption></figure>

# 7 Numerical experiments

## 7.1 Synthetic problems

**Experimental setting** Simulations involve Erdős-Renyi graphs with 9 nodes (see Figure 2) generated with connection probabilities $p$ ranging from 0.2 to 1. For each graph topology or connection probability, 30 Ising models (weighted graphs) are generated. An absence of edge means $J_{ij} = 0$, while existing edges had their weights sampled randomly according to $J_{ij} \sim N(0, 1)$ (spin-glass). As a reminder, $J_{ij}$ is associated to the *unoriented edge* $(i, j)$: $J_{i \to j} = J_{j \to i} \equiv J_{ij}$. For each weighted graph, external evidences $M_{\text{ext}}$ are generated according to $M_{\text{ext} \to i} \sim N(0, 1)$. 200 training examples, 100 validation examples and 100 test examples are generated, where an example is a vector $M_{\text{ext}}$.

The Appendices generalize this restrictive setup to structured rather than random graphs (see Figure S5) as grids and bipartite graphs, graphs with a different distribution of weights (uniform positive couplings; see Figure S7), and graphs with a higher number of nodes (100 nodes; see Figure S8).

**Results of supervised learning** Marginals obtained after fitting the CBP algorithm can be seen in Figure 2. Circular BP performs well qualitatively: it generalizes very well to test data for all connection probabilities. On the contrary, for BP, there is a sharp transition from carrying out rather good approximate inference for $p = 0.2$ to performing really poorly for $p \ge 0.3$ (at least on some of the 30 weighted graphs) because the system does not converge (frustration phenomenon).

We compare quantitatively Circular BP to various approximate inference algorithms in Figure 3. Circular BP outperforms Fractional BP, itself outperforming BP. Logically, all algorithms show a decreased performance for increased graph density, but at a higher rate for BP. Additionally, we compared Circular BP to two more complex approximate inference algorithms: Loop-corrected BP [Mooij and Kappen, 2007b], which approximates the cavity distributions for each variable in a two-step way, and Double-Loop BP [Heskes et al., 2002], which uses constrained minimization of the Bethe free energy itself in a double-loop procedure guaranteeing the convergence of the algorithm. We used the libDAI implementation of these algorithms [Mooij, 2010] (option “LCBP_FULLCAVin_SEQRND” for LC BP and “HAK_BETHE” for DL BP). Circular BP performs worse than Loop-corrected BP although reaches comparable performance for dense graphs, and strongly outperforms Double-Loop BP.

**Results of unsupervised learning** The unsupervised learning rule achieves good performance in all tested cases; see Figure 2. We used 5000 training examples, and learning rates starting from $\eta_1 = 0.03$ and $\eta_2 = 0.0003$ (which were both decreased by half after one third of the optimization,
<figure><img src="image_4.png" /><figcaption>Figure 3: **Comparison between various algorithms.** Circular BP strongly outperforms the algorithms from the same family (Fractional BP, BP, Tree-Reweighted BP and Mean-Field) and the more complex and Double-Loop BP. CBP has comparable performance for dense graphs with another more complex algorithm, Loop-Corrected BP. The score measure is given in Equation (19).</figcaption></figure>

and after two thirds of the optimization). Note that we add some damping to the algorithm (see Appendices), $\epsilon = 0.7$, contrary to the supervised learning case.

## 7.2 Application to computer vision

We tackle a denoising problem on the MNIST database of handwritten digits [Deng, 2012]. Initial images {$\xi^k$} are vectors of size $n = 784$ (28 by 28 pixels) composed of gray levels initially from 0 to 255, rescaled between -1 and +1. Noisy images {$\tilde{\xi}^k$} have half of their pixels set at value zero.

To reconstruct images, we first use the Hopfield network [Hopfield, 1982] with weights defined by the pseudoinverse rule [Personnaz et al., 1986] based on $N = 60000$ training images: $\theta = \frac{1}{n}\xi\xi^T$ where $\xi^I$ is the pseudo-inverse of matrix $\xi$ (whose columns are the training examples {$\xi^k$}). The Hopfield network consists of initializing $x_0 = \tilde{\xi}^k$ and running $x_{t+1} = \text{sgn}(\theta x_t)$ until convergence; see Figure 4B for 5 examples of reconstruction of noisy test data.

Interestingly, the energy function of a Hopfield network has the same form as the one of an Ising model. Therefore, probabilistic inference methods can be used for this task. We apply BP to an all-to-all connected graph of 784 nodes with interactions $J = 5\theta$ (rescaling of the Hopfield weights) and external inputs $M_{\text{ext}}^k = 3\tilde{\xi}^k$ (rescaling of the noisy images); Figure 4B shows the marginals rescaled between -1 and +1. Further, we apply CBP, trained using the unsupervised learning procedure described above (with damping $\epsilon = 0.8$ on CBP) with 1000 training examples consisting of pure noise: $\tilde{\xi}^k \sim N(0, \text{Id})$. CBP performs better than BP and the Hopfield network; see Figure 4B. This shows the applicability of the unsupervised learning method on a large dense graph (784 nodes).

# 8 Discussions and conclusion

In this paper, we present an approximate inference method on undirected graphical models, the *Circular Belief Propagation* algorithm, defined for general probability distributions on any type of variables (discrete or continuous) and possibly with higher-order interactions.

The algorithm is easily interpretable: CBP consists of disturbing the exchange of information in cycles of length two in order to partly compensate for the detrimental effect of cycles of length 3 and more, which create spurious correlations between messages in BP. One can also talk of “anti-Circular BP”. We motivate the Circular BP algorithm theoretically, by relating it to Fractional BP and other
<figure><img src="image_5.png" /><figcaption>Figure 4: Potential application of Circular BP: Denoising hand-written digits for computer vision. Noisy data is presented to the Hopfield network / BP / CBP, which reconstruct a denoised version of the signal.</figcaption></figure>

Reweighted BP algorithms based on a reweighting of terms in BP's Bethe free energy, although Circular BP itself does not come from such an reweighting, nor minimizes a free energy.

We propose a way of learning the parameters of Circular BP using a supervised learning procedure (see section 7). However, one limitation of this procedure is the need to generate exact marginals, which is infeasible in large graphs. The unsupervised learning method to learn the corrective multiplicative factors brings a solution to this issue. It is less accurate than the supervised learning rule, but it is not limited by the network size and complexity. The unsupervised learning rules are intuitive: they relate to the principles of efficient information transmission, by trying to remove redundant probabilistic information created by cycles. They maintain the gain of activity and remove the prediction of the incoming information, similar to predictive coding [Rao and Ballard, 1999]. Indeed, Circular BP can be written $M_{i→j} = f(B_i - \hat{B}_i^j, \beta_{ij}J_{ij})$ where $\hat{B}_i^j$ is the prediction of the state of $x_i$ by variable $x_j$.

We show the applications of the model to a variety of problems involving binary variables, both synthetic and from the real world. Overall, CBP improves the quality of probabilistic inference in cyclic graphs compared to BP. It also converges on problems on which BP fails. CBP achieves near-optimal performance, even if BP does not converge (see Figure 2) thanks to the existence of parameters ensuring the convergence of the algorithm. In all cases, Circular BP hugely outperforms existing message-passing algorithms operating on the probabilistic graph itself, like BP and Fractional BP. We show in simulations, by considering particular probability distributions (Ising models) that the algorithm outperforms not only BP but also other message-passing algorithms such as Fractional BP [Wiegerinck and Heskes, 2002] and Power EP [Minka, 2004]. Results of inference using Circular BP are impressively good when using supervised learning.

**Applications of the method** The paper presents results in the restricted binary case. However, BP is used in this case in many real-world applications. This includes the MIMO (multiple input, multiple output) detection problem for wireless communication where variants of BP are in use in Wi-Fi and 5G protocols [Som et al., 2010, Yang et al., 2015, Liu and Li, 2018], error-correcting codes [McEliece et al., 1998], and computer vision [George et al., 2017], as well as many applications of Restricted Boltzmann Machines like sampling and learning [Lazaro-Gredilla et al., 2021]. The improvement of Circular BP over BP and its guarantee of existence of convergence conditions may enable approximate probabilistic inference in a wide range of new applications. Note that some applications could be harmful to society: better inference systems could be used either for the best or for the worst.

**Inference in the brain** It has been proposed that the BP algorithm is implemented in the brain to carry out inferences based on partial sensory (and/or cognitive) information [Shon and Rao, 2005, Denève, 2008, Steimer et al., 2009]. Circular BP, as an extension of BP, could also be implemented in the brain. Jardri and Denève [2013] use a particular case of CBP to account for aberrant beliefs in general, and psychosis in particular [Jardri and Denève, 2013, Denève and Jardri, 2016, Jardri et al., 2017, Leptourgos et al., 2017, 2021, Bouttier et al., 2021], based on the fact that $\alpha < 1$ brings
overconfidence to the estimated marginals compared to BP. Importantly, the algorithm from Jardri and Denève [2013] is not motivated normatively but only at the intuitive level [Denève, 2005, 2008] as a way of allowing for information reverberations with the underlying idea of excitation-inhibition imbalance (by introducing positive feedback between excitatory populations).

**Machine learning research** This work paves the way for future research on approximate inference. Message-passing algorithms have recently been receiving a lot of attention from the machine learning community with the success of Graph Neural Networks (GNN), Message-Passing Neural Networks, and variants inspired from BP [Gilmer et al., 2017, Kuck et al., 2020, Garcia Satorras and Welling, 2021]. We believe that the idea behind Circular BP (removing only part of the opposite message) can improve further these message-passing algorithms. In fact, Equation (6) provides a general recipe to use correction in message-passing algorithms, including mean-field methods, which typically write $x_i = \sum_j f(x_j, J_{ji})$ and could be translated into $x_i = \sum_j m_{ji}$ together with $m_{ji} = f(x_j - \alpha_{ji}m_{ij}, J_{ji})$. This "corrective multiplicative factor trick" can be used on top of algorithms which already improve BP while conserving the same form, like Kuck et al. [2020] or the message-GNN of Yoon et al. [2018], potentially leading to further improvement of these methods.

## References

Ankur Ankan and Abinash Panda. pgmpy: Probabilistic graphical models using python. In *Proceedings of the 14th Python in Science Conference (SCIPY 2015)*. Citeseer, 2015.

R. J. Baxter. *Exactly solved models in statistical mechanics*. 1982. ISBN 978-0-486-46271-4.

Vincent Bouttier, Suhrit Duttagupta, Sophie Denève, and Renaud Jardri. Circular inference predicts nonuniform overactivation and dysconnectivity in brain-wide connectomes. *Schizophrenia Research*, 2021. ISSN 0920-9964. doi: https://doi.org/10.1016/j.schres.2020.12.045.

A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the em algorithm. *Journal of the Royal Statistical Society, Series B*, 39(1):1–38, 1977.

Sophie Denève. Bayesian inference in spiking neurons. In *Advances in neural information processing systems*, 2005.

Sophie Denève. Bayesian Spiking Neurons I: Inference. *Neural Computation*, 20(1):91–117, 2008. ISSN 0899-7667. doi: 10.1162/neco.2008.20.1.91.

Sophie Denève and Renaud Jardri. Circular inference: mistaken belief, misplaced trust. *Current Opinion in Behavioral Sciences*, 11:40–48, oct 2016. ISSN 2352-1546. doi: 10.1016/J.COBEHA. 2016.04.001.

Li Deng. The mnist database of handwritten digit images for machine learning research. *IEEE Signal Processing Magazine*, 29(6):141–142, 2012.

Víctor Garcia Satorras and Max Welling. Neural enhanced belief propagation on factor graphs. In Arindam Banerjee and Kenji Fukumizu, editors, *Proceedings of The 24th International Conference on Artificial Intelligence and Statistics*, volume 130 of *Proceedings of Machine Learning Research*, pages 685–693. PMLR, 2021.

Dileep George, Wolfgang Lehrach, Ken Kansky, Miguel Lázaro-Gredilla, Christopher Laan, Bhaskara Marthi, Xinghua Lou, Zhaoshi Meng, Yi Liu, Huayan Wang, Alex Lavin, and D. Scott Phoenix. A generative vision model that trains with high data efficiency and breaks text-based captchas. *Science*, 358(6368):eaag2612, 2017. doi: 10.1126/science.aag2612.

Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural Message Passing for Quantum Chemistry. *34th International Conference on Machine Learning*, ICML 2017, 3:2053–2070, apr 2017.

Tamir Hazan and Amnon Shashua. Convergent Message-Passing Algorithms for Inference over General Graphs with Convex Free Energies. In *Uncertainty in Artificial Intelligence*, 2008.
Tamir Hazan and Amnon Shashua. Norm-product belief propagation: Primal-dual message-passing for approximate inference. *IEEE Transactions on Information Theory*, 56(12):6294–6316, dec 2010. ISSN 00189448. doi: 10.1109/TIT.2010.2079014.

Tom Heskes. Stable Fixed Points of Loopy Belief Propagation Are Minima of the Bethe Free Energy. In *Advances in neural information processing systems*, 2002.

Tom Heskes. On the uniqueness of loopy belief propagation fixed points. *Neural Computation*, 16 (11):2379–2413, 2004. ISSN 08997667. doi: 10.1162/0899766041941943.

Tom Heskes, Kees Albers, and Bert Kappen. Approximate inference and constrained optimization. In *Proceedings of the Nineteenth Conference on Uncertainty in Artificial Intelligence*, UAI'03, page 313–320, San Francisco, CA, USA, 2002. Morgan Kaufmann Publishers Inc. ISBN 0127056645.

J J Hopfield. Neural networks and physical systems with emergent collective computational abilities. *Proceedings of the National Academy of Sciences*, 79(8):2554–2558, 1982. doi: 10.1073/pnas.79. 8.2554.

Alexander Ihler and David McAllester. Particle belief propagation. In David van Dyk and Max Welling, editors, *Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics*, volume 5 of *Proceedings of Machine Learning Research*, pages 256–263, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16–18 Apr 2009. PMLR.

Alexander T Ihler. Accuracy Bounds for Belief Propagation. In *Proceedings of the Twenty-Third Conference on Uncertainty in Artificial Intelligence*, UAI'07, pages 183–190, Arlington, Virginia, USA, 2007. AUAI Press. ISBN 0974903930.

Ernst Ising. Beitrag zur Theorie des Ferromagnetismus. *Zeitschrift fur Physik*, 31(1):253–258, 1925. doi: 10.1007/BF02980577.

Renaud Jardri and Sophie Denève. Circular inferences in schizophrenia. *Brain*, 136(11):3227–3241, 2013. ISSN 14602156. doi: 10.1093/brain/awt257.

Renaud Jardri, Sandrine Duverne, Alexandra S. Litvinova, and Sophie Denève. Experimental evidence for circular inference in schizophrenia. *Nature Communications*, 8:14218, jan 2017. ISSN 2041-1723. doi: 10.1038/ncomms14218.

Christian Knoll and Franz Pernkopf. On Loopy Belief Propagation - Local Stability Analysis for Non-Vanishing Fields. In *Proceedings of the conference on Uncertainty in Artificial Intelligence*, 2017.

D. Koller and N. Friedman. *Probabilistic Graphical Models: Principles and Techniques*. Adaptive computation and machine learning. MIT Press, 2009. ISBN 9780262013192.

F.R. Kschischang, B.J. Frey, and H.-A. Loeliger. Factor graphs and the sum-product algorithm. *IEEE Transactions on Information Theory*, 47(2):498–519, 2001. doi: 10.1109/18.910572.

Jonathan Kuck, Shuvam Chakraborty, Hao Tang, Rachel Luo, Jiaming Song, Ashish Sabharwal, and Stefano Ermon. Belief propagation neural networks. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, *Advances in Neural Information Processing Systems*, volume 33, pages 667–678. Curran Associates, Inc., 2020.

Miguel Lazaro-Gredilla, Antoine Dedieu, and Dileep George. Perturb-and-max-product: Sampling and learning in discrete energy-based models. In M Ranzato, A Beygelzimer, Y Dauphin, P S Liang, and J Wortman Vaughan, editors, *Advances in Neural Information Processing Systems*, volume 34, pages 928–940. Curran Associates, Inc., 2021.

Pantelis Leptourgos, Sophie Denève, and Renaud Jardri. Can circular inference relate the neuropathological and behavioral aspects of schizophrenia? *Current Opinion in Neurobiology*, 46:154–161, oct 2017. ISSN 0959-4388. doi: 10.1016/J.CONB.2017.08.012.

Pantelis Leptourgos, Vincent Bouttier, Sophie Denève, and Renaud Jardri. From hallucinations to synaesthesia: a circular inference account of unimodal and multimodal erroneous percepts in clinical and drug-induced psychosis. *PsyArXiv*, 2021.
Shai Litvak, Leonid Karlinsky, and Shimon Ullman. Properties of Cortical Networks Improve Inference in Highly Interconnected Graphical Models. 2009.

Dong Liu, Nima N. Moghadam, Lars K. Rasmussen, Jinliang Huang, and Saikat Chatterjee. α Belief Propagation as Fully Factorized Approximation. In *2019 IEEE Global Conference on Signal and Information Processing (GlobalSIP)*, pages 1–5, 2019. doi: 10.1109/GlobalSIP45357.2019.8969545.

Dong Liu, Minh Thành Vu, Zuxing Li, and Lars K. Rasmussen. α Belief propagation for approximate inference. *arXiv, 2020*. ISSN 23318422.

Xiangfeng Liu and Ying Li. Deep mimo detection based on belief propagation. In *2018 IEEE Information Theory Workshop (ITW)*, pages 1–5, 2018. doi: 10.1109/ITW.2018.8613336.

Po-Ling Loh and Andre Wibisono. Concavity of reweighted Kikuchi approximation. In Z Ghahramani, M Welling, C Cortes, N Lawrence, and K Q Weinberger, editors, *Advances in Neural Information Processing Systems*, volume 27. Curran Associates, Inc., 2014.

R.J. McEliece, D.J.C. MacKay, and Jung-Fu Cheng. Turbo decoding as an instance of pearl's "belief propagation" algorithm. *IEEE Journal on Selected Areas in Communications*, 16(2):140–152, 1998. doi: 10.1109/49.661103.

Thomas Minka. Power EP. Technical report, 2004.

Thomas Minka. Divergence measures and message passing. Technical report, 2005.

Thomas Minka and John Lafferty. Expectation-propagation for the generative aspect model. In *Proceedings of the Eighteenth Conference on Uncertainty in Artificial Intelligence*, UAI'02, page 352–359, San Francisco, CA, USA, 2002. Morgan Kaufmann Publishers Inc. ISBN 1558608974.

Thomas P. Minka. Expectation propagation for approximate Bayesian inference. In *Uncertainty in Artificial Intelligence*, volume 17, pages 362–369, 2001.

Gianluigi Mongillo and Sophie Deneve. Online Learning with Hidden Markov Models. *Neural Computation*, 20(7):1706–1716, 07 2008. ISSN 0899-7667. doi: 10.1162/neco.2008.10-06-351.

Joris M. Mooij. libDAI: A free and open source C++ library for discrete approximate inference in graphical models. *Journal of Machine Learning Research*, 11:2169–2173, August 2010.

Joris M Mooij and Hilbert Kappen. Bounds on marginal probability distributions. In D Koller, D Schuurmans, Y Bengio, and L Bottou, editors, *Advances in Neural Information Processing Systems*, volume 21. Curran Associates, Inc., 2009.

Joris M. Mooij and Hilbert J. Kappen. Sufficient conditions for convergence of the sum-product algorithm. *IEEE Transactions on Information Theory*, 53(12):4422–4437, 2007a. doi: 10.1109/TIT.2007.909166.

Joris M. Mooij and Hilbert J. Kappen. Loop corrections for approximate inference on factor graphs. *J. Mach. Learn. Res.*, 8:1113–1143, December 2007b. ISSN 1532-4435.

Joris Marten Mooij and Hilbert J Kappen. Validity estimates for loopy Belief Propagation on binary real-world networks. In *Advances in neural information processing systems*, 2004.

Kevin P. Murphy, Yair Weiss, and Michael I. Jordan. Loopy belief propagation for approximate inference: An empirical study. In *Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence*, UAI'99, page 467–475, San Francisco, CA, USA, 1999. Morgan Kaufmann Publishers Inc. ISBN 1558606149.

Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, *Advances in Neural Information Processing Systems 32*, pages 8024–8035. Curran Associates, Inc., 2019.
Judea Pearl. *Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference*. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988. ISBN 0934613737.

L. Personnaz, I. Guyon, and G. Dreyfus. Collective computational properties of neural networks: New learning mechanisms. *Phys. Rev. A*, 34:4217–4228, 1986. doi: 10.1103/PhysRevA.34.4217.

C. Peterson and J. R. Anderson. A mean field theory learning algorithm for neural networks. *Complex Systems*, 1:995–1019, 1987.

Rajesh P.N. Rao and Dana H. Ballard. Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects. *Nature Neuroscience*, 2(1):79–87, 1999. ISSN 1097-6256. doi: 10.1038/4580.

M. Riedmiller and H. Braun. A direct adaptive method for faster backpropagation learning: the rprop algorithm. In *IEEE International Conference on Neural Networks*, pages 586–591 vol.1, 1993. doi: 10.1109/ICNN.1993.298623.

Tanya G. Roosta, Martin J. Wainwright, and Shankar S. Sastry. Convergence analysis of reweighted sum-product algorithms. *IEEE Transactions on Signal Processing*, 56(9):4293–4305, 2008. doi: 10.1109/TSP.2008.924136.

Xiangqiong Shi, Dan Schonfeld, and Daniela Tuninetti. Message error analysis of loopy belief propagation for the sum-product algorithm. *CoRR*, abs/1009.2305, 2010.

Aaron P. Shon and Rajesh P.N. Rao. Implementing belief propagation in neural circuits. *Neurocomputing*, 65-66:393–399, 2005. ISSN 09252312. doi: 10.1016/j.neucom.2004.10.035.

Pritam Som, Tanumay Datta, A. Chockalingam, and B. Sundar Rajan. Improved large-mimo detection based on damped belief propagation. In *2010 IEEE Information Theory Workshop on Information Theory (ITW 2010, Cairo)*, pages 1–5, 2010. doi: 10.1109/ITWKSPS.2010.5503188.

Andreas Steimer, Wolfgang Maass, and Rodney Douglas. Belief propagation in networks of spiking neurons. *Neural Computation*, 21(219), 2009.

Erik B. Sudderth, Alexander T. Ihler, William T. Freeman, and Alan S. Willsky. Nonparametric belief propagation. In *Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition*, volume 1, 2003. doi: 10.1109/cvpr.2003.1211409.

Nobuyuki Taga and Shigeru Mase. Error bounds between marginal probabilities and beliefs of loopy belief propagation algorithm. In *Proceedings of the 5th Mexican International Conference on Artificial Intelligence*, MICAI'06, page 186–196, Berlin, Heidelberg, 2006. Springer-Verlag. ISBN 3540490264. doi: 10.1007/11925231\_18.

Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261–272, 2020. doi: 10.1038/s41592-019-0686-2.

M. J. Wainwright, T. S. Jaakkola, and A. S. Willsky. A new class of upper bounds on the log partition function. *IEEE Trans. Inf. Theor.*, 51(7):2313–2335, 2005. ISSN 0018-9448. doi: 10.1109/TIT.2005.850091.

Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. *Foundations and Trends in Machine Learning*, 1(1–2):1–305, 2008. ISSN 1935-8237. doi: 10.1561/2200000001.

Martin J. Wainwright, Tommi Jaakkola, and Alan S. Willsky. Tree-based reparameterization for approximate inference on loopy graphs. *Advances in Neural Information Processing Systems*, 2002. ISSN 10495258.
Martin J. Wainwright., Tommi Jaakkola, and Alan S. Willsky. Tree-reweighted belief propagation algorithms and approximate ML estimation by pseudo-moment matching. Workshop on Artificial Intelligence and Statistics, 21:97, 2003.

Martin J. Wainwright, Tommi Jaakkola, and Alan S. Willsky. Tree-based reparameterization framework for analysis of sum-product and related algorithms. IEEE Transactions on Information Theory, 49(5):1120–1146, 2003.

Yusuke Watanabe. Uniqueness of Belief Propagation on Signed Graphs. In Advances in neural information processing systems, 2011.

Yusuke Watanabe and Kenji Fukumizu. Graph zeta function in the bethe free energy and loopy belief propagation. In Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems, volume 22. Curran Associates, Inc., 2009.

Yair Weiss. Correctness of Local Probability Propagation in Graphical Models with Loops. Neural Computation, 2000.

Yair Weiss. Comparing the mean field method and belief propagation for approximate inference in mrfs. In Manfred Opper and David Saad, editors, Advanced Mean Field Methods: Theory and Practice, pages 229–239. MIT Press, 2001.

Adrian Weller, Kui Tang, David Sontag, and Tony Jebara. Understanding the bethe approximation: When and how can it go wrong? In Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI’14, page 868–877, Arlington, Virginia, USA, 2014. AUAI Press. ISBN 9780974903910.

Wim Wiegerinck and Tom Heskes. Fractional Belief Propagation. In Advances in Neural Information Processing Systems, volume 15, 2002.

Hugh R. Wilson and Jack D. Cowan. Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons. Biophysical Journal, 12(1):1, 1972. doi: 10.1016/S0006-3495(72)86068-5.

John Winn and Christopher M. Bishop. Variational message passing. Journal of Machine Learning Research, 6:661–694, 2005.

Junmei Yang, Chuan Zhang, Xiao Liang, Shugong Xu, and Xiaohu You. Improved symbol-based belief propagation detection for large-scale mimo. In 2015 IEEE Workshop on Signal Processing Systems (SiPS), pages 1–6, 2015. doi: 10.1109/SiPS.2015.7345035.

Jonathan S Yedidia, William T Freeman, and Yair Weiss. Generalized Belief Propagation. Advances in neural information processing systems, 2001.

Jonathan S. Yedidia, William T. Freeman, and Yair Weiss. Understanding Belief Propagation and Its Generalizations, page 239–269. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 2003. ISBN 1558608117.

Jonathan S. Yedidia, William T. Freeman, and Yair Weiss. Constructing Free Energy Approximations and Generalized Belief Propagation Algorithms. IEEE Transactions on Information Theory, 51(7): 2282–2312, 2005.

KiJung Yoon, Renjie Liao, Yuwen Xiong, Lisa Zhang, Ethan Fetaya, Raquel Urtasun, Richard Zemel, and Xaq Pitkow. Inference in probabilistic graphical models by graph neural networks. In ICLR Workshop, 2018.
# Appendices

## Table of Contents

<table><tr><td>A</td><td>BP and Circular BP: additional figures</td><td>16</td></tr><tr><td>B</td><td>Supervised learning of Circular BP - Additional analyses</td><td>17</td></tr><tr><td>B.1</td><td>Influence of the interaction strength</td><td>17</td></tr><tr><td>B.2</td><td>Generalization ability</td><td>17</td></tr><tr><td>B.3</td><td>Learning CBP on structured graphs</td><td>19</td></tr><tr><td>B.4</td><td>Learning CBP with positive couplings</td><td>19</td></tr><tr><td>B.5</td><td>Learning CBP in larger graphs</td><td>19</td></tr><tr><td>B.6</td><td>Breaking symmetry</td><td>20</td></tr><tr><td>B.7</td><td>Comparison between all algorithms</td><td>21</td></tr><tr><td>C</td><td>Unsupervised learning of Circular BP - Additional details</td><td>23</td></tr><tr><td>C.1</td><td>Justification for the unsupervised learning rule on α and κ</td><td>23</td></tr><tr><td>C.2</td><td>Unsupervised learning of β and γ</td><td>24</td></tr><tr><td>D</td><td>Convergence properties of Circular BP - Full proofs</td><td>25</td></tr><tr><td>D.1</td><td>Convergence results for Circular BP</td><td>25</td></tr><tr><td>D.2</td><td>Extension of the convergence results to other related algorithms</td><td>27</td></tr><tr><td>E</td><td>Theoretical background for Circular BP</td><td>28</td></tr><tr><td>E.1</td><td>Reweighted BP algorithms</td><td>28</td></tr><tr><td>E.2</td><td>Reweighted BP algorithms</td><td>32</td></tr><tr><td>E.3</td><td>Mathematical definitions of several Reweighted BP algorithms</td><td>33</td></tr><tr><td>E.4</td><td>Definitions</td><td>33</td></tr><tr><td>E.5</td><td>Relations between Reweighted BP algorithms</td><td>34</td></tr><tr><td>E.6</td><td>Circular BP as approximation of Reweighted BP</td><td>35</td></tr><tr><td>F</td><td>Damped Circular BP</td><td>39</td></tr><tr><td>G</td><td>Circular BP in general factor graphs</td><td>40</td></tr><tr><td>G.1</td><td>Inference in graphical models</td><td>40</td></tr><tr><td>G.2</td><td>Belief Propagation</td><td>40</td></tr><tr><td>G.3</td><td>Reweighted BP</td><td>41</td></tr><tr><td>G.4</td><td>Theoretical background for Reweighted BP</td><td>42</td></tr><tr><td>G.5</td><td>From Gibbs free energy approximation to messages</td><td>43</td></tr><tr><td>G.6</td><td>Comparison with related models: BP, Fractional BP, Power EP, alpha-BP, and Circular BP</td><td>45</td></tr></table>
A BP and Circular BP: additional figures

We illustrate in Figure S1 the Circular Belief Propagation algorithm applied on an Ising model, in the special case where (κ, β, γ) = (1, 1, 1), that is, with parameter α only. In this case, the Circular BP algorithms writes (special cases of Equations (6), (7), and (8)):

$$
M_{i \to j}^{\text{new}} = f(B_i - \alpha_{ij} M_{j \to i}, J_{ij}) \quad (16)
$$

$$
B_i = \sum_{j \in \mathcal{N}(i)} M_{j \to i} + M_{\text{ext} \to i} \quad (17)
$$

where

$$
f(x, J) = \operatorname{arctanh} \left( \tanh(J) \tanh(x) \right) \qquad (18)
$$

The Belief Propagation algorithm consists of setting the last parameter to the value 1: α = 1.

A

<figure><img src="image_6.png" /></figure>

B

<figure><img src="image_7.png" /></figure>

C

<figure><img src="image_8.png" /></figure>

D

<figure><img src="image_9.png" /></figure>

Figure S1: Running the Belief Propagation (BP) algorithm and its variant, the Circular Belief Propagation (CBP) algorithm. (A) Example of acyclic graph taken for the simulation. In the example, the probability distribution corresponds to an Ising model: pairwise potentials ψ$_{ij}$ ( x$_{i}$, x$_{j}$ ) ∝ exp( J$_{ij}$x$_{i}$x$_{j}$ ) and unitary potentials ψ$_{i}$ ( x$_{i}$ ) ∝ exp( M$_{ext}$$_{→}$$_{i}$x$_{i}$ ) where x$_{i}$ ∈ {− 1; +1 } (binary case). J$_{ij}$ is generated randomly ( ∼ N (0 , 3)), as well as M$_{ext}$$_{→}$$_{i}$ ( ∼ N (0 , 2)). (B) The update function f for both BP and CBP is a parametric sigmoidal function close to the hyperbolic tangent. The parameter J$_{ij}$ represents the level of trust between variables x$_{i}$ and x$_{j}$ . (C) Belief Propagation is a message-passing algorithm which consists of running the update equation (2) until convergence of the messages. The approximate marginals, or beliefs, are defined by Equation (3). Here the beliefs found by BP are exact as the graph has no cycles. (D) The Circular Belief Propagation algorithm is a parametric generalization of the Belief Propagation algorithm with parameter α$_{ij}$ assigned to each edge ( i, j ) of the graph. It is identical to BP for α = 1. In the simulation, α is taken uniformly over the edges, equal to 0.5.
# B Supervised learning of Circular BP - Additional analyses

In this section, we show additional analyses related to supervised learning of CBP parameters.

The supervised learning method for learning parameters of Circular BP $(\alpha, \kappa, \beta, \gamma)$ is described in the main text. Additionally, many other models than CBP are fitted, using the same learning procedure. This includes special cases of CBP (with one or more of the parameters being fixed) but different classes of algorithms as well, like Fractional BP, the classical tanh model, etc. Each model is initialized when possible with parameters which guarantee convergence of the algorithm (otherwise the learning procedure sometimes does not find a region of convergence, most of the time because there is none but not always); see Appendix D for results on the convergence of CBP depending on its parameters. Convergence can be ensured through the combination of low and positive $\alpha$ and $\kappa$ as for CBP, but also with low enough $\beta$ and lastly, with low enough $\kappa$ (even if there is no theoretical result about the spectral radius of matrix A being lower than one if $\kappa$ is sufficiently low, it is the case in practice with the generated graphs and thus guarantees convergence).

## B.1 Influence of the interaction strength

In the numerical experiments presented in the main text, $J_{ij} \sim N(0, 1)$. Figure S2 shows the effect of increasing the strength of interactions weights $J = \{J_{ij}\}$ from a given weighted graph. BP clearly shows a performance decrease with increased interaction strength, while CBP remains very good qualitatively.

<figure><img src="image_10.png" /></figure>

Figure S2: **Approximate marginals versus true marginals on the test set, for Belief Propagation and CBP.** We consider a given graph randomly generated (shown on the left) with normally distributed weights $J = \{J_{ij}\}$. By increasing the strength of the graph weights, BP gets worse, but CBP still performs inference with very good quality. In all cases, CBP outperforms BP. BP shows here overconfidence over the graph, while CBP is on average not overconfident nor underconfident.

## B.2 Generalization ability

Here we show that the Circular BP algorithm, using the learning procedure described in the main text, is able to generalize to new data. This shows the goodness of the learning procedure.

**Within-set and out-of-set generalization** A first necessary check for the supervised learning procedure is to make sure that the proposed model is able to learn properly the training data and to generalize to new inputs.

The training indeed learns to represent the training data and allows the model to generalize well to unseen inputs with identical statistics as in the training set; see Figure S3A. In other words, the Circular BP algorithm is able to generalize to unseen situations, where a “situation” is a vector $M_{\text{ext}}$ of size $n_{\text{nodes}}$. *Unseen* refers to data not included in the training examples. Note that only the external inputs $M_{\text{ext}}$ change in the test set: the graph weights $J$ are fixed. In other words, parameters $\alpha, \kappa, \beta$, and $\gamma$ are learnt given the weights $J$, although for any external inputs $M_{\text{ext}}$.

The model also generalizes well to out-of-set inputs, that is, to external inputs $M_{\text{ext}}$ with different statistics from the ones used for training; see Figure S3B.

**Ablation study: Circular BP outperforms its special cases** The Circular BP algorithm is expected to outperform all its special cases, including Circular BP from Jardri and Denève [2013] (i.e., $(\kappa, \beta, \gamma) = (1, 1, 1)$) and BP (i.e., $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$). This is indeed the case, as shown in
<figure><img src="image_11.png" /><figcaption>Figure S3: Circular BP generalizes to new data after training. The topology considered here is the Erdos-Renyi model with $p = 0.6$ and 9 nodes, with 30 random weighted graphs generated, each with 100 test examples. (A) The model learns and generalizes well to the test set (within-set generalization). (B) Generalization to examples with different statistics (out-of-set generalization). While $\sigma(M_{\text{ext}}) = 1$ on the training set, the model still performs well for lower and higher standard deviations of the input. The performance goes down for increased $\sigma(M_{\text{ext}})$ as expected (the system becomes highly non-linear, and the correction brought by CBP is linear) but only slightly. Inferences remain relatively good on such examples for instance compared to BP. (C, D) Same as above but for BP, which shows frustration in (at least) some of the randomly weighted graphs, that is, absence of convergence of the algorithm, in which case the beliefs have very little to do with the correct marginals.</figcaption></figure>

<caption>Figure S4 where models are ordered according to the following score (average over graphs of the log of the MSE averaged over examples) written "$-\log_{10}(\text{avg MSE})$" in figures:</caption>

$$ score = -\frac{1}{N_{\text{graphs}}} \sum_{\text{graphs}} \log_{10} \left( \frac{1}{n_{\text{examples}}} \sum_{\text{examples}} \text{MSE}_{\text{graph, example}} \right) \quad (19) $$

where the MSE loss is defined in the main text. The performance of all models decreases with the complexity of graphs. Several models like BP perform very poorly for dense graphs because the system becomes frustrated (see Figure 2). It is the true for Circular BP as well, which shows that parameter $\alpha$ alone is not enough to make the algorithm converge. On the contrary, other models (all models outperforming or equal to "CBP nodal") show no sign of frustration thanks to additional parameters, as expected by Theorem 5.3. This allows all these models, including CBP, to keep a rather good level of performance in the approximate inference task even for complete graphs.

More generally, and quite logically, models which were special cases of others performed comparatively worse on the test set. We used the libDAI [Mooij, 2010] implementation of Generalized BP ("GBP_MIN" algorithm, i.e., using minimal clusters: one outer region for each maximal factor), Tree-Reweighted BP (with 10000 sampled trees) and Variational message-passing, to compute the associated approximate marginals.

Among all tested algorithms, it is noteworthy that neither Circular BP from Jardri and Denève [2013] nor Fractional BP, for which the only only parameter fitted is $\alpha$, provides good results for moderately dense or highly dense graphs. By looking at the marginals produced by these algorithms for Erdos-Renyi graphs (see Figure 2, the reason why they do not perform well is the frustration of the system caused by cycles (strong oscillations between the two extremes $b_i(x_i) \approx 0$ and $b_i(x_i) \approx 1$ without convergence of the algorithm). As observed by Murphy et al. [1999], convergence of BP implies a good approximation of the correct marginals by the beliefs, and convergence can be forced by using a damped algorithm. Similarly, here, when Circular and Fractional BP converge (for graph with low density), they produce beliefs close to the correct marginals. However, for denser graphs, the algorithms do not converge and in this case beliefs have very little to do with the correct marginals, which Murphy pointed as well. As in Murphy et al. [1999], introducing damping to CBP and FBP
<figure><img src="image_12.png" /><figcaption>Figure S4: **From BP to Circular BP: additional parameters help generalize better.** Circular BP outperforms its special cases: indeed, adding parameters to the algorithms from BP (no parameter) to Circular BP (parameters α, κ, β, γ) increasingly helps generalize better to the test set. More generally, all models which are special cases of others perform worse comparatively on the test set. This indicates an absence of overfitting: parameters of CBP can be learnt despite their consequent number w.r.t. the amount of training data. Models are ordered w.r.t. their performance on Erdos-Renyi graphs with *p* = 0.6. Each point represents the log-MSE score on the test set, averaged over 30 weighted graphs (where the MSE is averaged over examples). Weighted graphs are randomly generated from a graph topology, with normally generated weights *J*<sub>ij</sub> ~ N(0, 1). Each weighted graph has 100 test examples.</figcaption></figure>

would help the algorithms to converge. However, the damped algorithm would become slower to
converge in cases where the undamped algorithm converged. This does not seem like an optimal
solution to the frustration problem, as the inference performed by the algorithm should be as fast
as possible (including when the inputs vary with time). The question whether the algorithms can
converge only by choosing the right α, without damping, has not been answered yet. It is probable
that the α parameter is not enough to prevent the frustration behavior from occurring. Alternately, it
is possible that the learning procedure does not find the optimal α. This second possibility as yet not
been ruled out.

B.3 Learning CBP on structured graphs

While the main text presents results on Erdos-Renyi graphs, Figures S5 and S6 shows that the
conclusions remain the same for structured graphs.

B.4 Learning CBP with positive couplings

We take here homogeneous positive couplings $J_{ij}$ and external fields $M_{ext \to i}$, a model designed to show off the weakness of BP; see also [Yedidia et al., 2001, Roosta et al., 2008]. For simplicity, we possibly in d-regular graphs like grid with toroidal boundary conditions). See Figure S7 for the results.

B.5 Learning CBP in larger graphs

Figure S8 shows that CBP (and its supervised optimization procedure on CBP parameters) scales up on graphs of moderate size, on which it is still possible to compute the true marginals. We took the example of 10 * 10 grid with normally distributed weights $J_{ij} \sim N(0, 1)$ (spin-glass). This contrasts with the main text, which only considers graphs with 9 nodes. Note also that the Hopfield network application from the main text involves unsupervised learning in a all-to-all connected graph with 784 nodes.
<figure><img src="image_13.png" /><figcaption>Figure S5: Estimated marginals of Circular BP and BP for graphs of various topologies, on the test set. For each graph topology, 30 randomly weighted graphs are considered with $J_{ij} \sim N(0, 1)$, each with 200 external input examples in the training set and 100 in the test set. Results are qualitatively the same as in Erdos-Renyi graphs. Graph structures are the same as in Yoon et al. [2018].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure S6: Quantitative comparison between various algorithms on different graph structures. The results are qualitatively the same as for randomly connected graphs.</figcaption></figure>

## B.6 Breaking symmetry

In Circular BP, parameters $\alpha = \{\alpha_{ij}\}$ and $\beta = \{\beta_{ij}\}$ are symmetric matrices as elements of the matrix depend on the *unoriented* edge (for instance, for $\alpha$, $\alpha_{i \to j} = \alpha_{j \to i}$). We tested whether removing this symmetry constraint would improve the inference algorithm. In this case, the message update equation of Circular BP becomes:

$$m_{i \to j}^{\text{new}}(x_j) \propto \sum_{x_i} \psi_{ij}(x_i, x_j)^{\beta_{i \to j}} (\psi_i(x_i)^{\gamma_i} \prod_{k \in N(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{i \to j}/\kappa_i})^{\kappa_i} \quad (20)$$

instead of Equation (4). Crucially, $\alpha_{i \to j}$ and $\beta_{i \to j}$ are assigned to a *directed* edge $(i, j)$. A consequence is that quantities $\alpha_{i \to j}$, $\alpha_{j \to i}$, $\beta_{i \to j}$ and $\beta_{j \to i}$ do not represent inverse overcounting numbers anymore, as these parameters weigh quantities which depend only on the pair $(i, j)$ but not the direction.

As Figure S9 shows, removing the symmetry constraint on $\alpha$ and $\beta$ does not significantly allows for better approximate inference, even though the newly created algorithm ("CBP + non-symmetric" in the figure) is a generalization of the previous one (CBP). The number of parameters to fit increases ($n_{edges}$ additional parameters for $\alpha$ and $n_{edges}$ additional parameters for $\beta$). The figure suggests a slight overfitting, together with the fact that the performance of "CBP + non-symmetric" is systematically (slightly) better than the one of CBP on the training set (not shown). An explanation among others could be the low amount of training data.
<figure><img src="image_15.png" /><figcaption>Figure S7: **Learning CBP with positive couplings.** Here couplings are taken homogeneous with $J_{ij} = \operatorname{arctanh}(0.6) \approx 0.7$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure S8: **Learning CBP in larger graphs.** Example of a grid with 100 nodes, on which supervised learning was applied. 4 different weighted graphs were generated.</figcaption></figure>

## B.7 Comparison between all algorithms

To complete the figures, which only illustrate specific points from the text by showing algorithms of the same class, we provide in Table 1 a ranked list containing algorithms of all types. The rank is defined by the performance of each algorithm for Erdos-Renyi graphs with connection probability $p = 0.6$. For each connection probability, 30 randomly weighted graphs of 9 nodes were generated, each with 100 examples (vectors $M_{\text{ext}}$) in the test set. More precisely, the performance is measured identically to previously, i.e., by computing the average over graphs of $-\log_{10}(\text{avg MSE})$ between the true marginals and the estimated marginals where “avg” of “avg MSE” is taken over examples.

Overall, among fitted models, the best algorithm was Reweighted BP, followed by Circular BP. Both significantly outperformed BP. Loop-Corrected BP is the best algorithm overall, although much more complex than Reweighted BP or Circular BP.

Figure 2 from the main text provides visual comparison for Erdos-Renyi graphs (including the case $p = 0.6$ considered in the above table) for most models listed in Table 1.
<figure><img src="image_17.png" /><figcaption>Figure S9: Breaking the symmetry constraint in the parameter matrices α and β does not provide significant improvement to the model, although it brings 2nedges new parameters to fit (it even performs worse sometimes, because of slight overfitting). However, it frees some constraints on the parameters. For clarity, several points are not shown, as all models perform extremely well for trees and graphs close to trees (here p = 0.2 and “cycle”).</figcaption></figure>

<table><caption>Table 1: Comparison between various classes of algorithms. Graphs considered here are Erdos-Renyi graphs (p = 0.6). The score measure is defined in Equation (19).</caption><thead><tr><th>Name</th><th>Score</th><th># parameters</th><th>Reference</th></tr></thead><tbody><tr><td>Loop-corrected BP</td><td>6.38</td><td>-</td><td>Mooij and Kappen [2007b]</td></tr><tr><td>Circular BP</td><td>4.83</td><td>2n<sub>edges</sub> + 2n<sub>nodes</sub></td><td>this paper</td></tr><tr><td>Classical tanh model</td><td>2.21</td><td>2n<sub>edges</sub> + n<sub>nodes</sub></td><td>Wilson and Cowan [1972]</td></tr><tr><td>Fractional BP</td><td>2.07</td><td>n<sub>edges</sub></td><td>Wiegerinck and Heskes [2002]</td></tr><tr><td>Double-loop BP</td><td>1.92</td><td>-</td><td>Heskes et al. [2002]</td></tr><tr><td>Circular BP with (κ, β, γ) = (1, 1, 1)</td><td>1.88</td><td>n<sub>edges</sub></td><td>Jardri and Denève [2013]</td></tr><tr><td>BP</td><td>1.49</td><td>-</td><td>Pearl [1988]</td></tr><tr><td>Tree-Reweighted BP</td><td>1.22</td><td>-</td><td>Wainwright. et al. [2003]</td></tr><tr><td>Mean Field / Variational MP</td><td>0.97</td><td>-</td><td>Winn and Bishop [2005]</td></tr></tbody></table>
## C Unsupervised learning of Circular BP - Additional details

In this section, we provide additional considerations related to the unsupervised learning rule proposed in the main text. Section C.1 proposes a mathematical intuition behind the learning rule on $\alpha$ and $\kappa$ proposed in the main text. Section C.2 proposes a way of learning $\beta$ and $\gamma$ in an unsupervised fashion (contrary to the main paper which investigates unsupervised learning of the other two parameters of CBP: $\alpha$ and $\kappa$).

### C.1 Justification for the unsupervised learning rule on $\alpha$ and $\kappa$

**Learning rule on $\alpha$** The unsupervised learning rule for $\alpha$ minimizes the total amount of information sent in the network $\sum_{(i,j)} M_{i \to j}^2$, or equivalently, minimizes the quantity of so-called prediction errors

$$ \sum_{(i,j)} (B_i - B_i^j)^2: $$

$$ L = \sum_{(i,j)} M_{i \to j}^2 \qquad (21) $$

$$ \approx \sum_{(i,j)} \left( \operatorname{arctanh} \left[ \tanh(J_{ij}) \tanh(B_i - \alpha_{ij} M_{j \to i}) \right] \right)^2 \qquad (22) $$

We compute the derivative of the loss function $L$ w.r.t. $\alpha$ by considering that $\frac{\partial(B_i - \alpha_{ij} M_{j \to i})}{\partial \alpha_{kl}} \approx -M_{j \to i} \delta_{jk} \delta_{il}$, i.e., by neglecting higher-order dependencies. Consequently,

$$ \frac{\partial L}{\partial \alpha_{ij}} \approx \frac{\partial M_{i \to j}^2}{\partial \alpha_{ij}} + \frac{\partial M_{j \to i}^2}{\partial \alpha_{ij}} \qquad (23) $$

where

$$ \frac{\partial M_{i \to j}^2}{\partial \alpha_{ij}} \approx -2M_{i \to j}(\operatorname{arctanh})'(\tanh(J_{ij})\tanh(B_i - \alpha_{ij}M_{j \to i})) \times \tanh(J_{ij})\tanh'(B_i - \alpha_{ij}M_{j \to i})M_{j \to i} \qquad (24) $$

which has the same sign as $-2M_{i \to j} \tanh(J_{ij})M_{j \to i}$ (symmetrical in $(i, j)$). We thus propose the following learning rule:

$$ \Delta \alpha_{ij} \propto -\frac{\partial L}{\partial \alpha_{ij}} \propto M_{i \to j} \tanh(J_{ij}) M_{j \to i} \qquad (25) $$

This learning rule tends to make the correlations between opposite messages disappear, i.e., $M_{i \to j} M_{j \to i} >_{\text{examples}} \approx 0$. In the simulations, we used an alternative learning rule (see Equation (14)):

$$ \Delta \alpha_{ij} \propto (B_i - \alpha_{ij} M_{j \to i}) M_{j \to i} + (B_j - \alpha_{ij} M_{i \to j}) M_{i \to j} \qquad (26) $$

(the right-hand terms of Equations (25) and (26) have identical signs).

**Learning rule on $\kappa$** The learning rule on $\kappa$ aims at avoiding reverberation of external information by cycles. $\kappa_i$ acts as a scaling factor on the log-likelihood ratio (see Equation (7)). The proposed learning rule aims at decorrelating the information truly received by node $i$ from the outside world ($M_{\text{ext} \to i}$), and all the information received by node $i$ except this true external information ($B_i - M_{\text{ext} \to i}$):

$$ \Delta \kappa_i \propto -M_{\text{ext} \to i}(B_i - M_{\text{ext} \to i}) \qquad (27) $$

Note that a learning rule consisting in from minimizing, as for $\alpha$, the quantity of messages sent throughout the graph, would not work here as $\kappa \to 0$ would be enough to have all messages go to zero. Finally, the learning rule $\Delta \kappa_i \propto -M_{\text{ext} \to i}(B_i - \kappa_i M_{\text{ext} \to i})$ would not work either as the sign of the correlation between $M_{\text{ext} \to i}$ and $B_i - \kappa_i M_{\text{ext} \to i} \propto M_{j \to i}$ does not depend on $\kappa_i$ but instead on the graph topology and weights $J$.
## C.2 Unsupervised learning of $\beta$ and $\gamma$

In section 6.2 of the main paper, we considered the case where $(\beta, \gamma) = (1, 1)$ in Circular BP and proposed a way of learning parameters $\alpha$ and $\kappa$ in an unsupervised fashion. Here we ask the question of the learning of the additional parameters $\beta$ and $\gamma$.

As a reminder, parameter $\beta$ simply appears in the formulation of CBP as a multiplicative factor to the true weights $J$ (defined by $\psi_{ij}(x_i, x_j) = \exp(J_{ij}x_i x_j)$). Similarly, parameter $\gamma$ simply appears as a multiplicative factor to the true external inputs $M_{\text{ext}}$ (defined by $\psi_i(x_i) = \exp(M_{\text{ext}} \rightarrow i x_i)$). Therefore, fitting $\beta$ and $\gamma$ is equivalent to fitting $J$ and $M_{\text{ext}}$, respectively the graph weights and inputs (where $J$ and $M_{\text{ext}}$ are not related anymore to the factors of the probability distribution).

Parameter $\beta$ can be learnt similarly to Mongillo and Deneve [2008], Jardri and Denève [2013] which proposes an expectation maximization (EM) algorithm [Dempster et al., 1977] and alternatively a Hebbian-like learning rule using stochastic gradient descent, to learn the factors $\{\psi_{ij}\}$. This approach is based on the fact that messages from BP and its variants allow not only to compute unitary marginal probabilities $\{p_i(x_i)\}$ but also pairwise marginal probabilities $\{p_{ij}(x_i, x_j)\}$. For instance, for CBP, one can assume that:

$$ b_{ij}(x_i, x_j) \propto \psi_{ij}(x_i, x_j)^{\alpha_{ij}\beta_{ij}} \psi_i(x_i)^{\kappa_i\gamma_i} \psi_j(x_j)^{\kappa_j\gamma_j} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \rightarrow i}(x_i)^{\kappa_i} \\ \times \prod_{k \in \mathcal{N}(j) \setminus i} m_{k \rightarrow j}(x_j)^{\kappa_j} m_{j \rightarrow i}(x_i)^{\kappa_i - \alpha_{ij}} m_{i \rightarrow j}(x_j)^{\kappa_i - \alpha_{ij}} \quad (28) $$

$$ \Rightarrow b_{ij}(x_i, x_j) \propto \psi_{ij}(x_i, x_j)^{\alpha_{ij}\beta_{ij}} \times \frac{b_i(x_i)}{m_{j \rightarrow i}(x_i)^{\alpha_{ij}}} \times \frac{b_j(x_j)}{m_{i \rightarrow j}(x_j)^{\alpha_{ij}}} \quad (29) $$

which has the same form as the expression of the message update equation for CBP: $m_{i \rightarrow j}(x_j)$ depends on $b_i(x_i)/m_{j \rightarrow i}(x_i)^{\alpha_{ij}}$ in Equation (4). In particular, for an Ising model, Equation (29) gives:

$$ \log \left( \frac{b_{ij}(1, 1)}{b_{ij}(1, 0)} \right) = 2\alpha_{ij}\beta_{ij}J_{ij} + (B_j - \alpha_{ij}M_{i \rightarrow j}) \quad (30) $$

which shows the importance of the quantity $B_j - \alpha_{ij}M_{i \rightarrow j}$, together with the message update equation of CBP $M_{j \rightarrow i} = f(B_j - \alpha_{ij}M_{i \rightarrow j}, \beta_{ij}J_{ij})$. Note, however, that Equations (28), (29) and (30) are only approximations for CBP. Indeed, Equation (28) is the expression of the pairwise beliefs for Reweighted BP and not Circular BP (see Equation (54)). Such an equation does not exist for Circular BP, as CBP does not come from any approximation of the Gibbs free energy (which is the starting point to yield such an expression). However, as seen previously, CBP can be defined as an approximation to Reweighted BP with a slightly simpler message update equation. By making the hypothesis that even though CBP has a different update equation for the messages, it eventually computes the (approximate) solution to Reweighted BP, then Equation (28) is (approximately) valid.

Jardri and Denève [2013] proposes, in the particular situation where $(\kappa, \beta, \gamma) = (1, 1, 1)$ (i.e., for the original Circular BP), an online update of the factor parameters (here, simply $J_{ij}$ or equivalently, $\beta_{ij}$) based on the pairwise beliefs $\{b_{ij}(x_i, x_j)\}$, with a Hebbian-like learning rule. Learning the external message $M_{\text{ext} \rightarrow i}$ (or equivalently, its weight $\gamma_i$) could be done similarly by defining an abstract “external unit” connected to $x_i$ through weight $\gamma_i$.

Overall, it might be possible (this still needs to be demonstrated numerically) to learn all parameters of Circular BP in an unsupervised way: $(\alpha, \kappa, \beta, \gamma)$. The learning on $\alpha$ and $\kappa$ is anti-Hebbian (or inhibitory Hebbian) and the learning on $\beta$ and $\gamma$ is Hebbian.
# D Convergence properties of Circular BP - Full proofs

In this section, we analyse the convergence properties of CBP by providing the full proofs of the theorems presented in the main text.

As seen in Figures 2 and S5, the convergence of the algorithm is crucial to carry out approximate inference: it can be observed a sharp transition between having beliefs very close to the correct marginals and very far from them. When one of the algorithms does not converge, beliefs have very little to do with the correct marginals (this observation was already made for BP in Murphy et al. [1999]). It is therefore important to control the convergence properties, and if possible, to ensure the convergence of our proposed algorithm.

Importantly, for all the variants of BP cited previously (Fractional BP, Circular BP from Jardri and Denève [2013] with $(\beta, \gamma, \kappa) = (1, 1, 1)$, Variational message-passing, etc.), there is no theoretical result about the existence of parameters that guarantee convergence. Absence of convergence of an algorithm means that it gets trapped in a limit cycle (frustration phenomenon) or alternatively wanders around chaotically without being able to produce approximate marginals. The typical behavior is a strong oscillation of the beliefs from one extreme ($b_i(x_i = 1) \approx 100\%$) to the other ($b_i(x_i = 1) \approx 0\%$). This is an serious obstacle to the use of such algorithms to perform approximate inference.

On the contrary, we prove here that in the case of Ising models, there exist parameters for which Circular BP is stable. More precisely, for any set of parameters $(\gamma, \beta)$, one can choose parameters $(\alpha, \kappa)$ such that the algorithm converges (see Theorem D.4).

Having a stable algorithm is important, for two reasons. First, it helps performing reasonably good approximate inference with the algorithm for any probability distribution. Indeed, whatever the probability distribution is (i.e., the weights $J$ and external inputs $M_{ext}$ are), it is possible to select parameters for which the algorithm will produce approximate marginals. Second, being in a stable regime allows to perform computations and ensures that estimations of probabilities at some point in time do not depend on the past state of the system (meaning that the system should have only one fixed point for any probability distribution, and the system should converge to this fixed point). In other words, the system should adapt quickly to new inputs and consider its past as little as possible.

## D.1 Convergence results for Circular BP

Here we state sufficient conditions for the convergence of Circular BP in an Ising model.

We start by defining matrix $A$ whose coefficients are:

$$A_{i \to j, k \to l} = |\kappa_i| \tanh|\beta_{ij} J_{ij}| \delta_{il} \mathbb{1}_{\mathcal{N}(i)}(k) \left|1 - \frac{\alpha_{ij}}{\kappa_i}\right|_{j=k} \quad (31)$$

where $\mathbb{1}_{\mathcal{N}(i)}(k) = 1$ if $k \in \mathcal{N}(i)$, otherwise $= 0$, $\delta_{il} = 1$ if $i = l$, otherwise $= 0$, and $\left|1 - \frac{\alpha_{ij}}{\kappa_i}\right|_{j=k} = 1 - \frac{\alpha_{ij}}{\kappa_i}$ if $j = k$, otherwise 1.

**Theorem D.1.** If for any induced operator norm $\|\cdot\|$ (sometimes called natural matrix norm), $\|A\| < 1$, then CBP converges to a unique fixed point and the rate of convergence is at least linear.

*Proof.* The proof of Theorem D.1 follows closely the one of Lemma 2 of Mooij and Kappen [2007a] for Belief Propagation, that is, the special case $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$. The message update equation for Circular BP, given by Equations (6) and (7), is:

$$M_{i \to j}^{\text{new}} = \tanh^{-1} \left[ \tanh(\beta_{ij} J_{ij}) \tanh \left( \kappa_i \gamma_i M_{\text{ext} \to i} + \kappa_i \sum_{k \in \mathcal{N}(i) \setminus j} M_{k \to i} + (\kappa_i - \alpha_{ij}) M_{j \to i} \right) \right] \quad (32)$$

Let $F$ be the update function for the messages: $M^{\text{new}} = F(M^{\text{old}})$ (in what follows, we drop the superscript old for clarity). The derivative of $F$ is:

$$F'(\mathbf{M})_{i \to j, k \to l} = \frac{\partial M_{i \to j}^{\text{new}}}{M_{k \to l}} \quad (33)$$

$$= \tilde{A}_{i \to j, k \to l} B_{i \to j}(\mathbf{M}) \quad (34)$$
where

$$ \tilde{A}_{i \to j, k \to l} = |\kappa_i| \tanh |\beta_{ij} J_{ij}| \delta_{il} \mathbb{1}_{\mathcal{N}(i)}(k) \left(1 - \frac{\alpha_{ij}}{\kappa_i}\right)_{j=k} \quad (35) $$

$$ B_{i \to j}(\mathbf{M}) = \operatorname{sgn}(\kappa_i \beta_{ij} J_{ij}) \frac{1 - \tanh^2 \left( \kappa_i \gamma_i M_{\text{ext} \to i} + \kappa_i \sum_{k \in \mathcal{N}(i) \setminus j} M_{k \to i} + (\kappa_i - \alpha_{ij}) M_{j \to i} \right)}{1 - \tanh^2 (M_{i \to j}^{\text{new}})} \quad (36) $$

Note that $\sup_M B_{i \to j}(\mathbf{M}) = 1$, as $\sup_x \frac{1 - \tanh^2(x)}{1 - \tanh^2(K) \tanh^2(x)} = 1$.

It comes that $\sup_M F'(\mathbf{M})_{i \to j, k \to l} \le A_{i \to j, k \to l}$ where we defined $A$ as $A = |\tilde{A}|$. $A$ does not depend on the external messages $M_{\text{ext}}$, nor does it depend on parameter $\gamma$.

If for any norm $\|\cdot\|$ on matrices, $\|A\| < 1$, $\sup_M F'(\mathbf{M})_{i \to j, k \to l} < 1$, then function $F$ is a $\|\cdot\|$-contraction (see Lemma 2 of Mooij and Kappen [2007a]) and the sequence $M, F(M), F \circ F(M), \dots$, that is, the Circular BP algorithm, converges to a unique fixed point with at least a linear rate. $\square$

Hence choosing in Theorem D.1 the spectral norm (induced by the $l_2$-norm), it comes straightforwardly:

**Corollary D.2.** If the largest singular value of $A$, $\sigma_{\max}(A) < 1$, then Circular BP converges.

**Theorem D.3.** If $\forall (i, j), \alpha_{ij}/\kappa_i \le 1$ and the spectral radius of $A$, $\rho(A) < 1$, then Circular BP converges to the unique fixed point.

*Proof.* The proof of Theorem D.3 follows closely the one of Corollary 3 of Mooij and Kappen [2007a]. According to Equation (34),

$$ F'(\mathbf{M})_{i \to j, k \to l} = \tilde{A}_{i \to j, k \to l} B_{i \to j}(\mathbf{M}) \quad (37) $$

Note that $A = |\tilde{A}| = \tilde{A}$ as $\tilde{A}$ is non-negative (because of the hypothesis $\alpha_{ij}/\kappa_i \le 1$). We directly conclude by using Theorem 2 of Mooij and Kappen [2007a]: $F'(\mathbf{M})$ is the product of $A$, a constant non-negative matrix, with $B(\mathbf{M})$ whose coefficients are bounded by 1 in absolute value, thus the sequence $M, F(M), F \circ F(M), \dots$, meaning, the Circular BP algorithm, converges to a fixed point and this fixed point does not depend on the initial condition $M$. $\square$

Theorem D.3 has the strong constraint that $\alpha_{ij}/\kappa_i \le 1$ for all $(i, j)$. However, when this constraint is verified, the spectral radius criterion is sharper than the norm criterion of Theorem D.1 because $\rho(A) \le \|A\|$ for all induced operator norms.

Finally, a consequence of Theorem D.3 is the following fundamental result, which distinguishes Circular BP from related approaches like Power EP, Fractional BP and $\alpha$-BP:

**Theorem D.4.** For a given weighted graph (defined by its weights $J_{ij}$), it is always possible to find parameters $\alpha$ and $\kappa$ such that Circular BP converges for any external input $M_{\text{ext}}$ and any choice of parameters $(\gamma, \beta)$.

*Proof.* Let us take $\alpha_{ij} = \kappa_i \equiv p \in \mathbb{R}_+$. In this case, $A_{i \to j, k \to l} = p \tanh |\beta_{ij} J_{ij}| \delta_{il} \mathbb{1}_{\mathcal{N}(i) \setminus j}(k)$. When $p \to 0$, all coefficients of $A$ go to zero. The spectral radius is a continuous application, and the null matrix has a spectral radius of zero, thus the spectral radius of $A$ goes to zero when $p \to 0$. We conclude by using Theorem D.3 as $\alpha_{ij}/\kappa_i = 1 \le 1$: there exists $p^*$ such that for all $p < p^*$, CBP converges. $\square$

Notably, proof of Theorem D.4 shows that choosing $\alpha$ and $\kappa$ uniformly, equal and large enough guarantees the convergence of CBP.

We use this result to initialize parameters in the supervised fitting procedure (gradient-descent based). $\alpha$ and $\kappa$ are first set at the BP value of $p = 1$, and we decrease $p$ until the spectral radius of matrix $A$ goes below 1, which ensures that Circular BP converges according to the theory. In practice, we
decrease $p$ by incrementing $1/p$ by steps of 1. Note that we are using a sufficient condition, thus Circular BP could converge for higher values of $p$ than the one chosen.

## D.2 Extension of the convergence results to other related algorithms

As shown above, Circular BP (with parameters $(\alpha, \kappa, \beta, \gamma)$) converges to the unique fixed point, whatever the probability distribution and parameters $(\beta, \gamma)$ are, given the right choice of parameters $(\alpha, \kappa)$.

In particular, for $(\beta, \gamma) = (1, 1)$, CBP with parameters $(\alpha, \kappa)$ has the same convergence properties as the general CBP. In fact, it is the combination of these two parameters that makes the convergence possible. Without parameter $\kappa$ (i.e., with $\kappa = 1$), it is not possible to guarantee that there exists $\alpha$ such that the algorithm converges. We go even further by stating the following conjecture: there exist weighted graphs (typically, with strong enough weights) for which, for any choice of parameter $\alpha$, Circular BP from Jardri and Denève [2013] (and similarly, Fractional BP, both having $\kappa = 1$) does not converge. As a reminder, parameter $\kappa_i$ scales the belief $B_i$. It is rather intuitive that it is involved in the stability of the system.

Very similarly, CBP converges to the unique fixed point, whatever the probability distribution and parameters $(\alpha, \kappa, \gamma)$ are, given the right choice of parameter $\beta$. This can be seen easily: in the demonstration of Theorem D.4, coefficients of $A$ go to zero as well if $\beta \to 0$, thus for $\beta$ uniform and sufficiently small, the algorithm converges to the unique fixed point.

As stated above, the convergence result for Circular BP distinguishes the algorithm from related approaches including Circular BP from Jardri and Denève [2013], Power EP, Fractional BP and $\alpha$-BP. Indeed, as seen above, $\alpha$ alone is not enough to ensure convergence of these algorithms: an additional parameter $\kappa$ would be needed but is not present in these algorithms. However, controlling parameter $\kappa$ alone (i.e., taking fixed $\alpha$) is not sufficient to guarantee the convergence of CBP but is enough for CBP with fixed $\alpha = 0$. Note that not being able to guarantee the convergence of an algorithm does not mean that the algorithm does not converge. Simulations seem to indicate that for strong enough weights and strong enough $\alpha$ values, $\kappa \to 0$ is not enough to make CBP converge. However, such an approach is successful on all the graphs used for which the weights are moderately strong ($J_{ij} \sim \mathcal{N}(0, 1)$).
# E Theoretical background for Circular BP

As explained in the main text, Circular BP highly relates to Reweighted BP algorithms, which are message-passing algorithms trying to minimize generalizations of the Bethe free energy (used to define the Belief Propagation algorithm). This section provides some theoretical background underlying these Reweighted BP algorithms (e.g., Fractional BP, Power EP, and α-BP), based on which the Circular BP algorithm can be defined. Importantly, we will see that the Circular BP algorithm is not associated to a particular reweighted Bethe free energy, but is very similar to algorithms trying to do so. We insist that the goal here is not to propose theoretical foundations of Circular BP per se, but instead to relate Circular BP to algorithms with solid theoretical foundations.

This section considers the special case of **Markov Random Fields with at most pairwise interactions** (and when stated, the even more special case of such distributions over binary variables). The general case follows closely the demonstration provided here, and is detailed in section G.

## E.1 Reweighted BP algorithms

### E.1.1 Gibbs free energy approximations for accurate inferences

**Problem** First, we start by writing the general variational problem: given the true probability distribution $p(\mathbf{x})$, we would like to find an approximating probability distribution $b(\mathbf{x})$ (variational distribution), whose marginals are easier to compute than those of $p(\mathbf{x})$. To do that in practice, we minimize the “difference” between $p(\mathbf{x})$ and $b(\mathbf{x})$, together with some hypotheses on $b(\mathbf{x})$ making the computation of its marginals feasible.

The difference chosen is the Kullback-Leibler divergence:

$$D_{KL}(b\|p) = \sum_{\mathbf{x}} b(\mathbf{x}) \log \left( \frac{b(\mathbf{x})}{p(\mathbf{x})} \right) \qquad (38)$$

Minimizing $D_{KL}(b\|p)$ will result in a probability distribution $b(\mathbf{x})$ “close” to the original distribution $p(\mathbf{x})$. The true probability distribution $p(\mathbf{x})$ is written $p(\mathbf{x}) = \frac{e^{-E(\mathbf{x})}}{Z}$ where Z is a normalization constant ensuring that $\sum_{\mathbf{x}} p(\mathbf{x}) = 1$.

$$\Rightarrow D_{KL}(b\|p) = \sum_{\mathbf{x}} b(\mathbf{x}) \log(b(\mathbf{x})) + \log(Z) + \sum_{\mathbf{x}} b(\mathbf{x}) E(\mathbf{x}) \qquad (39)$$

$$= -S_b - F + U_b \qquad (40)$$

where $S_b = -\sum_{\mathbf{x}} b(\mathbf{x}) \log(b(\mathbf{x}))$ is the variational entropy, $F = -\log(Z)$ is the Helmholtz free energy, and $U_b = \sum_{\mathbf{x}} b(\mathbf{x}) E(\mathbf{x})$ is the variational average energy. The **Gibbs free energy** (also called the variational free energy) is defined as $G_b = U_b - S_b$, which we want to minimize in order to minimize the KL divergence.

Given the probability distribution $p(\mathbf{x})$ and the variational distribution $b(\mathbf{x})$, the variational average energy $U_b$ can be decomposed very easily: $p(\mathbf{x}) \propto \prod_{(i,j)} \psi_{ij}(x_i, x_j) \prod_i \psi_i(x_i)$ or equivalently, $E(x) = -\sum_{(i,j)} \log(\psi_{ij}(x_i, x_j)) - \sum_i \log(\psi_i(x_i))$, and therefore,

$$\begin{aligned} U_b &= -\sum_{\mathbf{x}} b(\mathbf{x}) \sum_{(i,j)} \psi_{ij}(x_i, x_j) - \sum_{\mathbf{x}} b(\mathbf{x}) \sum_{i} \psi_i(x_i) \\ &= -\sum_{(i,j)} \sum_{x_i, x_j} b_{ij}(x_i, x_j) \psi_{ij}(x_i, x_j) - \sum_{i} \sum_{x_i} b_i(x_i) \psi_i(x_i) \end{aligned} \qquad (41)$$

However, it is most of the time very tricky to compute the entropy term $S_b = -\sum_{\mathbf{x}} b(\mathbf{x}) \log(b(\mathbf{x}))$: its formula does not simplify in the general case, and therefore involves to sum an exponential number of terms, which is infeasible in practice in graphs with a large number of nodes.

**The Bethe approximation** The Bethe approximation consists of estimating the entropy $S_b$ of $b(\mathbf{x})$ as if the probabilistic graph representing $b(\mathbf{x})$ were a tree (i.e., were cycle-free). Indeed, in the case of
trees, the probability distribution $b(\mathbf{x})$ can be written as a product of terms only involving its pairwise marginals $b_{ij}(x_i, x_j) \equiv \sum_{\mathbf{x} \setminus (x_i, x_j)} b(\mathbf{x})$ and its unitary marginals $b_i(x_i) \equiv \sum_{\mathbf{x} \setminus x_i} b(\mathbf{x})$):

$$b(\mathbf{x}) \approx \prod_{(i,j)} \frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)} \prod_i b_i(x_i) \quad (42)$$

This is equivalent to approximating the entropy $S_b$ of $b(\mathbf{x})$ as follows:

$$\begin{aligned} -S_b &\approx \sum_x b(x) \sum_{(i,j)} \log \left( \frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)} \right) + \sum_x b(x) \sum_i \log(b_i(x_i)) \\ &\approx \sum_{(i,j)} \sum_{(x_i, x_j)} b_{ij}(x_i, x_j) \log \left( \frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)} \right) + \sum_i \sum_{x_i} b_i(x_i) \log(b_i(x_i)) \end{aligned} \quad (43)$$

Note that $S_b = -\sum_{(i,j)} D_{KL}(b_{ij}, b_i b_j) + \sum_i S_{b_i}$ (which can also be written $S_b \approx \sum_{(i,j)} S_{b_{ij}} + \sum_i (1 - |\mathcal{N}(i)|)S_{b_i}$). This leads to the following approximation of the Gibbs free energy ($G_b$) known as the **Bethe approximation**: $G_b \approx G_b^{\text{Bethe}}$ where $G_b^{\text{Bethe}}$ is called the **Bethe free energy** and is defined by:

$$\begin{aligned} G_b^{\text{Bethe}} = & \sum_{(i,j)} \sum_{(x_i, x_j)} b_{ij}(x_i, x_j) \log \left( \frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)} \right) - \sum_{(i,j)} \sum_{(x_i, x_j)} b_{ij}(x_i, x_j) \log \left( \psi_{ij}(x_i, x_j) \right) \\ & + \sum_i \sum_{x_i} b_i(x_i) \log \left( b_i(x_i) \right) - \sum_i \sum_{x_i} b_i(x_i) \log \left( \psi_i(x_i) \right) \end{aligned} \quad (44)$$

**Generalizing the Bethe approximation** To allow for a better approximation of the Gibbs free energy, one can use a parametric approximation of the average energy $U_b$ and the variational entropy $S_b$, eventually leading to a generalization of the Bethe free energy. We will refer to this approximation as the **reweighted Bethe approximation**.

First, to compute the entropy term $S_b$, the variational distribution $b(x)$ is hypothesized to decompose into:

$$b(\mathbf{x}) \approx \prod_{(i,j)} \left( \frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)} \right)^{\hat{\beta}_{ij}} \prod_i (b_i(x_i))^{\hat{\gamma}_i} \quad (45)$$

with additional parameters $(\hat{\beta}$ and $\hat{\gamma})$ with respect to Equation (42). This leads to the following parametric approximation of the variational entropy $S_b$ (entropy of $b(\mathbf{x})$):

$$-S_b \approx \sum_{(i,j)} \hat{\beta}_{ij} \sum_{(x_i, x_j)} b_{ij}(x_i, x_j) \log \left( \frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)} \right) + \sum_i \hat{\gamma}_i \sum_{x_i} b_i(x_i) \log(b_i(x_i)) \quad (46)$$

Note that this is equivalent to approximating the entropy $S_b$ as $S_b \approx -\sum_{(i,j)} \hat{\beta}_{ij} D_{KL}(b_{ij}, b_i b_j) + \sum_i \hat{\gamma}_i S_{b_i}$, which can also be written as a weighted sum of local entropies $S_b \approx \sum_{(i,j)} \hat{\beta}_{ij} S_{b_{ij}} + \sum_i (\hat{\gamma}_i - \sum_{j \in \mathcal{N}(i)} \hat{\beta}_{ij}) S_{b_i}$.

Second, to compute the average energy term $U_b$, the true probability distribution $p(\mathbf{x})$ is hypothesized to decompose into:

$$p(\mathbf{x}) \approx \prod_{(i,j)} (\psi_{ij}(x_i, x_j))^{\beta_{ij}} \prod_i (\psi_i(x_i))^{\gamma_i} \quad (47)$$

with additional parameters $(\beta$ and $\gamma)$ with respect to Equation (1). Likewise, Equation (47) leads to the following parametric approximation of the average energy $U_b$:$^4$

$$U_b \approx - \sum_{(i,j)} \beta_{ij} \sum_{(x_i, x_j)} b_{ij}(x_i, x_j) \log \left( \psi_{ij}(x_i, x_j) \right) - \sum_i \gamma_i \sum_{x_i} b_i(x_i) \log \left( \psi_i(x_i) \right) \quad (48)$$

<footnote>$^4$Note that the correct expression of the average energy corresponds to $(\beta, \gamma) = (1, 1)$; see Equation (47). However, the approximation introduced in the computation of the average energy $U_b$ could compensate for some of the wrongness in the approximation of the entropy $S_b$, eventually leading to a better approximation of the Gibbs free energy $G_b = U_b - S_b$.</footnote>
Put together, these approximations of the variational entropy $S_b$ and average energy $U_b$ lead to the following (parametric) approximation of the Gibbs free energy $G_b = U_b - S_b$ which we will refer to as the **reweighted Bethe free energy**. $G_b \approx G_b^{approx}$, with

$$G_b^{approx} = \sum_{(i,j)} \frac{1}{\alpha_{ij}} \sum_{(x_i, x_j)} b_{ij}(x_i, x_j) \log \left( \frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)} \right) - \sum_{(i,j)} \beta_{ij} \sum_{(x_i, x_j)} b_{ij}(x_i, x_j) \log \left( \psi_{ij}(x_i, x_j) \right) \\ + \sum_i \frac{1}{\kappa_i} \sum_{x_i} b_i(x_i) \log \left( b_i(x_i) \right) - \sum_i \gamma_i \sum_{x_i} b_i(x_i) \log \left( \psi_i(x_i) \right) \quad (49)$$

where $\alpha := 1/\hat{\beta}$ and $\kappa := 1/\hat{\gamma}$. Parameters $(1/\alpha; \beta; 1/\kappa; \gamma)$ are called counting numbers (see Yedidia et al. [2005]), with entropic counting numbers $(1/\alpha; 1/\kappa)$ and average energy counting numbers $(\beta; \gamma)$. For $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$, one recovers the Bethe free energy used in BP [Yedidia et al., 2003].

### E.1.2 Message-passing algorithms to minimize Gibbs free energy approximations

**Theoretical idea behind Belief Propagation** The Belief Propagation algorithm consists of trying to minimize the Bethe free energy $G_b^{Bethe}$, coming from the Bethe approximation of the Gibbs free energy, with a fixed point equation. The minimization of the Bethe free energy results into a probability distribution $b(\mathbf{x})$. The underlying idea is the following: because the variational distribution $b(\mathbf{x})$ and the initial distribution $p(\mathbf{x})$ are ‘close’ according to the KL-divergence, they are thought to have ‘close’ marginals.

In more detail, as shown below in the paragraph “Demonstration: from Gibbs approximation to Reweighted BP”, we minimize the Bethe free energy $G_b^{Bethe}$ under the constraints $\sum_{x_i} b_i(x_i) = 1$ and $\sum_{x_i} b_{ij}(x_i, x_j) = b_j(x_j)$ (necessary conditions for $b(\mathbf{x})$ to be a probability distribution). This constrained minimization leads to an equation describing the stationary points of the Bethe free energy, involving quantities called “messages” ($m_{i \to j}$) which are related to the Lagrange multipliers of the constrained optimization problem. This (fixed-point) equation is used directly to define the Belief Propagation algorithm; see Equations (2) and (3), respectively message update equation and expression of the belief (approximate marginal probability).

BP is closely linked to the Bethe free energy $G_b^{Bethe}$: fixed points of Loopy BP are stationary points of the Bethe free energy [Yedidia et al., 2001], and stable fixed points of Loopy BP correspond to minima of the Bethe free energy [Heskes, 2002].

This leads to several theoretical considerations. First, the fixed points of BP correspond to stationary points of the Bethe free energy $G_b^{Bethe}$ [Yedidia et al., 2003]: the message update equation in BP is not a gradient descent procedure on the Bethe free energy, but instead a fixed-point equation. In fact, it was later shown that stable fixed points of BP are minima of the Bethe free energy [Heskes et al., 2002]. In the case where the system has one fixed point and this fixed point is stable, the BP procedure is guaranteed to converge to the unique fixed point, that is, to the global minimum of the Bethe free energy. However, in the general case, BP might converge only to a local optimum of the Bethe free energy, or not converge at all.

Note that in the case where the initial distribution can be represented by a tree, Equation (42) is not an approximation. This explains why BP is exact when applied to acyclic probabilistic graphs.

**Idea behind Reweighted Belief Propagation** As seen above, the Belief Propagation algorithm comes from approximating the entropy term $S_b$ in the Gibbs free energy, by making the hypothesis that the variational distribution $b(\mathbf{x})$ can be written as if the associated probabilistic graph were a tree, i.e., were cycle-free. Because of that, Belief Propagation performs exactly when applied to acyclic graphs, but performs poorly in many other cases (where the approximation becomes bad).

The idea of the Reweighted Belief Propagation algorithms is to provide a better approximation to the Gibbs free energy, by reweighting terms in the expression of the variational entropy, and possibly, as considered here, by also reweighting terms in the expression of the variational average energy (which is initially exact but its approximation can help compensate for some of the remaining approximation left by approximating only the variational entropy).

The message update equations of Reweighted BP are simply fixed-point equations of $G_b^{approx}$, meaning that the beliefs computed by Reweighted BP are stationary points of the reweighted Bethe energy
$G_b^{approx}$. This demonstration is similar to the one for BP [Yedidia et al., 2001, 2003], with the additional parameters $(\alpha, \kappa, \beta, \gamma)$. Just as BP is closely linked to the Bethe free energy $G_b^{Bethe}$, Reweighted BP is closely linked to its generalization $G_b^{approx}$. The idea (or rather, hope) of Reweighted BP is that a closer approximation to the Gibbs free energy will lead to closer marginal probabilities.

**Demonstration: from Gibbs approximation to Reweighted BP** In what follows, we derive the message-passing update equations of the Reweighted BP algorithm (which includes as special cases Fractional BP, Power EP and α-BP, among others) based on the reweighted Bethe approximation of the Gibbs free energy $G_b^{approx}$; see Equation (49). This of course includes the derivation of the BP algorithm (for which the counting numbers are equal to 1). We remind the reader that we consider the case where the probability distribution is a Markov Random Field with at most pairwise interactions (note that variables $x_1, x_2, \dots, x_n$ are not necessarily binary), and that Appendix G contains a similar demonstration for general factor graphs, in which interactions can be of higher order.

The goal is to minimize the Gibbs free energy approximation $G_b^{approx}$ of Equation (49) under some constraints (normalization constraint $\sum_{x_i} b_i(x_i) = 1$ and marginalization constraint $\sum_{x_i} b_{ij}(x_i, x_j) = b_j(x_j)$). Therefore, we form a Lagrangian to take these constraints into account by adding Lagrange multipliers $(\mu, \lambda)$ to $G_b^{approx}$, where $\{\mu_i\}$ correspond to the normalization constraints and $\{\lambda_{i \to j}(x_j)\}$ correspond to the marginalization constraints. The Lagrangian writes:

$$ \mathcal{L} = G_b^{approx} + \sum_i \mu_i (\sum_{x_i} b_i(x_i) - 1) + \sum_{(i,j)} \sum_{x_j} \lambda_{i \to j}(x_j) (\sum_{x_i} b_{ij}(x_i, x_j) - b_j(x_j)) + \sum_{(i,j)} \sum_{x_i} \lambda_{j \to i}(x_i) (\sum_{x_j} b_{ij}(x_i, x_j) - b_i(x_i)) \quad (50) $$

We compute the partial derivatives of the Lagrangian:

$$ \left\{ \begin{aligned} \frac{\partial \mathcal{L}}{\partial b_i(x_i)} &= - \sum_{j \in \mathcal{N}(i)} \frac{1}{\alpha_{ij}} + \frac{1}{\kappa_i} + \frac{1}{\kappa_i} \log(b_i(x_i)) - \gamma_i \log(\psi_i(x_i)) + \mu_i - \sum_{j \in \mathcal{N}(i)} \lambda_{j \to i}(x_i) \\ \frac{\partial \mathcal{L}}{\partial b_{ij}(x_i, x_j)} &= \frac{1}{\alpha_{ij}} + \frac{1}{\alpha_{ij}} \log\left(\frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)}\right) - \beta_{ij} \log(\psi_{ij}(x_i, x_j)) + \lambda_{j \to i}(x_i) + \lambda_{i \to j}(x_j) \end{aligned} \right. $$

and cancel these partial derivatives to obtain the following expression for the unitary beliefs:

$$ b_i(x_i) \propto \psi_i(x_i)^{\kappa_i \gamma_i} \prod_{k \in \mathcal{N}(i)} \exp(\kappa_i \lambda_{k \to i}(x_i)) \quad (51) $$

and for the pairwise beliefs:

$$ b_{ij}(x_i, x_j) \propto \psi_{ij}(x_i, x_j)^{\alpha_{ij} \beta_{ij}} \psi_i(x_i)^{\kappa_i \gamma_i} \psi_j(x_j)^{\kappa_j \gamma_j} \prod_{k \in \mathcal{N}(i) \setminus j} \exp(\kappa_i \lambda_{k \to i}(x_i)) \prod_{k \in \mathcal{N}(j) \setminus i} \exp(\kappa_j \lambda_{k \to j}(x_j)) \\ \times \exp(\lambda_{j \to i}(x_i)(\kappa_i - \alpha_{ij})) \exp(\lambda_{i \to j}(x_j)(\kappa_j - \alpha_{ij})) \quad (52) $$

Now defining the messages as a function of the Lagrange multipliers $m_{j \to i}(x_i) \equiv \exp(\lambda_{j \to i}(x_i))$, the approximate marginals $b_i(x_i)$ and approximate pairwise marginals $b_{ij}(x_i, x_j)$ can be written simply as:

$$ b_i(x_i) \propto \left( \psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i)} m_{k \to i}(x_i) \right)^{\kappa_i} \quad (53) $$

$$ b_{ij}(x_i, x_j) \propto \psi_{ij}(x_i, x_j)^{\alpha_{ij} \beta_{ij}} \psi_i(x_i)^{\kappa_i \gamma_i} \psi_j(x_j)^{\kappa_j \gamma_j} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i)^{\kappa_i} \\ \times \prod_{k \in \mathcal{N}(j) \setminus i} m_{k \to j}(x_j)^{\kappa_j} m_{j \to i}(x_i)^{\kappa_i - \alpha_{ij}} m_{i \to j}(x_j)^{\kappa_i - \alpha_{ij}} \quad (54) $$
Eventually, thanks to the constraint $\sum_{x_i} b_{ij}(x_i, x_j) = b_j(x_j)$, we obtain the relation between the messages $m$:

$$m_{i \to j}(x_j)^{\kappa_j} \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}\beta_{ij}} (\psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}/\kappa_i})^{\kappa_i} \right) (m_{i \to j}(x_j))^{\kappa_j - \alpha_{ij}} \quad (55)$$

$$\Leftrightarrow m_{i \to j}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}\beta_{ij}} (\psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}/\kappa_i})^{\kappa_i} \right)^{1/\alpha_{ij}} \quad (56)$$

The Reweighted BP algorithm consists of running iteratively the fixed-point equation (56):

$$m_{i \to j}^{\text{new}}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}\beta_{ij}} (\psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}^{\text{old}}(x_i) m_{j \to i}^{\text{old}}(x_i)^{1-\alpha_{ij}/\kappa_i})^{\kappa_i} \right)^{1/\alpha_{ij}} \quad (57)$$

Note that one could also use directly Equation (55) instead of (56) to define the Reweighted BP algorithm:

$$m_{i \to j}^{\text{new}}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}\beta_{ij}} (\psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}^{\text{old}}(x_i) m_{j \to i}^{\text{old}}(x_i)^{1-\alpha_{ij}/\kappa_i})^{\kappa_i} \right)^{1/\kappa_j} (m_{i \to j}^{\text{old}}(x_j))^{1-\alpha_{ij}/\kappa_j} \quad (58)$$

In fact, Equations (58) and (57) correspond respectively to the damped (with a particular damping value) versus undamped update equation; see section F. There is no absolute better choice: fixed points do not depend on the amount of damping, and damping might provide better convergence properties but might slow down the system in cases where the algorithm converges without damping (see section F). In the simulations, we take no damping, except stated otherwise, to distinguish between the effect of damping and the effect of the algorithm (e.g., Circular versus Fractional BP).

## E.2 Reweighted BP algorithms

**Definition** As shown above, minimizing the reweighted Bethe free energy $G_b^{\text{approx}}$ under the constraints $\sum_{x_i} b_i(x_i) = 1$ and $\sum_{x_i} b_{ij}(x_i, x_j) = b_j(x_j)$ leads to a general **Reweighted BP** algorithm defined by iterating the following (fixed-point) update equation:

$$m_{i \to j}^{\text{new}}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}\beta_{ij}} (\psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}/\kappa_i})^{\kappa_i} \right)^{1/\alpha_{ij}} \quad (59)$$

Approximate marginals or beliefs are computed using:

$$b_i(x_i) \propto \left( \psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i)} m_{k \to i}(x_i) \right)^{\kappa_i} \quad (60)$$

**Special cases of Reweighted BP** Special cases of the general reweighted Bethe approximation lead to different algorithms. We mention here several *Reweighted BP* algorithms: Fractional BP, Power-EP, $\alpha$-BP, Tree-Reweighted BP, and Variational message-passing (mean-field method); see below for more details (mathematical definitions and relationships between algorithms). In these algorithms, parameters $\beta$ and $\gamma$ take value 1, that is, can be related to Gibbs free energy approximations which only approximate the variational entropy part, not the average energy. **Fractional BP** [Wiegerinck and Heskes, 2002] comes as a result of approximating the entropy term of the Gibbs free energy by hypothesizing that the variational distribution $b(x)$ writes according to Equation (9) in the particular case $\kappa = 1$, that is, with the single parameter $\alpha$. **Reweighted BP** [Loh and Wibisono, 2014] is a message-passing algorithm based on approximating the variational entropy as a weighted sum of local entropies. This is equivalent to the Gibbs approximation above with parameters $\alpha, \kappa$ but as in previous algorithms with $(\beta, \gamma) = (1, 1)$. **Tree-reweighted BP** [Wainwright. et al., 2002, 2003, Wainwright et al., 2005] aims at computing an upper bound of the log-partition function ($\log(Z)$), by
using the same algorithm; the only difference is that $\alpha_{ij}$ symbolizes the inverse appearance probability of edge $(i, j)$ in the set of spanning trees, which imposes $\alpha_{ij} \ge 1$ [Minka, 2005]. **Convex BP** [Hazan and Shashua, 2008, 2010], which admits as special case Tree-reweighted BP, is a general message-passing algorithm based on convex free energies. **Power Expectation Propagation** [Minka and Lafferty, 2002, Minka, 2004] and **α-BP** [Liu et al., 2019, 2020] relate to the problem of minimizing the α-divergence between $p(\mathbf{x})$ and $b(\mathbf{x})$ instead of the KL-divergence of BP. α-BP and Power EP differ in that α-BP considers the projection into factorized distributions, while Power EP considered projections into exponential families. Despite the fact that these two last algorithms do not come directly from the reweighting of the Bethe free energy, their message update equation is the same as previously mentioned algorithms; see also [Minka, 2005]. In fact, Power EP, α BP, and Fractional BP are identical up to the amount of damping, where damping consists of taking partial message update steps [Murphy et al., 1999]; see below. They all correspond to $(\kappa, \beta, \gamma) = (1, 1, 1)$. **Variational message-passing** [Winn and Bishop, 2005], which is the message-passing version of the mean-field method [Peterson and Anderson, 1987], is the same as Fractional BP for $\alpha \to 0$ [Wiegerinck and Heskes, 2002, Minka, 2005]. Lastly, BP corresponds to $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$.

**Special case of Reweighted BP: Binary case, pairwise factors** With the additional hypothesis that variables $x_i$ are binary, Equations (59) and (60) can be written simply in the log domain:

$$M_{i \to j}^{\text{new}} = g(B_i - \alpha_{ij} M_{j \to i}, (\psi_{ij})^{\beta_{ij}}, \alpha_{ij}) \quad (61)$$

$$B_i = \kappa_i \left( \sum_{j \in \mathcal{N}(i)} M_{j \to i} + \gamma_i M_{\text{ext} \to i} \right) \quad (62)$$

where quantities $M$ and $B$ are given by $B_i \equiv \frac{1}{2} \log \left( \frac{b_i(x_i+1)}{b_i(x_i-1)} \right)$, $M_{i \to j} \equiv \frac{1}{2} \log \left( \frac{m_{i \to j}(x_j+1)}{m_{i \to j}(x_j-1)} \right)$, and $M_{\text{ext} \to i} \equiv \frac{1}{2} \log \left( \frac{\psi_i(x_i+1)}{\psi_i(x_i-1)} \right)$. Function $g$ is a sigmoidal function of its first argument, with parameters $\psi_{ij}$ (interaction between variables $x_i$ and $x_j$), $\alpha_{ij}$ and $\beta_{ij}$ (counting numbers associated to the edge $(i, j)):

$$g(x, \psi_{ij}^{\beta_{ij}}, \alpha_{ij}) = \frac{1}{2\alpha} \log \left( \frac{\left(\psi_{ij}(x_i = +1, x_j = +1)\right)^{\alpha_{ij}\beta_{ij}} e^{2x} + \left(\psi_{ij}(x_i = -1, x_j = +1)\right)^{\alpha_{ij}\beta_{ij}}}{\left(\psi_{ij}(x_i = +1, x_j = -1)\right)^{\alpha_{ij}\beta_{ij}} e^{2x} + \left(\psi_{ij}(x_i = -1, x_j = -1)\right)^{\alpha_{ij}\beta_{ij}}} \right) \quad (63)$$

In the even more particular case of the Ising model, for which $\psi_{ij}(x_i, x_j) \propto \exp(J_{ij}x_i x_j)$, then

$$g(x, \psi^{\beta}, \alpha) = \frac{1}{\alpha} \operatorname{arctanh} \left( \tanh(\alpha \beta J) \tanh(x) \right) \quad (64)$$

*Proof.* Use Equation (63), $\log(x) = 2 \operatorname{arctanh} \left(\frac{x-1}{x+1}\right)$, and $\tanh(x) = 2\sigma(2x) - 1 = \frac{e^{2x}-1}{e^{2x}+1}$. $\square$

## E.3 Mathematical definitions of several Reweighted BP algorithms

## E.4 Definitions

**BP** The message update equation associated to the Belief Propagation algorithm is:

$$m_{i \to j}^{\text{new}}(x_j) \propto \sum_{x_i} \psi_{ij}(x_i, x_j) \psi_i(x_i) \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) \quad (65)$$

**Power-EP** Power Expectation Propagation [Minka, 2004] is an extension of Expectation Propagation [Minka, 2001] which consists of minimizing the α-divergence between the true probability distribution $p(\mathbf{x})$ and the variational distribution $b(\mathbf{x})$, instead of the KL-divergence as in Expectation Propagation and Belief Propagation. The associated message update equation is:

$$m_{i \to j}^{\text{new}}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}} \psi_i(x_i) \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}} \right)^{1/\alpha_{ij}} \quad (66)$$
**Fractional BP and alpha-BP** The Fractional BP algorithm [Wiegerinck and Heskes, 2002] is defined as a result of approximating the entropy of the variational distribution b(x). As explained above, α-BP is derived differently, but has the same update equation as Fractional BP.

The message update equations of Fractional BP and α-BP are respectively defined in Equation (17) of Wiegerinck and Heskes [2002]) and in Equation (17) of Liu et al. [2020], and take the following form:

$$m_{i \to j}^{\text{new}}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}} \psi_i(x_i) \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}} \right) m_{i \to j}(x_j)^{1-\alpha_{ij}} \quad (67)$$

Note that Fractional BP is derived as this general Reweighted BP algorithm presented above (but with $(\kappa, \gamma, \beta) = (1, 1, 1)$), and is defined with damping, i.e., uses Equation (58) rather than Equation (57).

**Reweighted BP** The message update equation associated to the general Reweighted BP algorithm is recalled here:

$$m_{i \to j}^{\text{new}}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}\beta_{ij}} (\psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}/\kappa_i})^{\kappa_i} \right)^{1/\alpha_{ij}} \quad (68)$$

**Variational message-passing** Variational message-passing [Winn and Bishop, 2005], or mean-field inference, is the message-passing version of the mean-field method [Peterson and Anderson, 1987]. Its message update equation is:

$$m_{i \to j}^{\text{new}}(x_j) \propto \exp \left( \sum_{x_i} \log (\psi_{ij}(x_i, x_j)) \psi_i(x_i) \prod_{k \in \mathcal{N}(i)} m_{k \to i}(x_i) \right) \quad (69)$$

Mean-field is known to be overconfident and perform poorly compared to BP [Weiss, 2001, Mooij and Kappen, 2004].

## E.5 Relations between Reweighted BP algorithms

**Equivalence between damped Fractional BP, Power EP, and α-BP** The Fractional Belief Propagation algorithm closely relates to several approximate inference algorithms that themselves extended BP, such as Power EP [Minka and Lafferty, 2002, Minka, 2004], α-BP [Liu et al., 2019, 2020], Tree-reweighted BP [Wainwright. et al., 2002, 2003, Wainwright et al., 2005], and Variational message-passing [Winn and Bishop, 2005]. Here we show the equivalence between Fractional BP, Power EP, and α-BP, if we consider a damped version of the algorithms. These models are conceptually very similar to each other [Minka, 2005]. In fact, the three algorithms are identical up to the amount of damping, where damping consists of taking partial message update steps [Murphy et al., 1999]. Indeed, damped Fractional BP, defined similarly to damped BP (see section F) from the undamped Fractional BP message update equation in Equation (66), is written:

$$m_{i \to j}^{\text{new}}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\alpha_{ij}} \psi_i(x_i) \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}} \right)^{(1-\epsilon_{i \to j})/\alpha_{ij}} \times m_{i \to j}(x_j)^{\epsilon_{i \to j}} \quad (70)$$

and corresponds to Equation (67) (defining Fractional BP in Wiegerinck and Heskes [2002]) for the particular damping value $\epsilon_{i \to j} = 1 - \alpha_{ij}$. When it comes to Power EP [Minka, 2004], the damping does not appear from the derivations of the algorithm but Power EP is still presented with the possibility of having (any) damping $\epsilon$ (see Equations (22) and (23) of Minka [2004]). It is stated in the paper that with the particular value of damping $\epsilon_{i \to j} = 1 - \alpha_{ij}$, the algorithm is “convenient computationally and tends to have good convergence properties”. Interestingly, with this particular value of damping, it is shown that the Power EP (and equivalently, α-BP and Fractional BP) proposed implements a minimization of the α-divergence (see also the work on α-BP [Liu et al., 2020] as well as Minka [2005]).
Again, the damping values only alters the convergence properties of the algorithm but not the fixed points. A consequence is that if the system has a single fixed point and the undamped algorithm converges (respectively damped), then the damped (respectively undamped) algorithm converges and the convergence value is identical between the damped and undamped algorithms. In the rest of this work, we will consider these algorithms as one: all algorithms implemented in the simulations do not have damping in order not to confound the effect of having different algorithms (for instance, Circular BP versus Fractional BP) with the effects coming from using different damping values.

**Tree-reweighted BP as special case of Fractional BP** Tree-reweighted BP [Wainwright. et al., 2002, 2003, Wainwright et al., 2005] is a particular case of the algorithms mentioned above. The message update equation of Tree-reweighted BP is identical to Equation (66). The only difference is that in Tree-reweighted BP, $\alpha_{ij}$ symbolizes the inverse appearance probability of edge $(i, j)$ in the set of spanning trees and therefore it imposes the constraint $\alpha_{ij} \ge 1$ [Minka, 2005]. Note that Fractional BP does not impose any constraints on the value of $\alpha_{ij}$, which even could be negative.

**BP and Reweighted BP** The special case of Belief Propagation is recovered for $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$.

**Fractional BP and Reweighted BP** The Reweighted BP algorithm generalizes Fractional BP [Wiegerinck and Heskes, 2002], Power EP [Minka, 2004] and $\alpha$-BP [Liu et al., 2019] which all correspond to $(\kappa, \beta, \gamma) = (1, 1, 1)$, and more particularly use the damped message update equation (55) rather than its undamped version (56) (see damping in section F).

Note that the update equation for Power EP (Equation (66)) is not identical to the one given in the Fractional BP paper [Wiegerinck and Heskes, 2002], but is the same up the amount of damping in the algorithm.

**Variational message-passing as special case of Fractional BP** Variational message-passing [Winn and Bishop, 2005] is the same as Fractional BP for $\alpha = 0$ [Wiegerinck and Heskes, 2002, Minka, 2005]. More precisely, Fractional BP is not defined for $\alpha = 0$. However, taking $\alpha \to 0$ and $(\kappa, \beta, \gamma) = (1, 1, 1)$ in the Gibbs free energy approximation (Equation (49)) leads to the Variational mean-field free energy $G_b^{MF}$, where

$$G_b^{MF} = - \sum_{(i,j)} \sum_{(x_i, x_j)} b_i(x_i) b_j(x_j) \log (\psi_{ij}(x_i, x_j)) + \sum_i \sum_{x_i} b_i(x_i) \log (b_i(x_i)) - \sum_i \sum_{x_i} b_i(x_i) \log (\psi_i(x_i)) \quad (71)$$

Indeed, if $\alpha \to 0$, then terms in $\frac{b_{ij}(x_i, x_j)}{b_i(x_i)b_j(x_j)}$ dominate in the Gibbs free energy approximation, which forces the variational marginals to factorize: $b_{ij}(x_i, x_j) = b_i(x_i)b_j(x_j)$. Eventually, similarly to section E.1.2, one obtains the Variational message-passing update equation given in section E.4. Note that variational message-passing is different from Circular BP with $\alpha = 0$.

## E.6 Circular BP as approximation of Reweighted BP

**Definition of Circular BP** Circular BP consists of running the following message update equation:

$$m_{i \to j}^{\text{new}}(x_j) \propto \sum_{x_i} \psi_{ij}(x_i, x_j)^{\beta_{ij}} \left( \psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}/\kappa_i} \right)^{\kappa_i} \quad (72)$$

The approximate marginals or beliefs are computed using:

$$b_i(x_i) \propto \left( \psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i)} m_{k \to i}(x_i) \right)^{\kappa_i} \quad (73)$$

This algorithm admits as particular case BP (which corresponds to $(\alpha, \beta, \kappa, \gamma) = (1, 1, 1, 1)$) as well as the algorithm from Jardri and Denève [2013], also called Circular BP (which corresponds to $(\beta, \kappa, \gamma) = (1, 1, 1)$).
For probability distributions over binary variables, this translates into the following equation in the log-domain:

$$M_{i \to j}^{\text{new}} = f(B_i - \alpha_{ij} M_{j \to i}, \psi_{ij}^{\beta_{ij}}) \quad (74)$$

$$B_i = \kappa_i \left( \sum_{j \in \mathcal{N}(i)} M_{j \to i} + \gamma_i M_{\text{ext} \to i} \right) \quad (75)$$

where $B_i \equiv \frac{1}{2} \log (\frac{b_i(x_i=+1)}{b_i(x_i=-1)})$, $M_{i \to j} \equiv \frac{1}{2} \log (\frac{m_{i \to j}(x_j=+1)}{m_{i \to j}(x_j=-1)})$, and $M_{\text{ext} \to i} \equiv \frac{1}{2} \log (\frac{\psi_i(x_i=+1)}{\psi_i(x_i=-1)})$. $f$ is a sigmoidal function of its first argument given by:

$$f(x, \psi_{ij}^{\beta_{ij}}) = \frac{1}{2} \log \left( \frac{(\psi_{ij}(x_i = +1, x_j = +1))^{\beta_{ij}} e^{2x} + (\psi_{ij}(x_i = -1, x_j = +1))^{\beta_{ij}}}{(\psi_{ij}(x_i = +1, x_j = -1))^{\beta_{ij}} e^{2x} + (\psi_{ij}(x_i = -1, x_j = -1))^{\beta_{ij}}} \right) \quad (76)$$

in the general case, and

$$f(x, \psi_{ij}^{\beta_{ij}}) = \operatorname{arctanh} \left( \tanh(\beta_{ij} J_{ij}) \tanh(x) \right) \quad (77)$$

in the specific case of Ising models for which $\psi_{ij}(x_i, x_j) \propto \exp(J_{ij}x_i x_j)$.

**Proof.** Use Equation (76), $\log(x) = 2 \operatorname{arctanh} \left(\frac{x-1}{x+1}\right)$, and $\tanh(x) = 2\sigma(2x) - 1 = \frac{e^{2x}-1}{e^{2x}+1}$. $\square$

**Link to Jardri and Denève [2013]** Jardri and Denève [2013] defines an algorithm which corresponds to a special case of Circular BP: $\kappa, \beta, \gamma = (1, 1, 1)$. However, this algorithm was not designed for approximate inference, but instead to disturb the exactness of BP in graphs without cycles.

**Circular BP as approximation of Reweighted BP** Circular BP is not associated to any choice of $(\alpha, \kappa, \beta, \gamma)$ in Reweighted BP, thus does not come from any modification of the Bethe free energy (and does not relate either to any other approximation of the Gibbs free energy). Instead, as explained below, CBP can be seen as an approximation of Reweighted BP. Circular BP can be seen as an approximation to the general Reweighted BP algorithm defined above. This is not obvious from the general message update equations (see Equation (59) for Reweighted BP and Equation (72) for Circular BP), but can be justified mathematically in the case where the distribution $p(x)$ considered is an Ising model, that is, a probability distribution over binary variables, involving at most pairwise interactions, and with particular pairwise interactions $(\psi_{ij}(x_i, x_j) \propto \exp(J_{ij}x_i x_j))$. In the case of Ising models, Reweighted BP and Circular BP are indeed written identically, except Reweighted BP uses function $g$ and Circular BP uses $f$ (see Equations (61)-(62) versus (74)-(75)). However, functions $f$ and $g$ are approximately equal:

$$g(x, \beta J, \alpha) \equiv \frac{1}{\alpha} \operatorname{arctanh} \left( \tanh \left( \alpha \beta J \right) \tanh(x) \right) \approx \operatorname{arctanh} \left( \tanh \left( \beta J \right) \tanh(x) \right) \equiv f(x, \beta J) \quad (78)$$

In summary, the Circular BP algorithm strongly relates to the Reweighted BP algorithm; neither of these two algorithms generalizes the other one, but their message update equations are approximately the same.

**Circular BP and Reweighted BP** As seen in section E.6, the update equation of Circular BP algorithm approximates the one of Reweighted BP: indeed, in the Ising model case, update equations are identical up to the function used ($g$ for Reweighted BP and $f$ for Circular BP) and functions are approximately equal ($f(x, J) \approx g(x, J, \alpha)$).

**Circular BP from Jardri and Denève [2013] and Fractional BP** Likewise, Circular BP from Jardri and Denève [2013] (which corresponds to the case $(\kappa, \beta, \gamma) = (1, 1, 1)$ of the Circular BP algorithm defined in this paper) and Fractional BP (which corresponds to the case $(\kappa, \beta, \gamma) = (1, 1, 1)$ of the Reweighted BP algorithm of Equation (68) up to the damping value).
<figure><img src="image_18.png" /><figcaption>Figure S10: Relations between algorithms. Algorithms on the first line (Reweighted and Circular BP) have parameters (α, κ, β, γ). Algorithms on the second line have parameter α. Algorithms on the third line do not have parameters.</figcaption></figure>

Circular BP does not correspond to any choice of (α, κ, β, or γ). However, because functions g and f (both in the case β = 1) are similar:

$$g(x, J, \alpha) \equiv \frac{1}{\alpha} \operatorname{arctanh} \left( \tanh(\alpha\beta J) \tanh(x) \right) \approx \tanh^{-1} \left( \tanh(\beta J) \tanh(x) \right) \equiv f(x, J) \qquad (79)$$

, then the message update equation of Circular BP approximates the update equation corresponding to (κ, β, γ) = (1, 1, 1).

**Comparing Fractional BP and Circular BP** In this paragraph, we show the similarity between Fractional BP (Reweighted BP in the particular case (κ, β, γ) = (1, 1, 1)) and Circular BP (also in the particular case (κ, β, γ) = (1, 1, 1)). Note that this particular case of Circular BP corresponds to the algorithm from Jardri and Denève [2013].

The message-update equations for the two algorithms (respectively Equation (66) and Equation (4)) are very similar. It becomes even more obvious when considering an Ising model: both write as:

$$M_{i \to j}^{\text{new}} = g(B_i - \alpha_{ij} M_{j \to i}, J_{ij}, \alpha_{ij}) \qquad (80)$$

$$B_i = \sum_{j \in \mathcal{N}(i)} M_{j \to i} + M_{\text{ext} \to i} \qquad (81)$$

where $g(x, J, \alpha) \equiv \frac{1}{\alpha}f(x, \alpha J)$ for Fractional BP, and $g(x, J, \alpha = 1) = f(x, J)$ for Circular BP. Interestingly, functions g and f are close to each other:

$$g(x, J, \alpha) \equiv \frac{1}{\alpha} \operatorname{arctanh} \left( \tanh(\alpha J) \tanh(x) \right) \approx \tanh^{-1} \left( \tanh(J) \tanh(x) \right) \equiv f(x, J) \qquad (82)$$

which can be easily justified mathematically for α or J small enough, and can be seen in practice for various J and α in Figure S11A. This means that the message update equation of Fractional BP approximates the update equation of Circular BP. A consequence of the similarity between their update equations is that Circular BP and Fractional BP yield similar results; see Figure S11C and D.
A

<figure><img src="image_19.png" /></figure>

B

<figure><img src="image_20.png" /></figure>

C

<figure><img src="image_21.png" /></figure>

D

<figure><img src="image_22.png" /><figcaption>Figure S11: Comparison between Fractional BP and Circular BP. (A) Visual comparison between the update functions for Fractional BP (function *g*, which depends on the strength of the relationship *J* and on *α*) and Circular BP (function *f* which depends only on *J* and corresponds to *g* for *α* = 1) for various *α* and given *J* = 1. All functions have the identical sigmoidal shape. (B) The update functions for Circular BP (*f*) and for Fractional BP (*g*) are approximately identical, for various *J* and *α*. (C) Example of marginals produced by the Fractional BP algorithm versus produced by the Circular BP algorithm for a given random graph with randomly generated weights and random *α*. (D) Influence of the interaction strength *σ*(*J*) on the approximation: lower weights lead to a better fit between Circular BP and Fractional BP. One point represents different randomly generated *α* for the same graph.</figcaption></figure>
# F Damped Circular BP

Section F gives a definition of damped CBP.

**Notion of damping** Let a discrete system:

$$x_{t+1} = f(x_t)$$

A common technique to improve convergence of the system is to take partial (or damped) update steps:

$$x_{t+1} = (1 - \epsilon)f(x_t) + \epsilon x_t$$

with $\epsilon \in [0; 1[$. This procedure, called damping, does not modify the fixed points of the system: a fixed point of the original system is a fixed point of the damped system, and reciprocally.

A parallel can be drawn between the damped discrete system and the following continuous system:

$$\tau \dot{x}(t) = -x(t) + f(x(t)) \quad (83)$$

Indeed, this continuous system can be discretized (Euler approximation) into:

$$x_{t+\delta t} = (1 - \epsilon)f(x_t) + \epsilon x_t \quad (84)$$

with $\epsilon = 1 - \delta t/\tau$ (i.e., $\tau = \delta t/(1 - \epsilon)$). This corresponds to the discrete damped system.

**Damped Circular BP** We start by rewriting the message update equation for Circular BP (Equation (4)) without damping:

$$m_{i \to j}^{\text{new}}(x_j) \propto \sum_{x_i} \psi_{ij}(x_i, x_j)^{\beta_{ij}} \left( \psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}/\kappa_i} \right)^{\kappa_i} \quad (85)$$

and its equivalent in the log-domain for probability distributions with binary variables (Equation (6)):

$$M_{i \to j}^{\text{new}} = f(B_i - \alpha_{ij} M_{j \to i}, \beta_{ij} J_{ij}) \quad (86)$$

where $M_{i \to j} \equiv \frac{1}{2} \log(\frac{m_{i \to j}(x_j=+1)}{m_{i \to j}(x_j=-1)})$, $M_{\text{ext} \to i} \equiv \frac{1}{2} \log(\frac{\psi_i(x_i=+1)}{\psi_i(x_i=-1)})$, $B_i \equiv \frac{1}{2} \log(\frac{b_i(x_i=+1)}{b_i(x_i=-1)}) = \kappa_i (\sum_{j \in \mathcal{N}(i)} M_{j \to i} + \gamma_i M_{\text{ext} \to i})$, and function f is given in Equation (76).

Damping is a technique commonly used while running Belief Propagation (special case $(\alpha, \kappa, \gamma, \beta) = (1, 1, 1, 1)$), and consists of taking partial message update steps in the log-space: $M_{i \to j}^{\text{new}} = f(B_i - M_{j \to i}, \beta_{ij} J_{ij})$ becomes $M_{i \to j}^{\text{new}} = (1 - \epsilon)f(B_i - M_{j \to i}, \beta_{ij} J_{ij}) + \epsilon M_{i \to j}$. Damping improves the convergence properties of BP. We talk about damped Belief Propagation in this case [Murphy et al., 1999].

Similarly to damped BP, a damped Circular BP algorithm can be defined, by taking partial message update steps in the log-space. The message update equation for the damped algorithm is:

$$m_{i \to j}^{\text{new}}(x_j) \propto \left( \sum_{x_i} \psi_{ij}(x_i, x_j)^{\beta_{ij}} \left( \psi_i(x_i)^{\gamma_i} \prod_{k \in \mathcal{N}(i) \setminus j} m_{k \to i}(x_i) m_{j \to i}(x_i)^{1-\alpha_{ij}/\kappa_i} \right)^{\kappa_i} \right)^{1-\epsilon_{i \to j}} \times m_{i \to j}(x_j)^{\epsilon_{i \to j}} \quad (87)$$

where $\epsilon_{i \to j}$ is the damping factor associated to the oriented edge $i \to j$ ($\epsilon = 0$ means no damping, i.e., standard CBP) and is often taken uniformly over the edges. It is equivalent in the log-domain to the following equation:

$$M_{i \to j}^{\text{new}} = (1 - \epsilon_{i \to j}) f(B_i - \alpha_{ij} M_{j \to i}, \beta_{ij} J_{ij}) + \epsilon_{i \to j} M_{i \to j} \quad (88)$$

where

$$B_i = \kappa_i \left( \sum_{j \in \mathcal{N}(i)} M_{j \to i} + \gamma_i M_{\text{ext} \to i} \right) \quad (89)$$
# G Circular BP in general factor graphs

In the main paper, we only considered Ising models, which are a particular class of pairwise factors graphs (that is, probability distributions $p(x)$ which could be written as the product of unitary and pairwise potentials: $p(\mathbf{x}) \propto \prod_{(i,j)} \psi_{ij}(x_i, x_j) \prod_i \psi_i(x_i)$), acting on binary variables where potentials take a particular form ($\psi_{ij}(x_i, x_j) = \exp(J_{ij}x_i x_j)$). In this section, we present the general case: we formulate the Circular BP algorithm in general factor graphs (instead of pairwise graphs) $p(\mathbf{x}) \propto \prod_c \psi_c(\mathbf{x}_c)$ where variables can be continuous.

Circular BP algorithm can be defined on general factor graphs similarly to the pairwise factor case:

$$m_{\psi_c \to x_j}^{\text{new}}(x_j) \propto \sum_{\mathbf{x}_c \setminus x_j} \psi_c(\mathbf{x}_c)^{\beta_c} \left( \prod_{x_i \in \mathcal{N}(\psi_c) \setminus x_j} \psi_i(x_i)^{\gamma_i} \prod_{x_i \in \mathcal{N}(\psi_c) \setminus x_j} \left[ \prod_{\psi_d \in \mathcal{N}(x_i) \setminus \psi_c} m_{\psi_d \to x_i}(x_i) \times m_{\psi_c \to x_i}(x_i)^{1-\alpha_c/\kappa_i} \right] \right)^{\kappa_i} \quad (90)$$

and produces the following approximate marginals:

$$b_i(x_i) \propto \psi_i(x_i)^{\gamma_i \kappa_i} \prod_{\psi_c \in \mathcal{N}(x_i)} m_{\psi_c \to x_i}(x_i)^{\kappa_i} \quad (91)$$

Note that Equation (90) can be rewritten as a system of equations involving both factor-to-variable messages and variable-to-factor messages, instead of factor-to-variables messages alone:

$$m_{\psi_c \to x_i}^{\text{new}}(x_i) \propto \sum_{\mathbf{x}_c \setminus x_i} \psi_c(\mathbf{x}_c)^{\beta_c} \prod_{x_j \in \mathcal{N}(\psi_c) \setminus x_i} m_{x_j \to \psi_c}(x_j) \quad (92)$$

$$m_{x_j \to \psi_d}^{\text{new}}(x_j) \propto \psi_j(x_j)^{\gamma_j \kappa_j} \left( \prod_{\psi_c \in \mathcal{N}(x_j) \setminus \psi_d} m_{\psi_c \to x_j}(x_j) \right)^{\kappa_j} m_{\psi_d \to x_j}(x_j)^{\kappa_j - \alpha_d} \quad (93)$$

## G.1 Inference in graphical models

A general Markov random field $p(\mathbf{x})$ can be written:

$$p(\mathbf{x}) \propto \prod_c \psi_c(\mathbf{x}_c) \quad (94)$$

Potentials $\psi_c$ are called "factors" and are associated to a clique $\mathbf{x}_c$ ($\mathbf{x}_c$ is a group of variables $\{x_i\}$ or simply one variable $x_i$). As Figure S12 shows, the probability distribution can be represented graphically as a factor graph composed of variable nodes $x_i$ and factor nodes $\psi_c$, with links between $\psi_c$ and all the variable nodes present in (vector) $\mathbf{x}_c$.

## G.2 Belief Propagation

Belief Propagation can be defined on a factor graph [Kschischang et al., 2001]. At every iteration, the algorithm updates messages from factor nodes to variable nodes and from variable nodes to factor nodes:

$$m_{\psi_c \to x_i}^{\text{new}}(x_i) \propto \sum_{\mathbf{x}_c \setminus x_i} \psi_c(\mathbf{x}_c) \prod_{x_j \in \mathcal{N}(\psi_c) \setminus x_i} m_{x_j \to \psi_c}(x_j) \quad (95)$$

$$m_{x_j \to \psi_d}^{\text{new}}(x_j) \propto \psi_j(x_j) \prod_{\psi_c \in \mathcal{N}(x_j) \setminus \psi_d} m_{\psi_c \to x_j}(x_j) \quad (96)$$

where $\mathcal{N}(x)$ are the neighbors of node $x$ in the factor graph. Neighbors of variable nodes are factor nodes, and reciprocally.

Once messages have converged (or at some given maximum iteration), approximate marginal probabilities or beliefs are computed as:

$$b_i(x_i) \propto \psi_i(x_i) \prod_{\psi_c \in \mathcal{N}(x_i)} m_{\psi_c \to x_i}(x_i) \quad (97)$$
<figure><img src="image_23.png" /><figcaption>Figure S12: Belief Propagation and Circular Belief Propagation on a general factor graph. (A) Example of a factor graph with higher-order interactions (i.e., not only unitary and pairwise factors), here representing the probability distribution $p(\mathbf{x}) = \psi_{125}(x_1, x_2, x_5)\psi_{23}(x_2, x_3)\psi_{14}(x_1, x_4)\psi_{24}(x_2, x_4)\psi_1(x_1)\psi_3(x_3)\psi_4(x_4)\psi_5(x_5)$. (B) Belief Propagation updates a given variable-to-factor message (from $x$ to $\psi$, dotted black line) according to the messages received by node $x$ from other factor nodes than $\psi$ (full black lines). For Circular BP, the opposite message (coming from $\psi$ to $x$, full red line) is also taken into account. (C) Both BP and Circular BP update factor-to-variable messages (from $\psi$ to $x$, dotted black line) according to the messages collected by $\psi$ from other variable nodes than $x$ (full black lines), and the interaction function $\psi$.</figcaption></figure>

A message $m_{\psi_c \to x_i}(x_i)$ from factor node to variable node correspond to probabilistic information about variable $x_i$ collected by the factor node $\psi_c$. The message is based on the information available elsewhere in the graph (observed variables, prior distribution over variables) received by $\psi_c$, and takes into account the probabilistic interactions between $x_i$ and its neighbors (i.e., the interaction factor $\psi_c$).

A message $m_{x_j \to \psi_d}(x_j)$ from variable node to factor node is simply the sum (in the log-domain) of the local information at $x_j$ (e.g., noisy observation or prior) with the messages received by $x_j$ from all factors neighboring $x_j$ except $\psi_d$.

Note that if factors are all pairwise (case considered in the main text: $\psi_c = \psi_{ij}$) then BP equations (95) and (96) can be written with messages going from factor to variable node only. We recover Equations (2) and (3) from the main text by defining $m_{i \to j}(x_j) \equiv m_{\psi_c \to x_j}(x_j)$.

## G.3 Reweighted BP

We consider here a modification of Belief Propagation based on a parametric approximation of the entropy of the approximating distribution $b(\mathbf{x})$, with the parameters $(\kappa, \alpha, \beta, \gamma)$. $\kappa_i$ and $\gamma_i$ are assigned to the variable node $x_i$, while $\alpha_c$ and $\beta_c$ are assigned to the factor $\psi_c$. We obtain the following modified update equations (see following paragraph for the demonstration):

$$ m_{\psi_c \to x_i}^{\text{new}}(x_i) \propto \left( \sum_{\mathbf{x}_c \setminus x_i} \psi_c(\mathbf{x}_c)^{\alpha_c} \beta_c \prod_{x_j \in \mathcal{N}(\psi_c) \setminus x_i} m_{x_j \to \psi_c}(x_j) \right)^{1/\alpha_c} \quad (98) $$

$$ m_{x_j \to \psi_d}^{\text{new}}(x_j) \propto \psi_j(x_j)^{\gamma_j \kappa_j} \left( \prod_{\psi_c \in \mathcal{N}(x_j) \setminus \psi_d} m_{\psi_c \to x_j}(x_j) \right)^{\kappa_j} m_{\psi_d \to x_j}(x_j)^{\kappa_j - \alpha_d} \quad (99) $$

and beliefs (approximate marginal probabilities) are computed using:

$$ b_i(x_i) \propto \psi_i(x_i)^{\gamma_i \kappa_i} \left( \prod_{\psi_c \in \mathcal{N}(x_i)} m_{\psi_c \to x_i}(x_i) \right)^{\kappa_i} \quad (100) $$
The special case where $(\kappa, \alpha, \beta, \gamma) = (1, 1, 1, 1)$ corresponds to BP.

Similarly to above, if factors are all pairwise (case considered in the main text) then the Reweighted BP equations (98) and (99) can be written with messages going from factor to variable node only. We recover Equations (59) and (60) from the main text by defining $m_{i \to j}(x_j) \equiv m_{\psi_c \to x_j}(x_j)$.

## G.4 Theoretical background for Reweighted BP

Here we provide the theoretical foundations underlying the modification of BP given in Equations (98), (99), and (100).

Note that in the demonstration, we cover all the special cases of Reweighted BP, including BP (for which $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$), as well as Fractional BP, Power EP and $\alpha$-BP (which all correspond to $(\kappa, \beta, \gamma) = (1, 1, 1)$).

As stated in section E of the Appendix, the approach taken by BP to compute marginals $p(\mathbf{x})$, whose formula is known but whose marginals are hard to compute, is to approximate $p(\mathbf{x})$ with distribution $b(\mathbf{x})$ (called the variational distribution) whose marginals are easier to compute. The Gibbs free energy (that we would like to minimize) is given by:

$$G_b = U_b - S_b \qquad (101)$$

where the variational average energy $U_b$ can be computed easily:

$$ \begin{aligned} U_b &= \sum_{\mathbf{x}} b(\mathbf{x}) E(\mathbf{x}) \\ &= -\sum_{\mathbf{x}} b(\mathbf{x}) \sum_{\text{cliques } c} \log(\psi_c(\mathbf{x}_c)) - \sum_{\mathbf{x}} b(\mathbf{x}) \sum_{i} \log(\psi_i(x_i)) \\ &= -\sum_{\text{cliques } c} \sum_{\mathbf{x}_c} b_c(\mathbf{x}_c) \log(\psi_c(\mathbf{x}_c)) - \sum_{i} \sum_{\mathbf{x}_i} b_i(x_i) \log(\psi_i(x_i)) \end{aligned} \qquad (102) $$

, contrary to the variational entropy $S_b$.

As it is not possible to easily compute $S_b$, Belief Propagation estimates it as if the factor graph representing $b(\mathbf{x})$ was a tree (i.e., was acyclic). This means that:

$$b(\mathbf{x}) \approx \prod_{\text{cliques } c} \left( \frac{b_c(\mathbf{x}_c)}{\prod_{i \in \mathcal{N}(c)} b_i(x_i)} \right) \prod_{\text{nodes } i} b_i(x_i) \qquad (103)$$

where $b_c(\mathbf{x}_c) \equiv \sum_{x \setminus \mathbf{x}_c} b(\mathbf{x})$ (where for instance $\mathbf{x}_c = (x_1, x_2)$) and $b_i(x_i) \equiv \sum_{x \setminus x_i} b(\mathbf{x})$.

The equation above can also be written:

$$b(\mathbf{x}) \approx \prod_{\text{cliques } c} b_c(\mathbf{x}_c) \prod_{\text{nodes } i} b_i(x_i)^{1-|\mathcal{N}(i)|}$$

where $|\mathcal{N}(i)|$ is the number of neighbors of node $i$ in the graph representation of the distribution.

The approximation of $b(\mathbf{x})$ given in Equation (103) is equivalent to approximating the entropy $S_b$ of $b(\mathbf{x})$ as follows:

$$ \begin{aligned} -S_b &= \sum_{\mathbf{x}} b(\mathbf{x}) \log(b(\mathbf{x})) \\ &\approx \sum_{\mathbf{x}} b(\mathbf{x}) \sum_{\text{cliques } c} \log(b_c(\mathbf{x}_c)) + \sum_{\mathbf{x}} b(\mathbf{x}) \sum_{\text{nodes } i} (1 - |\mathcal{N}(i)|) \log(b_i(x_i)) \\ &\approx \sum_{c} \sum_{\mathbf{x}_c} b_c(\mathbf{x}_c) \log(b_c(\mathbf{x}_c)) + \sum_{i} (1 - |\mathcal{N}(i)|) \sum_{\mathbf{x}_i} b_i(x_i) \log(b_i(x_i)) \end{aligned} \qquad (104) $$

In contrast, the Reweighted BP algorithm consists of approximating the variational distribution $b(\mathbf{x})$ as:

$$b(\mathbf{x}) \approx \prod_{\text{cliques } c} \left( \frac{b_c(\mathbf{x}_c)}{\prod_{i \in \mathcal{N}(c)} b_i(x_i)} \right)^{1/\alpha_c} \prod_{\text{nodes } i} b_i(x_i)^{1/\kappa_i} \qquad (105)$$
which can also be written:

$$b(\mathbf{x}) \approx \prod_{\text{cliques } c} b_c(\mathbf{x}_c)^{1/\alpha_c} \prod_{\text{nodes } i} b_i(x_i)^{1/\kappa_i - |\mathcal{N}(i)|/\alpha_i} \quad \text{where} \quad \frac{1}{\alpha_i} = \frac{1}{|\mathcal{N}(i)|} \sum_{c; i \in \mathcal{N}(c)} \frac{1}{\alpha_c} \quad (106)$$

This leads to the following parametric approximation of the variational entropy:

$$\begin{aligned} -S_b &= \sum_{\mathbf{x}} b(\mathbf{x}) \log(b(\mathbf{x})) \\ &\approx \sum_{\mathbf{x}} b(\mathbf{x}) \sum_{\text{cliques } c} \frac{1}{\alpha_c} \log(b_c(\mathbf{x}_c)) + \sum_{\mathbf{x}} b(\mathbf{x}) \sum_{\text{nodes } i} \left( \frac{1}{\kappa_i} - \frac{|\mathcal{N}(i)|}{\alpha_i} \right) \log(b_i(x_i)) \\ &\approx \sum_{\text{cliques } c} \frac{1}{\alpha_c} \sum_{\mathbf{x}_c} b_c(\mathbf{x}_c) \log(b_c(\mathbf{x}_c)) + \sum_{i} \left( \frac{1}{\kappa_i} - \frac{|\mathcal{N}(i)|}{\alpha_i} \right) \sum_{x_i} b_i(x_i) \log(b_i(x_i)) \end{aligned} \quad (107)$$

Additionally, the initial probability distribution $p(\mathbf{x})$ is approximated as follows:

$$p(\mathbf{x}) \approx \prod_{\text{cliques } c} \psi_c(\mathbf{x}_c)^{\beta_c} \prod_{\text{nodes } i} \psi_i(x_i)^{\gamma_i} \quad (108)$$

which is equivalent to the following approximation of the average energy $U_b$:

$$\begin{aligned} U_b &= \sum_{\mathbf{x}} b(\mathbf{x}) E(\mathbf{x}) \\ &\approx -\sum_{\mathbf{x}} b(\mathbf{x}) \sum_{\text{cliques } c} \beta_c \log(\psi_c(\mathbf{x}_c)) - \sum_{\mathbf{x}} b(\mathbf{x}) \sum_{i} \gamma_i \log(\psi_i(x_i)) \\ &\approx -\sum_{\text{cliques } c} \beta_c \sum_{\mathbf{x}_c} b_c(\mathbf{x}_c) \log(\psi_c(\mathbf{x}_c)) - \sum_{i} \gamma_i \sum_{x_i} b_i(x_i) \log(\psi_i(x_i)) \end{aligned} \quad (109)$$

Altogether, we obtain the following approximation of the Gibbs free energy $G_b^{\text{approx}} \approx G_b = U_b - S_b$:

$$\begin{aligned} G_b^{\text{approx}} = & \sum_{\text{cliques } c} \frac{1}{\alpha_c} \sum_{\mathbf{x}_c} b_c(\mathbf{x}_c) \log \left( \frac{b_c(\mathbf{x}_c)}{\prod_{i \in \mathcal{N}(c)} b_i(x_i)} \right) - \sum_{\text{cliques } c} \beta_c \sum_{\mathbf{x}_c} b_c(\mathbf{x}_c) \log (\psi_c(\mathbf{x}_c)) \\ & + \sum_{i} \frac{1}{\kappa_i} \sum_{x_i} b_i(x_i) \log (b_i(x_i)) - \sum_{i} \gamma_i \sum_{x_i} b_i(x_i) \log (\psi_i(x_i)) \end{aligned} \quad (110)$$

## G.5 From Gibbs free energy approximation to messages

In what follows, we derive the message-passing update equations of Reweighted BP (98), (99), and (100) (and its special cases BP, Fractional BP, Power EP and $\alpha$-BP, among others) based on the approximation of the Gibbs free energy $G_b^{\text{approx}}$ given in Equation (110). The message update equations are simply fixed-point equations of $G_b^{\text{approx}}$, meaning that the beliefs computed by Reweighted BP are stationary points of the approximate Gibbs free energy $G_b^{\text{approx}}$. This demonstration is similar to the one for BP [Yedidia et al., 2001, 2003], with additional parameters $\alpha$, $\kappa$, $\beta$, and $\gamma$.

We form the Lagrangian by adding Lagrange multipliers to $G_b^{\text{approx}}$. Lagrange multiplier $\mu_i$ (resp. $\mu_c$) corresponds to the normalization constraint $\sum_{x_i} b_i(x_i) = 1$ (resp. $\sum_{\mathbf{x}_c} b_c(\mathbf{x}_c) = 1$), while $\lambda_{cj}(x_j)$ corresponds to the marginalization constraint $\sum_{\mathbf{x}_{c \setminus j}} b_c(\mathbf{x}_c) = b_j(x_j)$. We obtain:

$$\begin{aligned} \mathcal{L} = G_b^{\text{approx}} &+ \sum_{i} \mu_i \left( \sum_{x_i} b_i(x_i) - 1 \right) + \sum_{\text{cliques } c} \mu_c \left( \sum_{\mathbf{x}_c} b_c(\mathbf{x}_c) - 1 \right) \\ &+ \sum_{\text{cliques } c} \sum_{j \in \mathcal{N}(c)} \sum_{x_j} \lambda_{cj}(x_j) \left( \sum_{\mathbf{x}_{c \setminus j}} b_c(\mathbf{x}_c) - b_j(x_j) \right) \end{aligned} \quad (111)$$
The partial derivatives of the Lagrangian are:

$$ \frac{\partial \mathcal{L}}{\partial b_i(x_i)} = - \sum_{c \in \mathcal{N}(i)} \frac{1}{\alpha_c} + \frac{1}{\kappa_i} + \frac{1}{\kappa_i} \log(b_i(x_i)) - \gamma_i \log(\psi_i(x_i)) + \mu_i - \sum_{c \in \mathcal{N}(i)} \lambda_{ci}(x_i) \quad (112) $$

$$ \frac{\partial \mathcal{L}}{\partial b_c(\mathbf{x}_c)} = \frac{1}{\alpha_c} + \frac{1}{\alpha_c} \log \left( \frac{b_c(\mathbf{x}_c)}{\prod_{i \in \mathcal{N}(c)} b_i(x_i)} \right) - \beta_c \log(\psi_c(\mathbf{x}_c)) + \sum_{i \in \mathcal{N}(c)} \lambda_{ci}(x_i) + \mu_c \quad (113) $$

It comes, by cancelling the partial derivatives of the Lagrangian:

$$ b_i(x_i) \propto \psi_i(x_i)^{\gamma_i \kappa_i} \prod_{c \in \mathcal{N}(i)} \exp(\kappa_i \lambda_{ci}(x_i)) \quad (114) $$

and

$$ b_c(\mathbf{x}_c) \propto \left( \prod_{i \in \mathcal{N}(c)} b_i(x_i) \right) \psi_c(\mathbf{x}_c)^{\alpha_c \beta_c} \prod_{i \in \mathcal{N}(c)} \exp(-\alpha_c \lambda_{ci}(x_i)) \quad (115) $$

$$ \implies b_c(\mathbf{x}_c) \propto \psi_c(\mathbf{x}_c)^{\alpha_c \beta_c} \left( \prod_{i \in \mathcal{N}(c)} \psi_i(x_i)^{\gamma_i \kappa_i} \left( \prod_{d \in \mathcal{N}(i) \setminus c} \exp(\kappa_i \lambda_{di}(x_i)) \right) \exp(\lambda_{ci}(x_i)(\kappa_i - \alpha_c)) \right) \quad (116) $$

We obtain, for $m_{c \to i}(x_i) \equiv \exp(\lambda_{ci}(x_i))$, the following expression of the (approximate) marginal and pairwise beliefs:

$$ b_i(x_i) \propto \psi_i(x_i)^{\gamma_i \kappa_i} \prod_{c \in \mathcal{N}(i)} m_{c \to i}(x_i)^{\kappa_i} \\ b_c(\mathbf{x}_c) \propto \psi_c(\mathbf{x}_c)^{\alpha_c \beta_c} \left( \prod_{i \in \mathcal{N}(c)} \psi_i(x_i)^{\gamma_i \kappa_i} \right) \left( \prod_{i \in \mathcal{N}(c)} \prod_{d \in \mathcal{N}(i) \setminus c} m_{d \to i}(x_i)^{\kappa_i} \right) \left( \prod_{i \in \mathcal{N}(c)} m_{c \to i}(x_i)^{\kappa_i - \alpha_c} \right) \quad (117) $$

(118)

Eventually, thanks to the constraint $\sum_{\mathbf{x}_c \setminus x_j} b_c(\mathbf{x}_c) = b_j(x_j)$, we obtain the fixed point equations for the messages:

$$ m_{c \to j}(x_j)^{\kappa_j} \propto \left( \sum_{\mathbf{x}_c \setminus x_j} \psi_c(\mathbf{x}_c)^{\alpha_c \beta_c} \prod_{i \in \mathcal{N}(c) \setminus j} \left[ \psi_i(x_i)^{\gamma_i} \left( \prod_{d \in \mathcal{N}(i) \setminus c} m_{d \to i}(x_i) \right) m_{c \to i}(x_i)^{1-\alpha_c/\kappa_i} \right]^{\kappa_i} \right) m_{c \to j}(x_j)^{\kappa_j - \alpha_c} \quad (119) $$

$$ \Leftrightarrow m_{c \to j}(x_j) \propto \left( \sum_{\mathbf{x}_c \setminus x_j} \psi_c(\mathbf{x}_c)^{\alpha_c \beta_c} \prod_{i \in \mathcal{N}(c) \setminus j} \left[ \psi_i(x_i)^{\gamma_i} \left( \prod_{d \in \mathcal{N}(i) \setminus c} m_{d \to i}(x_i) \right) m_{c \to i}(x_i)^{1-\alpha_c/\kappa_i} \right]^{\kappa_i} \right)^{1/\alpha_c} \quad (120) $$

The Reweighted BP algorithm consists of running iteratively the fixed-point equation (120). This single equation, which involves only messages from factor node to variable node, can be rewritten into the two following equations by introducing messages from variable node to factor node:

$$ m_{\psi_c \to x_i}^{\text{new}}(x_i) \propto \left( \sum_{\mathbf{x}_c \setminus x_i} \psi_c(\mathbf{x}_c)^{\alpha_c \beta_c} \prod_{x_j \in \mathcal{N}(\psi_c) \setminus x_i} m_{x_j \to \psi_c}(x_j) \right)^{1/\alpha_c} \quad (121) $$

$$ m_{x_j \to \psi_d}^{\text{new}}(x_j) \propto \psi_j(x_j)^{\gamma_j \kappa_j} \left( \prod_{\psi_c \in \mathcal{N}(x_j) \setminus \psi_d} m_{\psi_c \to x_j}(x_j) \right)^{\kappa_j} m_{\psi_d \to x_j}(x_j)^{\kappa_j - \alpha_d} \quad (122) $$
The expression of the unitary $b_i(x_i)$ and clique beliefs $b_c(\mathbf{x}_c)$ comes from Equations (117) and (118):

$$ b_i(x_i) \propto \psi_i(x_i)^{\gamma_i \kappa_i} \prod_{\psi_c \in \mathcal{N}(x_i)} m_{\psi_c \to x_i}(x_i)^{\kappa_i} \quad (123) $$

$$ b_c(\mathbf{x}_c) \propto \psi_c(\mathbf{x}_c)^{\alpha_c \beta_c} \prod_{x_i \in \mathcal{N}(\psi_c)} m_{x_i \to \psi_c}(x_i) \quad (124) $$

Note that one could also use directly Equation (119) instead of (120) to define the message update equation of the Reweighted BP algorithm. In fact, Equations (119) and (120) correspond to the damped versus undamped update equation; see section F about damping. There is no absolute better choice: fixed points obtained are identical in both cases, and damping provides better convergence properties but slows down the system (see section F).

## G.6 Comparison with related models: BP, Fractional BP, Power EP, alpha-BP, and Circular BP

The special case of Belief Propagation is recovered for $(\alpha, \kappa, \beta, \gamma) = (1, 1, 1, 1)$.

Fractional BP [Wiegerinck and Heskes, 2002], Power EP [Minka, 2004] and α-BP [Liu et al., 2019] use the damped message update equation (55) (with $(\kappa, \beta, \gamma) = (1, 1, 1)$) rather than its undamped version (56) (see section F).

Circular BP of Jardri and Denève [2013] also considers $(\kappa, \beta, \gamma) = (1, 1, 1)$, and modifies the message update equation from variable to factor (Equation (96) for BP) into:

$$ m_{x_j \to \psi_d}(x_j) \propto \psi_j(x_j) \left( \prod_{\psi_c \in \mathcal{N}(x_j) \setminus \psi_d} m_{\psi_c \to x_j}(x_j) \right) m_{\psi_d \to x_j}(x_j)^{1-\alpha_{x_j \to \psi_d}} \quad (125) $$

Note that Circular BP was only defined on pairwise factor graphs, meaning that all cliques c have at most 2 elements: the product of Equation (122) does not appear in this case.

Notably, Reweighted BP, BP, Fractional BP and Circular BP are similar but differ in the number of degrees of liberty. BP has no degree of liberty. Fractional BP (as well as Power EP and α-BP) has $n_{factors}$ degrees of liberty. Circular BP has $n_{factors} + n_{variables}$ degrees of liberty. Circular BP has as many degrees of liberty as the number of edges in the factor graph.

**Circular BP on general factor graphs** Similarly to the pairwise factors case, Circular BP can be seen as an approximation of (an extension of) Fractional BP, which includes reweighting of both the variational entropy and the variational average energy components.

As in the pairwise factor case presented in the main text, the Circular BP algorithm on general factor graphs can be seen as an approximation of reweighted BP algorithms, based on the modification of Equation (120) into:

$$ m_{c \to j}(x_j) \propto \sum_{\mathbf{x}_c \setminus x_j} \psi_c(\mathbf{x}_c)^{\beta_c} \prod_{i \in \mathcal{N}(c) \setminus j} \left[ \psi_i(x_i)^{\gamma_i} \left( \prod_{d \in \mathcal{N}(i) \setminus c} m_{d \to i}(x_i) \right) m_{c \to i}(x_i)^{1-\alpha_c / \kappa_i} \right]^{\kappa_i} \quad (126) $$

**Learning to carry out approximate inference in the general case** The experiments reported in section 7 demonstrated successes of supervised learning on rather small and binary graphical models with pairwise interactions. However, many situations where humans perform inference nearly optimally involve much more generative models. A follow-up of this work using the theory developed in the current section could focus on more general graphical models, with higher-order interactions and/or with continuous variables, and investigate whether parameters can be learnt (in a supervised but also unsupervised manner) in these cases.

In particular, Circular BP applied to a probability distribution composed of at most pairwise interactions but any type of variables (discrete or continuous) can be written very simply, similarly to the
binary case (see also Equations 10-12 of the supplementary material of Jardri and Denève [2013]):

$$ M_{i \to j}^{\text{new}}(X_j) = f(B_i(X_i) - \alpha_{ij} M_{j \to i}(X_i), J_{ij}) \quad (127) $$

$$ B_i(X_i) = \kappa_i \left( \sum_{k \in \mathcal{N}(i)} M_{k \to i}(X_i) + \gamma_i M_{\text{ext} \to i}(X_i) \right) \quad (128) $$

with log-messages $M_{i \to j}(X_j) = \log(m_{i \to j}(X_j))$, $M_{\text{ext} \to i}(X_i) = \log(m_{\text{ext} \to i}(X_i))$, and log-beliefs $B_i(X_i) = \log(p(X_i)) + \log(Z)$ where Z is the normalization constant of BP. Finally, f is a function of function (it applies to $g(X_i)$) which, among other things, performs an integration over $X_i$ and returns a function of $X_j$:

$$ f(g(x_i), J_{ij}) = \log \left( \int_{x_i} \psi_{ij}(x_i, x_j)^{\beta_{ij}} \exp(g(x_i)) dx_i \right) \quad (129) $$

Instead of encoding for $\log(p(X = +1)/p(X = -1))$ as in the binary case, the quantity computed is $\log(p(X))$, the log of the whole probability distribution of variable X. Likewise, the external input to the nodes is no longer a scalar ($M_{\text{ext} \to i}(X_i = +1) - M_{\text{ext} \to i}(X_i = -1)$) but a function (that is, $M_{\text{ext} \to i}(X_i)$).