Cyril Furtlehner

Aurélien Decelle

Cycle-based Cluster Variational Method for Direct and Inverse Inference

Cycle-based Cluster Variational Method for Direct and Inverse Inference

Cyril Furtlehner, Aurélien

Introduction

Markov random fields [START_REF] Lauritzen | Graphical models[END_REF] (MRF) are widely used probabilistic models, able to represent multivariate structured data in order to perform inference tasks. They are at the confluence of probability, statistical physics and machine learning [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF]. From the formal probabilistic viewpoint they express the conditional independence properties of a collection of n random variables x = {x 1 , . . . , x n }, in the form of a factorized probability measure, where each factor involves a subset of x. In statistical mechanics the Gibbs measure takes the form of an MRF, to express the thermodynamic equilibrium probability of a system of n interacting degrees of freedom. The practical use of MRF appears also in various applied fields, like image processing, bioinformatics, spatial statistics or information and coding theory. Recent breathtaking successes in artificial intelligence have been obtained by learning deep neural networks whose building blocks are so-called restricted Boltzmann machines i.e. bipartite networks of interacting Ising spins. By stacking them into deep architectures some high level features can be learned recursively [START_REF] Lecun | Deep learning[END_REF] using schematically Monte-Carlo based learning algorithms in combination with the Bragg-Williams mean-field method within a Gibbs-sampling loop. The use of more advanced mean-field methods like the cavity approach could possibly be helpful in this context [START_REF] Gabri É | Training restricted Boltzmann machine via the Thouless-Anderson-Palmer free energy[END_REF]. The main difficulty resides in the fact that these MRF are of practical use in a domain of parameters which clearly corresponds to an ordered phase with strong couplings, which is usually not the most favorable one for applying mean-field methods. Putting aside this potential difficulty, let us simply state the two main generic problems that have to be commonly dealt with when using MRF in practical applications: Direct inference problems:

• computation of marginal probabilities, also called marginalization problem:

p i (x i ) = x\xi P (x i ),
which involves in general an exponential cost with respect to N to be done exactly;

• computing the mode, also referred to as the maximum a posteriori probability (MAP) x = argmax

x P (x), which is generally an NP hard problem [START_REF] Cooper | The computational complexity of probabilistic inference using bayesian belief networks (research note)[END_REF][START_REF] Shimony | Finding MAPs for belief networks is NP-hard[END_REF].

This two problems are of different nature and involve generally distinct techniques which can share sometimes some similarities. The former can be addressed e.g. by Monte-Carlo sampling or by mean-field methods which boils down to some approximation of the entropy contribution to the free energy; the latter is a combinatorial optimization problem which corresponds to the search for the ground state of a system at zero temperature. In this paper we are primarily interested in tackling the first one.

Inverse problem: learning the parameters of the model, given for example by sufficient statistics when the MRF is in the exponential family. For instance the inverse Ising problem [START_REF] In | Efficient structure learning of Markov networks using L 1 -regularization[END_REF][START_REF] Öfling | Estimation of sparse binary pairwise Markov networks using pseudo-likelihood[END_REF][START_REF] Kappen | Efficient learning in Boltzmann machines using linear response theory[END_REF][START_REF] Welling | Approximate inference in Boltzmann machines[END_REF][START_REF] Ézard | Constraint satisfaction problems and neural networks: a statistical physics perspective[END_REF][START_REF] Yasuda | Approximate learning algorithm in Boltzmann machines[END_REF][START_REF] Cocco | Adaptive cluster expansion for the inverse Ising problem: Convergence, algorithm and tests[END_REF][START_REF] Nguyen | Bethe-Peierls approximation and the inverse Ising model[END_REF] consists of finding the set of couplings {J ij } and external fields {h i } of an Ising model

P (s) = 1 Z(h, J) exp i,j J ij s i s j + i h i s i ,
which maximize the associated log likelihood (LL), given data in form of sequences s (k) , k = 1 . . . M or of empirical marginals Ê(s i ), Ê(s i s j ). Generally the partition function Z(h, J) requires an exponential cost with respect to N to be computed exactly.

In order to be useful, any approach based on MRF modeling relies therefore strongly on efficient approximate algorithms, since both direct and inverse problems have potentially an exponential cost in the system size. Belief propagation (BP) and its generalization GBP [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF] have opened the possibility for using MRF in large scale problems even though many restrictions stand in the way of a systematic use, either from convergence problems or from precision issues.

Concerning the direct problem, the situation is quite different when considering the marginalization problem or the MAP problem. Belief propagation algorithms defined for the marginalization problem have a zero temperature counterpart, the max-product or min-sum algorithm which can be used to solve MAP problems. In this context, methods based on refined message passing techniques have been developed rather successfully for addressing some of the shortcomings of belief propagation. First, a treereweighted algorithm with guaranteed convergence, corresponding to solving a problem dual to the linear programming (LP) relaxation of the MAP has been proposed in [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF][START_REF] Wainwright | MAP estimation via agreement on trees: message-passing and linear programming[END_REF], yielding exact solutions on submodular functions [START_REF] Kolmogorov | On the optimality of treereweighted max-product message-passing[END_REF]. Then for tightening LP relaxations, many strategies similar in spirit with GBP, able to scale with large size systems have been proposed. A region pursuit strategy has been set up in [START_REF] Sontag | Tightening LP-relaxations for MAP using message passing[END_REF] in combination with a dual LP message passing solver [START_REF] Globerson | Fixing max-product: Convergent message passing algorithms for MAP LP-relaxations[END_REF], in order to take into account consistency constraints of marginals from higher order clusters of variables. This being limited in practice to small clusters, a strategy able to incorporate frustrated cycles of arbitrary sizes, supposedly responsible for large integral gaps, has been proposed in different contexts [START_REF] Sontag | New outer bounds on the marginal polytope[END_REF][START_REF] Sontag | Efficiently searching for frustrated cycles in MAP inference[END_REF][START_REF] Kudekar | Improved linear programming decoding using frustrated cycles[END_REF].

For the marginalization problem the situation is less favorable. Firstly the use of GBP is hampered by notoriously difficult convergence problems, which have led some authors [START_REF] Yuille | CCCP algorithms to minimize the Bethe and Kikuchi free energies: convergent alternatives to belief propagation[END_REF][START_REF] Heskes | Approximate inference and constrained optimization[END_REF] to consider double loop algorithms, at the price of some computational costs [START_REF] Pelizzola | Cluster variation method in statistical physics and probabilistic graphical models[END_REF]. In addition the choice to be made for region definition is rather open in general, except that a bad choice may lead to poor precision and lack of convergence [START_REF] Welling | On the choice of regions for generalized belief propagation[END_REF], and too large regions are excluded, as computational cost grows exponentially in the size of the largest regions. In particular, feeding GBP with more regions do not guarantee a monotonic increase in precision on the contrary to what the region pursuit in the MAP context is doing. Indeed, with the dual LP setting, the duality bound allows one to directly check whether a new constrained region may improve a solution or not, while there is no such option for GBP. For regular graphs, regions are straightforwardly identified for example with square plaquettes or cells of 2-D and 3-D lattices, as in the Kikuchi cluster variational methods [START_REF] Kikuchi | A theory of cooperative phenomena[END_REF][START_REF] Pelizzola | Cluster variation method in statistical physics and probabilistic graphical models[END_REF] (CVM). But for general graphs, a systematic choice of regions is more difficult to define and also some complexity issues may occur if the size of regions is not controlled. As suggested in [START_REF] Welling | Structured region graphs: Morphing EP into GBP[END_REF] a good choice for the regions to run GBP might be provided by a cycle basis and possibly a weakly fundamental [START_REF] Gelfand | Generalized belief propagation on tree robust structured region graphs[END_REF] cycle basis. An alternative line of research which has also been followed over recent years consists of estimating loop corrections to the Bethe-Peierls approximation in order to improve its accuracy, by addressing directly the errors caused by the presence of loops on multiply connected factor graphs [START_REF] Montanari | How to compute loop corrections to the Bethe approximation[END_REF][START_REF] Chertkov | Loop series for discrete statistical models on graphs[END_REF][START_REF] Parisi | Loop expansion around the Bethe-Peierls approximation for lattice models[END_REF][START_REF] Mooij | Loop corrections for approximate inference on factor graphs[END_REF][START_REF] Xiao | Partition function loop series for a general graphical model: free-energy corrections and message-passing equations[END_REF][START_REF] Ramezanpour | Computing loop corrections by message passing[END_REF][START_REF] Furtlehner | Approximate inverse Ising models close to a Bethe reference point[END_REF]. Note that the frustrated cycle constraints mentioned previously in the MAP context and the loop corrections considered now correspond to two distinct considerations: the first are additional constraints to impose on a collection of pairwise beliefs in order that they can originate from a true probability distribution; the second aims at refining the approximation made on the variational entropy.

In the present work, we investigate further along this direction by generalizing in some way previous considerations [START_REF] Lage-Castellanos | A very fast inference algorithm for finite-dimensional spin glasses: belief propagation on the dual lattice[END_REF][START_REF] Furtlehner | Approximate inverse Ising models close to a Bethe reference point[END_REF] concerning the random Ising model in the absence of local fields. Firstly we analyze in this context convergence problems emerging from canonical definitions of the region graph. This leads us to propose a specific construction of the factor graph, which to some extent solves the convergence issue, as is observed experimentally. Secondly, we exploit a property of the minimizer of the Kikuchi free energy functional associated with certain cycle bases, such that the message to be sent from one region to another can be computed efficiently with help of an internal BP routine to be performed within each (cycle) region, allowing for arbitrary loop sizes to be considered. For binary variables in particular, it is worth exploiting the fact that BP has one single fixed point on a circle [START_REF] Weiss | Correctness of local probability propagation in graphical models with loops[END_REF], and that the loop correction can be computed explicitly on this geometry. These views apply as well to the inverse problem, which consists of learning the model. We show that the aforementioned property of the Kikuchi free energy minimizer can also be exploited, for the inverse Ising problem in particular, in order to learn efficiently the parameters of the model.

The paper is organized as follows: in Section 2 we give a brief introduction of CVM and related GBP algorithms. In Section 3 we specify GBP and the Kikuchi approximation using regions defined with a cycle basis, analyze the Lagrange multiplier structure and propose a mixed region graph, which discards all unnecessary constraints. Section 3.5 details how this framework can be adapted to the maximum a posteriori probability estimation (MAP) context. The problem of choosing a relevant cycle basis is discussed in Section 4. Then Section 5 is devoted to an efficient computation of messages exchanged between cycle and links regions which completes our generalized cycle based belief propagation (GCBP) formulation for direct inference. Some properties of the free energy functional are also discussed at the end of this section. In Section 6 we reverse the equations of Section 5 to address the inverse Ising problem. Finally some numerical tests are presented in Section 7 both for the direct and inverse inference problem.

Cluster variational method and generalized BP

In this Section we give all the necessary material concerning the relation between BP, generalized BP and mean-field approximations in statistical physics. Further details and references can be found e.g. in [START_REF] Pelizzola | Cluster variation method in statistical physics and probabilistic graphical models[END_REF].

Belief propagation and the Bethe approximation

As far as large scale inference is concerned, Pearl's belief propagation [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Network of Plausible Inference[END_REF] and related algorithms constitute central tools in MRF-based inference approaches. The BP algorithm is an iterative algorithm designed to solve a set of fixed point equations. Given an

(x i ) = 1 Z i φ i (x i ) a i m a→i (x i ), b(x a ) = 1 Z a ψ a (x a ) i∈a n i→a (x i ),
where the messages m a→i and n i→a relating factor to variables and variables to factors satisfy the following set of self-consistent equations

m a→i (x i ) = xa\xi ψ a (x a ) j∈a\i n j→a (x j ), (2.1
)

n j→a (x j ) = φ j (x j ) b j\a m b→j (x j ). (2.2)
This algorithm as sketched on Figure 2.1 is exact on a tree, but only approximate on multiply connected factor graphs. When it converges, it does so empirically in O(N log(N )) steps on a sparse random graphs, yielding often rather good approximate marginals.

In [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF] was first established the connection between the BP algorithm of Pearl with a standard mean-field method -the Bethe approximation [START_REF] Bethe | Statistical theory of superlattices[END_REF] -used in statistical physics. As is well known in statistical physics, the Gibbs distribution associated with the energy function E(x) and inverse temperature β, is obtained as a minimizer of the free energy functional of a trial distribution b(x)

βF[b] = βE[b] -S[b] = β x b(x)E(x) + x b(x) log b(x) = -log(Z Gibbs ) + x b(x) log b(x) p Gibbs (x) = -log(Z Gibbs ) + D KL b p Gibbs
as is explicitly seen in the last equality from the non-negativity property of the Kullback Leibler divergence D KL . The mean energy term E[b] can be expressed exactly in terms of marginal distributions obtained from b, like e.g. single and pairwise marginals if E(x) decomposes over pairwise terms. On the other hand, the entropy term S[b] is in general intractable and mean field methods in statistical physics correspond to different ways of approximating this term. The Bethe approximation for instance corresponds to

S[b] ≈ S Bethe def = - i b i (x i ) log b i (x i ) - a b a (x a ) log b a (x a ) i∈a b i (x i ) = i S i + a ∆S a ,
i.e. a sum of individual entropy of each variable, corrected by the mutual information among the group of variables indexed by a. The connection with BP is precisely that a BP fixed point of (2.1,2.2) corresponds to a stationary point of the approximate Bethe free energy complemented with compatibility constraints among marginal probabilities

βF Bethe [b] = βE[b] -S Bethe [b] + a∈F ,i∈a x i λ ai (x i ) b i (x i ) - xa\xi b a (x a )
with help of Lagrange multipliers λ ai (x i ). The BP algorithm actually corresponds to performing the dual optimization with log messages in (2.1,2.2) corresponding to an invertible linear transformation of the Lagrange multipliers,

λ ai (x i ) = log n i→a (x i ) , (2.3
)

log m a→i (x i ) = 1 d i -1 b i λ bi (x i ) -λ ai (x i ), (2.4)
with d i the number of factors containing i. Moreover, as shown in [START_REF] Heskes | Stable fixed points of loopy belief propagation are minima of the Bethe free energy[END_REF] a stable fixed point corresponds to a local minimum of the free energy functional.

The Kikuchi approximation and associated message passing algorithms

In fact as observed in [START_REF] Kikuchi | A theory of cooperative phenomena[END_REF][START_REF] Morita | Cluster variation method and Möbius inversion formula[END_REF], the Bethe approximation is only the first stage of a systematic entropy cumulant expansion over a poset {α} of clusters

S = α ∆S α ,
where ∆S α is the entropy correction delivered by the cluster α with respect to the entropy of all its subclusters. The decomposition is actually valid at the level of each cluster, such that with help of some Möbius inversion formulae, the corrections

∆S β = α⊆β µ(α, β) S α
and subsequently the full entropy can be expressed as a weighted sum

S = α κ α S α of individual cluster entropies S α = - xα b α (x α ) log b α (x α ),
where κ α ∈ Z are a set of counting numbers. For example on the 2D square lattice, one popular Kikuchi approximation amounts to retain as cluster the set of nodes v ∈ V, of links ∈ E and of square plaquettes c ∈ C such that on a periodic lattice the corresponding approximate entropy reads

S = c S c - S + v S v .
In the CVM, the choice of constraints may be arbitrary, as long as the clusters hierarchy is closed under intersection. Once identified, the connection between the Bethe approximation and BP led Yedidia et al. to propose in [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF] a generalization to BP as an algorithmic counterpart to CVM. In fact they introduced a notion of region, relaxing the notion of cluster used in CVM. In their formulation, any region R containing a factor a should contain all variable nodes attached to a in order to be valid. The approximate free energy functional associated with a set of regions is given by

F(b) = R∈R κ R F R (b R ) + R ⊆R x R λ RR (x R ) b R (x R ) - x R \x R b R (x R ) ,
where b R (x R ) and κ R are respectively the marginal probability and counting number associated with region R. The λ RR are again Lagrange multipliers enforcing the constraints among regions beliefs. The only constraint for the counting numbers is that for any variable i or node a

R i κ R = R a κ R = 1.
This ensures the exactness of the mean energy contribution E(b) to the free energy in general as well as the entropy term for uniform distributions in particular. By comparison, there is no freedom in the CVM on the choice of the counting numbers once the set of cluster is given. Additional desirable constraints on the counting numbers are (i) the maxent-normal constraint and (ii) a global unit sum rule for counting numbers,

R∈R κ R = 1.
(2.5)

Condition (i) means that the approximate region based entropy reaches its maximum for the uniform distribution. Condition (ii) insures exactness of the entropy estimate for perfectly correlated distributions. As for belief propagation, a set of compatibility constraints among beliefs are introduced with help of Lagrange multipliers and generalized belief propagation again amounts to solving the dual problem after a suitable linear transformation of Lagrange multipliers hereby defining the messages. Once a fixed point is found a reparameterization property of the joint measure holds:

P (x) ∝ R∈R b R (x R ) κ R .
When the region graph has no cycle, this factorization involves the true marginals probabilities of each region and is exact.

There is some degree of freedom both in the initial choice of Lagrange multipliers and messages leading to different algorithms without changing the free energy and associated variational solutions. A canonical choice is to connect regions only to their direct ancestor or direct child regions leading to the parent-to-child algorithm. With this choice the constraints are however redundant, some linear dependencies are present and this can potentially affect the convergence of the algorithm by adding unnecessary loops in the factor graph. This problem has been addressed in [START_REF] Pakzad | Estimation and marginalization using the Kikuchi approximation methods[END_REF] where for a given region set a construction for a minimal factor graph is proposed.

Main contributions

GBP is a framework corresponding to a wide class of algorithms, which upon a good choice of regions can lead to much accurate results than basic BP. Its systematic use is however made delicate by the following unsolved issues as far as large scale inference is concerned for the marginalization problem:

• there is no automatic and efficient procedure of choosing the regions able to scale with large scale problems for non-regular factor graphs, despite proposals like the region pursuit algorithm [START_REF] Welling | On the choice of regions for generalized belief propagation[END_REF] whose potential use seems however limited to small size systems.

• without special care the computational cost grows exponentially with respect to region size.

• there are difficult convergence problems associated with GBP which have led some to consider double loop algorithms [START_REF] Yuille | CCCP algorithms to minimize the Bethe and Kikuchi free energies: convergent alternatives to belief propagation[END_REF][START_REF] Heskes | Approximate inference and constrained optimization[END_REF] at the price of additional computational burden.

Concerning inverse problems, we are not aware of any method in the family of region based approximation of the log likelihood, going beyond the Bethe approximation at the exception of the exact method proposed in [START_REF] Cocco | Adaptive cluster expansion for the inverse Ising problem: Convergence, algorithm and tests[END_REF], which is however limited to small systems size from the practical point of view. The idea of constructing the region graph from a cycle basis is not new, it is already present as a special case of CVM in [START_REF] Kikuchi | A theory of cooperative phenomena[END_REF] and was first formally proposed in [START_REF] Welling | Structured region graphs: Morphing EP into GBP[END_REF] and refined in [START_REF] Gelfand | Generalized belief propagation on tree robust structured region graphs[END_REF], regarding the choice for the cycle basis, without however explicitly addressing large scale issues listed above. Our contributions in this context is to settle a certain number of technical problems regarding this construction in order to address the above restrictions such that large scale problems can be treated by means of two algorithms GCBP and KIC respectively for direct and inverse pairwise MRF inference. More specifically,

• we address convergence problems by proposing a specific construction of the factor graph in Section 3.4 based on some decomposition of single variable counting numbers unraveled in Section 3.2;

• our construction leads to a linear cost with respect to region size i.e. large cycles, instead of exponential in general as detailed in 5;

• our region graph construction as discussed in Section 4 relies on a minimal cycle basis optimization, which to some extent and thanks to some approximate algorithm can scale-up to relatively large size as seen experimentally in Section 7;

• we propose in Section 6 a general inverse pairwise MRF method based on the Kikuchi approximation which scales linearly with system size, once a cycle basis is given or properly guessed, again without any limitation in the cycles' sizes.

Generalized cycle based BP (GCBP)

The first motivation for attaching regions to the elements of a cycle basis originates in the observation that the Bethe approximation violates the "global unit sum rule" (2.5) for counting numbers, except on singly connected graphs, precisely by an amount corresponding to the cyclomatic number of the graph. Completing the regions set with elements of a cycle basis restores the unit sum rule property [START_REF] Welling | Structured region graphs: Morphing EP into GBP[END_REF].

A different motivation comes from statistical physics considerations associated with the duality transformation [START_REF] Savit | Duality in field theory and statistical systems[END_REF] which can be performed with certain restrictions on the models like e.g. the Ising model without external fields. In such cases, one is naturally led to consider a dual belief propagation on the dual graph whose nodes correspond to the element of a cycle basis [START_REF] Furtlehner | Approximate inverse Ising models close to a Bethe reference point[END_REF]. The extension of such considerations to arbitrary pairwise models led us to consider GBP based on such cycle basis.

Cycle based Kikuchi approximation

To set up notations, we consider a pairwise MRF of n random variables valued in some arbitrary subset x = {x 1 , . . . , x n } ∈ I 1 × . . . I n ⊂ R n , specified by some undirected graph G = (V, E), with vertex set V = {1, . . . , n} and edge set E ⊂ V × V. To simplify we also assume G to be connected. The reference distribution, considered to be pairwise, is of the form

P (x) = ∈E ψ 0 (x ) v∈V φ v (x v ). (3.1)
By definition a cycle of G is an unoriented subgraph where each node has an even degree. The set of cycles is a vector space over Z 2 of dimension |E| -|V| + 1 for a graph with one single component which is assumed from now on. This means that when two cycles are combined, edges are counted modulo 2. Examples of cycle bases are shown on Figure 3.1. For heterogeneous graphs, a simple way to generate a basis consists first in selecting a spanning tree of the graph and associating a cycle with each of the |E| -|V| + 1 remaining links of the graph, by adding to each one the path on the spanning tree joining the two ends of the link. This yields by definition a fundamental cycle basis, associated with the considered spanning tree. Let us assume that a cycle basis of G is given with cycles indexed by c ∈ C = {1, . . . |C|}. |C| = |E| -|V| + 1 also called the cyclomatic number represents the number of independent loops of G. In the Kikuchi CVM approximation that we consider, the maximal clusters are associated with each element of the cycle basis and possibly links which are not contained in any basic cycle. We assume also that one cycle has at most one edge in common with any other cycle. If this is not the case then one edge and one cycle can be added to G in order to restore this property, for each set of cycles having a common group of edges in common (see Figure 3.2). Disconnected intersections can be eliminated by a proper choice of cycle basis. As explained in Section 2 all mean-field type approximations underlying BP or GBP, consist in assuming a factorized form of the joint measure in term of some of its marginal distributions. Within the CVM and given our choice for the maximal cluster this leads to assuming the following factorization of the joint measure:

P GBP (x) = c∈C p c (x c ) ∈E p (x ) κ v∈V p v (x v ) κv , (3.2)
where p c , p and p v are marginal probabilities respectively associated with cycles, links and single variables. As we shall see, and this is an important observation for what follows, the probability p c associated with a cycle can be itself expressed as a pairwise MRF, with each factor corresponding to one edge of the cycle:

p c (x c ) = ∈c ϕ (x ). (3.3)
In (3.2) counting numbers respectively of cycles, edges and vertices are set to

κ c = 1, κ = 1 -d and κ v = 1 -c v κ c - v κ .
d is the number of cycles in C containing edge . This choice is in accordance to general CVM prescriptions, as being dictated by the constraint that each degree of freedom is counted exactly once in the Kikuchi free energy. As already said, thanks to these rules the global unit sum rule for counting numbers is automatically satisfied:

c∈C κ c + ∈E κ + v∈V κ v = |C| -|E| + |V| = 1.
A dual bipartite graph G = (V c , V t , E ) can be defined, where V c indexes the cycle basis, and elements of V t represent connected intersection between cycles, i.e. either single nodes, links or sub-trees corresponding to bridges connecting distant cycles.

Elements of E connect intersecting elements of V c and V t . Under this assumption

Proof. See Appendix A.

The variational problem that GBP aims at solving, is to find the closest distribution of the form (3.2) to the reference distribution (3.1). For later convenience we define and also introduce for any c ∈ C:

ψ (x ) def = ψ 0 (x ) v∈ φ v (x v ),
Ψ c (x c ) def = ∈c ψ (x ) v∈c φ v (x v ). (3.4)
For a candidate measure p, the variational free energy functional reads

F(P GBP ||P ) = c∈C,xc p c (x c ) log p c (x c ) Ψ c (x c ) + ∈E, x κ p (x ) log p (x ) ψ (x ) + v∈V, xv κ v p v (x v ) log p v (x v ) φ v (x v ) + ,c ,x λ c (x ) p (x ) - xc\x p c (x c ) + v, v,xv λ v (x v ) p v (x v ) - x \xv p (x ) + v,c v,xv λ cv (x v ) p v (x v ) - xc\xv p c (x c ) (3.5)
after introducing three sets of Lagrange multipliers, λ c (x ), λ v (x v ) and λ cv (x v ) to enforce respectively cycle-edge, edge-variable and cycle-variable marginals compatibility. The minimum of the free energy is then obtained as:

                   p c (x c ) ∝ Ψ c (x c ) exp ∈c λ c (x ) + v∈c λ cv (x v ) p (x ) ∝ ψ (x ) exp 1 κ v∈ λ v (x v ) - c λ c (x ) p v (x v ) ∝ φ v (x v ) exp - 1 κ v c v λ cv (x v ) + v λ v (x v )
As direct consequence of these expressions we have Corollaire 3.2. p c has the form 3.3.

Single variable counting numbers and dual loops

The counting number κ v contains some information about the local structure of the dual graph. In order to unravel it we define the local dual graph

G v ⊂ G attached to v as G v = (V v;c , V v;t , E v )
, where V v;c are dual vertices corresponding to cycles containing v; V v;t are dual vertices corresponding to all edges containing v with nonzero counting number; E v is the set of dual edges connecting -nodes in

V v;t to their v G G G v
κ v = 1 -d v + C v . (3.6)
Proof. By definition, we have

C v = |E v | -|V v;c | -|V v;t | + d v = v d - c v 1 - v 1 + d v = κ v + d v -1.
where between the first and second line it is remarked that for any parent of v, any c parent of necessarily contains v.

Qualitatively C v represents the number of dual cycles "centered" on v. This decomposition will prove useful for building our cycle based region graph.

Let us give a few examples: for nodes in the bulk of a planar graph we have

C v = 1, on a cubic lattice C v = 3 which generalizes to C v = d(d -1)/2 on a d-dimensional square lattice. On a N/2 + N/2 bipartite graph we have C v = 3N/2 -1 while on a complete graph of size N , using a cycle basis {(1ij), 1 < i < j ≤ N } rooted on node 1, C 1 = (N -2)(N -3)/2 and C v = 0 ∀v = 1.

Parent-to-child algorithm and minimal graphical representation

At this point, following the region-based algorithm [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF] prescriptions, a message passing algorithm can be set-up which rules are associated with the Hasse diagram of the regions hierarchy. Regions are associated with all terms having non-vanishing counting number in (3.2), and directed edges are associated with each Lagrange multiplier added in (3.5), corresponding to direct parent to child relationship, hence discarding the λ cv .

The message rules which are obtained are then based on the existence of a certain linear transformation of the Lagrange multipliers, which allows one to parameterize the beliefs as follows

p v (x v ) = φ v (x v ) v m →v (x v ), p (x ) = ψ (x ) c m c→ (x ) v∈ n v→ (x v ), p c (x c ) = Ψ c (x c ) ∈c n →c (x ) v∈c n v→c (x v ), (3.7)
with

n →c (x ) def = c \c m c → (x ) n v→ (x v ) def = v\ m →v (x v ), n v→c (x v ) def = v, / ∈c m →v (x v ).
From this we get the following message passing rules:

m c→ (x ) v∈ m vc\ →v (x v ) ←- xc\x Ψ c (x c ) ψ (x ) ∈c\ n →c (x ) × v∈c\ n v→c (x v ), (3.8)
m →v (x v ) ←- x \xv ψ (x ) φ v (x v ) c m c→ (x ) v ∈ \v n v → (x v ), (3.9)
where in the first rule the shorthand notation vc\ is used to denote the link in c containing v different from . As noticed in [START_REF] Pakzad | Estimation and marginalization using the Kikuchi approximation methods[END_REF], dependences between Lagrange multipliers are present in the parent-to-child algorithm. This results in a more complex factor graph with more feedback loops than necessary which in turn may cause convergence failures of GBP. As a matter of fact we observe experimentally, both on grids and on heterogeneous graphs tested in Section 7 that the parent-to-child algorithm fails to converge for systems sizes exceeding a few hundreds of nodes whatever damping coefficient is inserted into the message passing equations. In [START_REF] Pakzad | Estimation and marginalization using the Kikuchi approximation methods[END_REF] a minimal graphical representation construction is proposed to settle such problems, in order to eliminate all redundant Lagrange multipliers. In our setting this leads in particular to having any (non-bridge) variable node to be attached to at most one link node and to have therefore at most one ancestor cycle node in the factor graph. As a consequence we have always

n v→c (x v ) = m vc\ →v (x v ) = 1.
As shown in Appendix B this leads to an essentially unstable algorithm for graphs containing at least one single dual loop. So in short we have

• poor global convergence properties of the parent to child algorithm;

• local convergence problems for the minimal region graph based algorithm caused by dual loops.

This problem of redundant Lagrange multipliers has actually also been discussed in the context of the 2-D Edward Anderson (EA) model in [START_REF] Domínguez | Characterizing and improving generalized belief propagation algorithms on the 2d Edwards-Anderson model[END_REF]. In this context the authors propose a solution based on a specific gauge choice for the message definition in order to regularize GBP. Our approach to this problem is different. As we shall see in the next Section it is solely based on topological properties of the graph of interactions. This yields a generic method independent of the graph or the type of interactions.

Mixed factor graph and associated message passing rules

We introduce here a specification of the region graph which on the one hand eliminates all unnecessary feed-back loops present in the parent-to-child algorithm, but on the other hand prevents instabilities associated with dual loops. In this formulation first a minimal set of Lagrange multipliers are taken into account as proposed in [START_REF] Pakzad | Estimation and marginalization using the Kikuchi approximation methods[END_REF]; but additional "clone variables" need to be introduced for variables at the center of dual loops, i.e. for which C v = 0, as defined in Section 3.2, to prevent some instability which we have identified (see Appendix B). Before explaining it in detail let us give the specification of the region graph which we refer to as the mixed factor graph (MFG) for reasons which will soon be clear: -(a) if d v > 1 one v-node is associated with v, which connects exactly to one single arbitrary -node of each components of G v , its degree being therefore d v and a counting number of 1 -d v is attributed to it. If necessary an -node with zero counting number can be inserted into the MFG in order to ensure that this v-node is properly connected to all components it needs to be.

• (i) Each
-(b) if C v > 0,
to each containing v we associate one v -node that is singly connected to as long as this -node is in a component of G v containing at least one dual loops. Each clone is attributed a counting number

κ v = C v /q if q is the number of clones.
This set of rules is illustrated on Figure 3.5. Rule (iii) ensures that all marginal probabilities of cycles are compatibles at link intersections. Rule (iv)(a) is applied to cutvertices, i.e. vertices which separate G in multiple components when removed as shown on the example of Figure 3.5. Rule (iv)(b) is there to take into account dual loop corrections. The prescription (iv)(b) is there to ensure a better convergence of GCBP by making use of replicas of v-nodes, while preserving the minimal use of Lagrange multipliers. The number of constraints is still minimal in the sense that the number of independent loops of the MFG is equal to the number of independent loops of the dual graph G . From the Lagrangian formulation κ v is constrained by

v ≈v κ v = C v
where ≈ indicates the correspondence between v -node and variable v. The choice made in rule (iv)b for κ vs satisfies this constraint, albeit other ones are possible. The reason for introducing clone variables becomes clearer when trying to write down message passing equations. In fact a direct generalization of the change of variables (2.3,2.4) used to define ordinary BP from the Lagrange multipliers can be ob-tained as follows:

λ v (x v ) = log v\ m →v (x v ) def = log n v→ (x v ), λ v (x v ) = -κ v log m →v (x v ) def = log n v → (x v ), λ c (x ) = log n →c (x ) + v∈ log n v→ (x v ),
where v∈ is taken over all types of v-nodes and with

n →c (x ) def = c \c m c → (x )
Note that λ cv have disappeared by definition of the MFG. We get the following expression for the beliefs

p v (x v ) = φ v (x v ) exp - 1 1 -d v v λ v (x v ) = φ v (x v ) v m →v (x v ), p v (x v ) = φ v (x v ) exp - 1 κ v λ v v (x v ) = φ v (x v )m v →v (x v ), p (x ) = ψ (x ) exp 1 κ v∈ λ v (x v ) - c λ c (x ) = ψ (x ) c m c→ (x ) v∈ n v→ (x v ), p c (x c ) = Ψ c (x c ) exp ∈c λ c (x ) = Ψ c (x c ) ∈c n →c (x ) v∈ n v→ (x v ) , (3.10)
where v denotes the -node connected to v . From this we get the following message passing rules:

m c→ (x ) ←- xc\x Ψ c (x c ) ψ (x ) ∈c\ n →c (x ) v∈ n v→ (x v ) , (3.11)
m →v (x v ) ←- x \xv ψ (x ) φ v (x v ) × c m c→ (x ) v ∈ \v n v → (x v ), (3.12
)

m →v (x v ) ←- x \xv ψ (x ) φ v (x v ) × c m c→ (x ) v ∈ \v n v → (x v ) 1/(1+κ v )
.

(3.13)

The difference between factor graph of the standard parent-to-child algorithm, the minimal one proposed in [START_REF] Pakzad | Estimation and marginalization using the Kikuchi approximation methods[END_REF] and the one associated with MFG is illustrated on Figure 3.6.

3 v 2 1 c3 c1 c2 3 c1 c2 c3 2 3 v 3 1 v 1 v 2 Parent to Child v 2 1 c3 c1 c2
Minimal Region graph Mixed Factor Graph With this formulation GCBP can be seen mainly as an ordinary belief propagation defined on the MFG, where (3.11,3.12) are direct generalization on a MFG of ordinary BP update rules (2.1,2.2), with an additional peculiarity given by dual loop corrections carried by clone variables in (3.13).

MAP estimation

The general inference schema proposed in the previous sections can be straightforwardly adapted to the optimization context, the same way as the min-sum algorithm also called belief revision [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Network of Plausible Inference[END_REF] is derived from BP, by simply replacing " " by "min" (see e.g. [START_REF] Ruozzi | Message Passing Algorithms for Optimization[END_REF]) . First, adding some specific notations, the messages are parameterized in terms of log probability ratio:

m c→ (x ) ∝ exp -µ c→ (x ) and m →v (x v ) ∝ exp -µ →v (x v ) .
The counterparts to "n" messages are in turn defined as:

ν →c (x ) def = c \c µ c → (x ), ν v→ (x v ) def = v\ µ →v (x v ), ν →c (x v ) def = -κ v µ →v (x v ),
where again clone variables are distinguished from ordinary ones using notation. Correspondingly, let

E c (x c ) def = -log Ψ c (x c ) , E (x ) def = -log ψ (x ) and E v (x v ) def = -log φ v (x v ) .
To the generalized belief propagation rules (3.11,3.12,3.13) correspond the following min-sum update rules:

µ c→ (x ) ←-min xc\x E c (x c ) -E (x ) + ∈c\ ν →c (x ) + v∈ ν v→ (x v ) , (3.14)
µ →v (x v ) ←-min x \xv E (x ) -E v (x v ) + c µ c→ (x ) + v ∈ \v ν v → (x v ) , (3.15)
µ →v (x v ) ←- 1 1 + κ v min x \xv E (x ) -E v (x v ) + c µ c→ (x ) + v ∈ \v ν v → (x v ) . (3.16)
As a result the beliefs associated with the various family of nodes, expressing log marginal probabilities, are given by

E v (x v ) = E v (x v ) + v µ →v (x v ), E (x ) = E (x ) + c µ c→ (x ) + v∈ ν v→ (x v ), E c (x c ) = E c (x c ) + ∈c ν →c (x ) + v∈ ν v→ (x v ) . (3.17)
When the messages correspond to a fixed point, the usual compatibility between beliefs is expressed as

min xc\x E c (x c ) = E (x ), ∀ ∈ c, min x \xv E (x ) = E v (x v ), ∀v ∈ .
In addition, if the joint probability measure is given in a Gibbs form,

P (x) = e -E(x) ,
these beliefs provide us with the following decomposition, up to a constant, of the energy function :

E(x) = c E c (x c ) + κ E (x ) + v κ v E v (x v ),
and the approximate minimizer of E(x), given by

x min i = argmin xi E i (x i ), ∀i ∈ V, verifies E(x min ) = c min xc [E c (x c )] + κ min x [E (x )] + v κ v min xv [E v (x v )],
by virtue of the belief's compatibility. Next, as will be also the case for inference, we exploit the ring geometry in order to compute efficiently the c-node to -node messages 3.14. This can be done in O(nq 3 ) time complexity per message. Indeed, the c-node to -node message update simply reads:

µ c→ (x ) = min xc\x E c (x c ) -E (x ) -ν →c (x ) - v∈ µ →v (x v ).
Running a min-sum algorithm associated with the energy function E c (x c ) given x on the loop c for each ∈ c yields immediately µ c→ .

Cycle basis determination 4.1 Various criteria

At this point, nothing has been said concerning the choice of the cycle basis. In [START_REF] Gelfand | Generalized belief propagation on tree robust structured region graphs[END_REF] it is argued that a good choice of basis ensures the algorithm of being tree-robust (TR), namely that GBP converges to an exact fixed point when the underlying graph G is singly connected. They provide a characterization for cycle bases ensuring this property. First it has to be a weak fundamental cycle basis (WFCB), ensuring in particular the maxent property to be satisfied. By definition a cycle basis is fundamental if each cycle contains an edge that is not included in any other basis cycle. For a WFCB, this constraint is relaxed, it is a cycle basis for which there is an ordering such that each cycle contains a link which is absent of all preceding cycles in this ordering. In addition the WFCB is TR, if it is such that any subset of the cycle basis contains a set of links, each one pertaining to a unique cycle in this subset, and altogether forming at least one loop. The reason behind this can be understood quite simply in the special context of CVM approximation (3.2), where a simple reduction rule as the ones given in [START_REF] Welling | On the choice of regions for generalized belief propagation[END_REF] is at work. Suppose the MRF is such that the set of non trivial links ψ

(0) ij (x i , x j ) = f (x i )g(x j
) in (3.1) forms a tree T .

Proposition 4.1. (i) if a trivial link pertains to a single cycle, the factorized joint measure (3.2) coincides with the same CVM approximation defined on a reduced graph, where link has been removed and c is discarded.

(ii) if the cycle basis is a WFCB based on a series of trivial links, the factorized joint measure (3.2) is reduced to the Bethe joint measure associated with the underlying tree T .

Proof. (ii) is the direct consequence of (i) by induction. See Appendix C.

As already stated in Proposition 3.1, GCBP is exact when the dual graph G and henceforth the MFG are acyclic. It could be tempting to push the logic to the end and try to impose a "dual-tree robust" condition for the cycle basis, i.e. that GCBP be exact if there exists a cycle basis of G such that G be singly connected. Clearly this is a dead end, as can already be seen by considering the simple example of a planar graph: the natural cycle basis given by the faces of the graph cannot fulfill such property, when all links at the border of the graph are non-trivial. Nevertheless, let us simply notice that in the case where the underlying graph of non trivial links noted T 2 has an acyclic dual graph T 2 , we have the following Proposition 4.2. GCBP will converge to the exact fixed point if (i) the cycle basis has a subset which is a cycle basis of T 2 , (ii) the complementary set of cycles defines a graph for which it is a WFCB based on trivial links.

Proof. The argument is the same as before, applying the reduction property (i) of the preceding Proposition to the complementary set of cycles, until reaching the core sub-graph T 2 , for which GCBP is exact.

TR cycle bases are easily identified in special cases like planar or complete graphs [START_REF] Gelfand | Generalized belief propagation on tree robust structured region graphs[END_REF], but searching for such a basis in general is difficult, its existence being not always guaranteed. Instead there is yet another feature that could be even more desirable, namely that the cycle basis be such that the number of independent dual cycles, i.e. the cyclomatic number of G , be minimal. Recall that GCBP is similar to an ordinary BP on the MFG. Consequently, as for an ordinary BP, we expect these (dual) loops to be a source of problems. As observed in [START_REF] Furtlehner | Approximate inverse Ising models close to a Bethe reference point[END_REF], the dual cyclomatic number depends on the sum of cycle sizes noted |c|:

C(G ) = C(G) c=1 |c| -C(G) -|E| + P(G ),
with P(G ) the number of connected components of G . As a result, a good choice for the cycle basis could be the minimal cycle basis (MCB) and there exists polynomial time algorithms for finding it [START_REF] Horton | A polynomial-time algorithm to find the shortest cycle basis of a graph[END_REF]. Furthermore if one wants to remain close to the TR prescription, one could even search for a minimal WFCB, which is an APX-hard problem but for which efficient heuristic do exist [START_REF] Rizzi | Minimum weakly fundamental cycle bases are hard to find[END_REF].

Heuristic algorithm

Exact algorithms for solving the MCB problem have a polynomial time complexity, scaling typically like O(N L 2 ) up to logarithmic corrections [START_REF] Kavitha | Cycle bases in graphs characterization, algorithms, complexity, and applications[END_REF]. Making use of these would spoil the efficiency of GCBP, whose main virtue is to scale linearly with systems size. We have therefore to resort to some approximate procedure. It is guided by the empirical assumption that most important loops to be taken care of are the smallest ones. The main steps of the method are the following:

• (i) build a subset of candidate cycles which contains the most important ones.

This step can be made linear with system size for sparse graphs with bounded degree d max ; typically O(N d n max ) for finding cycles with sizes ≤ n.

• (ii) complete this set in order to have a complete set containing the MCB. This step can be done exactly in O(N L) time complexity [START_REF] Kavitha | Cycle bases in graphs characterization, algorithms, complexity, and applications[END_REF].

• (iii) Extract an independent set of shortest sizes. Exact methods typically use Gaussian elimination which is the main source of time consumption.

This strategy is basically the one which is followed by the most efficient exact algorithms. In order not to be a limiting speed factor for GCBP steps (ii) and (iii) have to be approximated. Note that step (ii) is not mandatory. Since the goal is to take into account most important loop corrections, an independent set of short cycles, not necessarily complete can be used. Concerning step (iii) we replace the Gaussian elimination procedure by an approximate one which has the additional virtue of respecting as much as possible the WFCB criteria explained in the previous section. Our algorithm goes as follows:

• (S0) Initialization: weight all the links with the number n of cycles in the candidate set they belong to and extract with respect to these weights a maximum spanning tree from G called G 0 . Create a double ordered list {c 0 (n, s)} of candidate cycles indexed by their number n of links not already present in G 0 and their sizes s. Create an empty list of cycle elements B 0 .

• (S1) cycle selection: At step t select in c t the cycle c with smallest n and then with smallest size s and update B t+1 ←-B t + {c}.

• (S2) update (c t , G t ) -→ (c t+1 , G t+1 ):
if n = 1: insert the corresponding link into G t to obtain G t+1 and update c t in c t+1 . All cycles with n = 0 have a linear decomposition in B t+1 and are eliminated.

if n > 1: insert one of the n free links of c into G t to obtain G t+1 . Update c t in c t+1 as if all the n links where selected. For each of the n -1 nonselected links of c create a new cycle by joining this link to the path on G t connecting its two ends point, using the Dijkstra algorithm 1 . Insert these new cycles into c t+1 .

if c t+1 = ∅ go back to (S1) else exit().

Note that if by chance the new added cycle at each step corresponds to n = 1 we would get a WFCB. The procedure followed in the case n > 1 is there to ensure that we get a complete set at the end. As already said this is not mandatory in practice, so if this constraint is relaxed, then the n links can be directly inserted into G t .

Cycle basis cleaning

Once a cycle basis has been obtained some adjustments have to be performed to cope with GCBP. First the basis can be optimized further by a local random greedy shuffling procedure, which consists in looking for combination of pairs of cycles sharing some links, from which smaller cycles can be generated (see Figure 4.1). Secondly, S = 5 + 6 S = 5 + 5 as already stated in the MFG prescriptions any pair of cycles must have at most one single link in common. Note in passing that this requirement seems actually difficult if not impossible in general to conciliate with the search for TR cycle basis advocated in [START_REF] Gelfand | Generalized belief propagation on tree robust structured region graphs[END_REF]. In contrast, the smaller the aggregated cycle's size is, the less cleaning is to be expected. By cleaning we mean the operation shown on Figure 3.2. This consists in to add one link relating the two ends of a path common to two or more cycles and formed by at least two links. In this operation a new cycle composed of this path and of the new added link is created which, when combined with the other cycles containing that path leaves all these cycles intersecting on this single link. This cleaning operation is done greedily by treating in order the intersection paths with largest sizes until intersections only composed of one single link remain. Finally in some cases, cycles remain which have non-connected intersections with other cycles. This kind of situation occur sometimes but rarely, so in practice the adopted cleaning procedure consists simply to discard the largest cycle involved in such pathological intersections.

As we observed in practice, these cleaning operations take a small if not negligible part in the overall computation time needed to determine the cycle basis. The complete workflow is shown on the example of Figure 4.2 leading to the MFG starting from a bipartite graph. 5 Loop corrections: c-node to -node messages

General case

We exploit now the specific structure of the cycle-based region definition to propose an efficient method (see Figure 5.1) for computing the messages (3.11), with a cost at most linear in the size of the cycles per message. c-node to -node messages are computed using,

m c→ (x ) = p c (x ) ψ (x )n →c (x ) v∈ n v→ (x v ) , (5.1)
where

p c (x ) def = xc\x p c (x c ).
is the pairwise marginal associated with any link ∈ c, obtained from distribution (3.10). We wish to bypass the summation over x c \x , which has an exponential cost in the size of the loop. Variables x ∈ {1, . . . q} are assumed to have q possible states and p c is a product of pairwise factors along the cycle

p c (x c ) = ∈c ψ c (x ).
On the ring geometry, the partition function as well as any correlation function can be expressed as the trace of a product of transition matrices:

Z ring = Tr n =1 M ( ) ,
where M ( ) is a q 2 matrix with elements given by M ( ) xy = ψ c (x, y).

Upon introducing the following matrices

U def = n i=1 M (i) , U (i) def = n j=i M (j) i-1 j=1 M (j) , V (i) def = n j=i+1 M (j) i-1 j=1 M (j) ,
the expression for the exact marginals are given by

p c i (x) = 1 Z ring Tr δ xx U (i) , p c i (x, y) = 1 Z ring Tr δ xy V (i) .
In this form the cost for computing each c-node to -node message is O nq 3 ). As shown in [START_REF] Weiss | Correctness of local probability propagation in graphical models with loops[END_REF], running BP on a single cycle always converges and there is a linear relation between single variable beliefs and the exact marginals given by the largest eigenvalue of some product of matrices taken from the factors along the loop. In fact, somewhat simpler relations can be established, valid also for pairwise marginals, by applying to a single loop the general loop corrections [START_REF] Chertkov | Loop series for discrete statistical models on graphs[END_REF][START_REF] Sudderth | Loop series and Bethe variational bounds in attractive graphical models[END_REF] expansion to BP. First factorize p c (x c ) with help of BP,

p c (x c ) = 1 Z BP n i=1 b c i (x i , x i+1 ) b c i (x i ) (5.2)
by means of a set of single and pairwise beliefs b c i (x i ) and b c i (x i , x i+1 ), where i = 1, . . . n indexes the variables along the cycle. We define the following q 2 matrices in operator form:

B (i) xy def = b c i (x, y) -b c i (x)b c i+1 (y) b c i (x)
, and associated product of matrices

U def = n i=1 B i , U (i) def = n j=i B j i-1 j=1 B j , (5.3)
V (i) def = n j=i+1 B j i-1 j=1 B j .
Proposition 5.1. The relations between beliefs and exact marginals are then given by

p c i (x) = b c i (x) + U (i) xx Z BP with Z BP = 1 + Tr(U ) p c i (x, y) = b c i (x, y) + V (i) yx b i (x) + B (i) xy V (i) yx Z BP Proof. See Appendix D for details.
The c-nodes messages 3.11 can then be computed from these exact marginals. From these expressions, we see that the cost for computing each message is still O nq 3 ). The only benefit of using the BP factorization at this point resides in the fact that B (j) and therefore U (i) and V (i) have an obvious zero eigenmode:

y B (j) xy b c i (y) = 0.
Trying to find the other modes is not advantageous in general except if some symmetries are present or when q is small. In particular for the binary case (q = 2) we end up with a scalar problem for expressing loop corrections, as is detailed in the next section.

Binary case

For binary variables this relationship can be made even more explicit as we show now. Using of the standard Ising spin notation, each node i ∈ 0, . . . n-1 is associated with a binary variable s i ∈ {-1, 1} and the joint measure of s def = {s 1 , . . . , s n } is exponential and given by

P c (s) = 1 Z c exp n i=1 h c i s i + n-1 i=1 J c i s i s i+1 , (5.4)
where h c i ∈ R is the local field exerted on variable i and J c i ∈ R denotes the coupling between s i and s i+1 . Running BP on this measure leads to the following factorization of the joint measure:

P (s) = 1 Z BP n i=1 b c i (s i , s i+1 ) b c i (s i )b c i+1 (s i+1 ) n i=1 b c i (s i ), (5.5)
where the b c i (•) and b c i (•, •) are the single and pairwise approximate marginals delivered by BP. These can be parameterized as follows

b c i (s i ) = 1 2 (1 + mi s i ), (5.6) b c i (s i , s i+1 ) = 1 4 (1 + mi s i + mj s j + ( mi mj + χi )s i s j ), (5.7)
where m i def = E(s i ) represents the "magnetization" of spin s i and

χ i def = E(s i s i+1 ) - E(s i )E(s i+1
) is the covariance, also named "susceptibility" coefficient, between s i and s i+1 . We use the sign˘to denote a BP estimate, which is to be distinguished from the exact value. The relation between BP values and exact ones can be made explicit in the following form.

Proposition 5.2. Let Q def = n i=1 χi » (1 -m2 i )(1 -m2 i+1 ) , (5.8)
then the BP normalization constant, the exact magnetization and susceptibility coefficients read:

Z BP = 1 + Q, (5.9
)

m i = 1 -Q 1 + Q mi (5.10
)

χ i = χi 1 + Q + Q 1 + Q (1 -m2 i )(1 -m2 i+1 ) χi + 4 mi mi+1 1 + Q . (5.11)
Proof. The proof is based on the following identity

b i (s i , s i+1 ) b i (s i )b i+1 (s i+1 ) = 1 + χi (s i -mi )(s i+1 -mi+1 ) (1 -m2 i )(1 -m2 i+1 )
, and follows the same lines as the proof of Proposition 5.1.

Section 6 will be based on these identities. The corresponding loop corrected marginals p i and p ii+1 are expressed from the loop corrected quantities (m i , m i+1 , χ i ) through the same relations (5.6) and (5.7) and allow one to obtain all messages 3.11 sent by the c-node at once from the BP beliefs, so the cost per-message in this special case is now O(1) instead of O(n) if there are n messages to be sent.

In addition to this slight but non-crucial reduction in computational cost, one should note the scalar characterization in terms of Q ∈] -1, 1] of the cycle which shows up. First from the matrix formulation 5.8, Q is the non-zero eigenvalue of U . It is the product of "BP correlations" along the loop and characterizes its strength.

• Q 0 corresponds to weak loop correction, BP is nearly exact.

• Q → 1 corresponds to a strongly correlated loop.

• Q → -1 corresponds to a strongly correlated frustrated loop.

Loop corrections to the Bethe Free Energy

The formalism used previously suggests reconsidering the cycle based Kikuchi approximate free energy by rewriting it in an appealing form where loop corrections are made more explicit. Indeed using the BP factorization of each independent cycle marginal (5.2) yields the following decomposition of the entropy term for any pairwise MRF in terms of single and pairwise marginals {p i , i ∈ V} and {p , ∈ E} and associated cycle beliefs {b c i , (i, c) ∈ V × C} and {b c , ( , c) ∈ E × C}. Starting from the cluster expansion we have:

S Kikuchi = i∈V S i + ∈E ∆S + c∈C ∆S c .
The first two terms represent the Bethe entropy,

S Bethe = i S i + ∆S ,
as a sum of individual variables entropy S i corrected by mutual information of variables

-∆S = x p (x ) log p (x ) p 1 (x 1 )p 2 (x 2 )
≥ 0, counted for each link ∈ E. The corrections induced by each cycle c has the following expression:

∆S c = S c - i∈c S i - ∈c ∆S = log(Z c BP ) - i∈c D KL (p i b c i ) + ∈c D KL (p b c ), (5.12) = F Bethe p c p c ],
where Z c BP is the normalizing factor of the BP factorization (5.2) associated with the cycle marginal distribution p c compatible with the p i 's and p 's. The cycle beliefs b c i and b c are implicitly and uniquely determined from the p 's. F Bethe is the Bethe approximation to the free energy functional:

F p p 0 ] = D KL (p p 0 ) + F 0 ,
F 0 being the free energy associated with p 0 . This has the following immediate consequence. Let us consider an auxiliary measure, build from the exact marginals:

pc (x c ) def = 1 Zc BP ∈c p (x ) i∈c p i (x i )
with normalization constant Zc BP .

Lemma 5.3.

log Z c BP ≤ ∆S c ≤ log Zc BP .
(5.13)

Proof. Recall that on the loop geometry BP has one single stable fixed point which corresponds to a global minimum of the approximate Bethe free energy functional [START_REF] Heskes | Stable fixed points of loopy belief propagation are minima of the Bethe free energy[END_REF]. Consequently, the minimum for fixed b is obtained for p i = b i and p = b in (5.12)

F Bethe p c p c ] ≥ log(Z c BP ),
which proves the left hand side inequality. Next consider the following quantity:

D KL (p c pc ) = log Zc BP Z c BP + i∈c D KL (p i b c i ) - ∈c D KL (p b c ) = log Zc BP -∆S c ≥ 0,
since the Kullback-Liebler divergence is non-negative, we get the right hand side inequality of (5.13).

As a consequence of (5.13), if the stochastic operator U defined by ( 5.3) has a positive trace then the loop correction has a counter effect to the Bethe correction ∆S . In particular for binary variables in the ferromagnetic case, log(Z c BP ) = log 1 + Q c with Q c ≥ 0, leading therefore to negative loop corrections to the Bethe free energy. Since the Kikuchi correction is exact in absence of dual loops, i.e. when C i = 0, ∀i ∈ V, we may expect that the correction is overestimated in presence of dual loops, i.e. that we should have a bounding of the free energy

F Kikuchi ≤ F ≤ F Bethe , (5.14)
for ferromagnetic like systems, when F Bethe and F Kikuchi are given in terms of the exact single and pairwise beliefs {p i , i ∈ V} and {p , ∈ E}. Note that the inequality F ≤ F Bethe only proved in some special ferromagnetic cases [START_REF] Sudderth | Loop series and Bethe variational bounds in attractive graphical models[END_REF], involves the approximate marginals given by BP instead of the exact ones in our case. The conditions under which the bounding (5.14) might be relevant is left aside to future investigations. All this also suggests that in presence of dual loops some appropriate correction terms proportional to local dual loop counting numbers C v could be inserted into the free energy functional in order to compensate for the kind of "overcounting" of loop corrections which occurs in such cases. This possibility which would potentially lead to a new family of approximate and hopefully more precise mean field schema is left as a side remark for the moment and will be investigated in the near future.

Kikuchi cycle-based (KIC) inverse inference

From the explicit expression of the Kikuchi type approximation (3.2) it should be in principle possible to find a set of fields and couplings corresponding to a given input of single and pairwise empirical marginals. Assuming first we know the graph structure and have a cycle basis, it remains to determine the marginal probabilities p c , p and p v associated with each region. We expect the p 's and p v 's to be given from the data, but the p c 's have to be constructed. This means that the global inverse problem is decomposed into |C| small inverse problems. In the Ising case, if we denote by h c i and J c the local field and coupling associated as in (5.4) with the marginal representing cycle c, ĥ i , Ĵ associated with p and finally ĥi to p i , then from (3.2) the corresponding Kikuchi cycle based (KIC) approximate inverse Ising solution reads

h (KIC) i = κ i ĥi + c i h c i + i (1 -d ) ĥ i , J (KIC) = (1 -d ) Ĵ + c J c .
When the graph structure is unknown, one possibility is to select a set of candidate links, the one carrying the largest amount of mutual empirical information among all possible edges. Then on the graph defined by those links an algorithm is run in order to find the minimal cycle basis, with respect to the weights given by minus the mutual information. More refined strategies could then be used like the one based on iterative proportional scaling proposed in [START_REF] Martin | GMRF Estimation under Topological and Spectral Constraints[END_REF] in the context of Gaussian MRF.

In the following we concentrate on how to invert equations (5.10,5.11) in order to compute h c i and J c for any cycle c ∈ |C|.

Fixed point method

Consider a single loop of size n. Assume we are given a set of empirical marginals pi (s i ) and pi (s i , s i+1 ), for i = 1, . . . n or equivalently a set of magnetizations mi and susceptibilities χi . First note that the change of variables {h i , J i , i = 1, . . . n} to { mi , χi } is a one to one mapping: on the one hand h i and J i can be explicitly written in terms of the { mi , χi } (see below); on the other hand, on a loop there is a unique BP fixed point yielding factorization (5.5), so through relations (5.6,5.7) { mi , χi } are uniquely determined. Finding a joint-measure of highest likelihood to model the empirical marginals is therefore equivalent to find a set of parameters mi and χi defining the joint-measure (5.5) which satisfy χ i = χi and m i = mi in equations (5.10,5.11). The problem is therefore to find the unique value of Q for which all these relations are satisfied. Note also that these relations could be as well obtained by writing down the gradient of the log likelihood, which in the (h, J) variables is a convex function. Hence these equations must anyway have a unique valid solution. The reason for not working in these (h, J) variables is that the LL is not given explicitly in these variables but in the m and χ variables (see below). By rewriting equations (5.10,5.11) in terms of the spin-spin correlation

Θ i def = χ i » (1 -m 2 i )(1 -m 2 i+1 ) , (6.1)
letting Q simply read

Q = n i=1 Θi , (6.2)
we arrive at the following fixed-point equation:

Proposition 6.1. The solution ( m, χ) satisfying equations (5.10,5.11) for a given set {m i = mi def = tanh( ĥi ), i = 1, . . . n} and {χ i = χi , i = 1, . . . n} of empirical magnetizations and susceptibilities is determined by the n-dimensional vector Θ obeying

Θ = f ( Θ), with f i ( Θ) def = A i (Q) Θi - Q Θi , (6.3)
where

A i (Q) def = (1 + Q)(1 -Q) 2 Θi -4Q(1 + Q) sinh( ĥi ) sinh( ĥi+1 ) » (1 -2Q cosh( ĥi ) + Q 2 )(1 -2Q cosh( ĥi+1 ) + Q 2 ) , (6.4)
Proof. Expressing all the magnetizations mi in equation (5.11), in terms of Q and tanh( ĥi ) with help of (5.10), after performing the change of variables χi -→ Θi yields the desired result.

Let us specify the domain D ⊂ [-1, 1] n of validity for this iteration schema. For arbitrary magnetizations and susceptibilities there are some basic constraints. The first one is that mi ∈ [-1, 1], for all i ∈ {1, . . . , n} which entails

Q ≤ Q max def = max i 1 -mi 1 + mi .
The second set of constraints is that probabilities b(s i , s i+1 ) are in [0, 1]:

∀(s i , s i+1 ) ∈ {-1, 1} 2 , 0 ≤ (1 + mi s i )(1 + mi+1 s i+1 ) + χi s i s i+1 ≤ 4. (6.5)
We may rewrite these constraints in a more convenient form. We denote by hi the local fields corresponding to mi = tanh( hi ). In these notations the constraints now read:

0 ≤ e hisi+ hi+1si+1 + Θi s i s i+1 ≤ 4 cosh( hi ) cosh( hi+1 ).
Considering all possible cases for (s i , s j ) we end up with the following somewhat simpler constraints:

-e -| hi+ hi+1| ≤ Θi ≤ e -| hi-hi+1| , (6.6)
which combined with Q ∈ [-1, Q max ] entirely defines the domain D and which prove useful in practice to restrict efficiently the search for a fixed point in a valid domain.

Stability analysis: In order to remain inside D the iteration schema is defined as follows:

g :D -→ D (6.7) X -→ Y =    f X , if f X ∈ D, U (D), if f X / ∈ D. (6.8)
where f coincide with (6.3) for any χ such the image is in the domain D and is otherwise replaced by a random function U : D -→ D. This function consists first of drawing Q uniformly between ] -1, Q max ], and then drawing Θi for each i = 1 . . . n, uniformly between the bounds given in (6.6). Finally an overall scaling is applied to each Θi if the product exceeds Q max . Defined as it is g is an iterative map on a compact domain with no other guaranty than there exists one unique fixed point solution.

Let us examine the conditions under which this solution corresponds to a stable fixed point. The Jacobian of this iterative map, when it coincides with f reads

J ij def = ∂f i ∂Θ j = Q Θ j A i (Q) -(1 -δ ij ) 1 Θ i .
Denoting Θ 1,2 min the two lowest absolute values of Θ i and

B(Q) def = max i |A i (Q)Θ i |,
we get the following sufficient condition of local convergence: Proposition 6.2. The fixed point is stable in general if

|Q| < Θ (1) min Θ
(2)

min n -1 + B(Q) , (6.9)
and in particular if

|Q| < Θ (1) min Θ (2) min n . (6.10)
in absence of magnetization.

Proof. See Appendix E

When some of the magnetizations mi are non zero, the coefficient B(Q) can become arbitrarily large when Q approaches Q max so clearly there exists a value of |Q| above which the condition 6.9 will be violated. For small Q we have

B(0) = max i Θi -4 sinh( ĥi ) sinh( ĥi+1 ) + cosh( ĥi ) + cosh( ĥi+1 ) , (6.11)
which as well diverges when one of the magnetization mi approaches ±1, which means that convergence problems are likely to occur in this domain. Instead, for small magnetizations B(Q) can get smaller to 1, lim maxi mi→0

B(Q) = max i Θi ≤ 1.
The inequality (6.10) becomes relevant in this regime and the iterative schema can converge for small Q, in particular if the largest correlation Θ is not greater than n -1/(n-2) which is close to 1 for n 1.

Line search optimization

The preceding conditions are not always met to guarantee the convergence of the fixed point method. Therefore we develop an alternative method which directly maximizes the log likelihood, this latter being an explicit function LL( Θ) of the Θ i 's,

LL( Θ) def = -log 1+Q( Θ) + i w i ( Θ)+h i ( Θ) mi +J i ( Θ)( χi + mi mi+1 ) (6.12)
By convention we have

LL( Θ) = -∞, ∀ Θ / ∈ D.
The corresponding Ising fields and couplings of the cycle are given by

w i = 1 4 log b i (-1, -1)b i (-1, 1)b i (1, -1)b i (1, 1) b 2 i (-1)b 2 i (1)
h i = 1 2 log b i (-1) b i (1) + 1 4 j∈{i-1,i} log b j (1, 1)b j (s i = 1, s j = -1) b j (s i = -1, s j = 1)b j (-1, -1) J i = 1 4 log b i (-1, -1)b i (1, 1) b i (-1, 1)b i (1, -1) ,
in addition to the weighting exponents w i which show up. All these parameters are given through (5.6,5.7) as function of the magnetizations mi and susceptibilities χi which in turn are fully determined by the Θi 's through (6.1) and (5.10,6.2) given m i = mi . Let D Q ⊂ [-1, Q max ] denote the domain of possible values for Q. In order to find the optimal point we show the following Proposition 6.3. There exists two functions

h : D Q -→ R Θ : D Q -→ D such that argmax Θ∈D LL( Θ) = Θ(Q ) with Q = argmax Q∈D Q h(Q).
Proof. To prove this we explicitly construct these functions, which in turn will be used to run a line search algorithm.

First note that taking the gradient of LL( Θ) with respect to the mi 's and χi 's in order to find the stationary points leads to equations (5.10) and (5.11). After doing the change of variables and manipulations given in Proposition 6.1, the set of equations to be solved reads:

Θ2 i -A(Q) Θi + Q = 0, for i = 1, . . . n,
where Q depends implicitly on the solutions. A first consequence is that, given Q, there is the constraint that the quadratic equation have solutions, i.e. that

A i (Q) 2 -4Q ≥ 0, ∀ i = 1, . . . n,
which depends only on the empirical values mi and χi . This further constraints the domain

D Q ⊂ [-1, Q max ] of possible values of Q.
If this condition is fulfilled, for each i = 1, . . . n, there are two solutions,

Θi (Q, σ i ) = A(Q) + σ i A(Q) 2 -4Q 2 ,
where σ i ∈ {-1, 1} is introduced by convenience. Unfortunately, in general both solutions can be valid, as long as they satisfy the constraints (6.6). At the fixed point, which is unique, the mi 's and Θi 's are uniquely given by Q, therefore among the 2 n possible choices, the correct one will satisfy (6.2) and corresponds to the lowest likelihood. The function h can now be defined as follows:

h : D Q -→ R Q -→ LL Θ(Q)
where Θ(Q) in turn is given as

Θ(Q) = argmax σ LL Θ (Q, σ) (6.13) with Θ i (Q, σ i ) = Q n j=1 Θj (Q, σ j ) Θi (Q, σ i ).
This last normalization is there to ensure that Θ(Q) effectively corresponds to Q.

Combined method and MRF inference

The two methods can be combined by selecting the solution with highest LL (6.12), after running each one with a fixed computational budget. The line search method has a combinatorial step present in (6.13), which can be solved by simple enumeration for small loops, but may become problematic for large ones, n 1. However, for larger cycles, already typically for n > 5, Q is usually very small and the iterative schema of Section 6.1 is converging. Even though some specific optimization might well be possibly developed to solve (6.13), we leave this question aside, as being non-critical from what is seen experimentally on Figure 6.1.

To infer an MRF, a set of candidate cycles is either given either pre-processed from the data e.g. using mutual information scores. As already mentioned, in such case we look for a minimal cycle basis, which in practice, can be approximately obtained at low computational cost as in experiments of the next Section, by a simple stochastic heuristic of loop mixing. For general pairwise MRF, with non-binary variables no specific method is proposed at the cycle level, but at least a gradient descent could be used to solve each cycle independently. If necessary, a posterior selection procedure, based on the generated solution, could be used to refine the cycle basis, with various possible heuristics. Concerning the overall computational cost needed to generate an approximate MRF solution, assuming a "low-cost" method for fixing the cycle basis, it is linear in the number of candidate cycles i.e. in the number of potential links. Therefore the method can in principle cope with large scale problems when a sparse graph is to be expected.

Experiments

We have run various experiments to see how this approach to direct and inverse inference works in practice. With a damping factor up to .5 inserted in the c-node to -node messages needed at low temperature, GCBP always converge on these small grids instances to a fixed point corresponding to a paramagnetic state. At larger scale Figure 7.1, thanks again to a damping factor up to .6, the algorithm is also always converging on the considered range of temperature and sizes but two dynamical regimes are observed. At high temperature, for β ≤ 1.5 the computational time grows like N α with a slight departure from linear complexity as β increases, α = 1.05 for β = 0.5 and α = 1.15 at β = 1.5.

Direct inference

In that case all the fixed points correspond to paramagnetic states. Instead at β > 1.5 and no external fields, the occurrence of non-paramagnetic states is observed at sufficiently large scale, N ≥ 10 5 for β > 1.5 and N ≥ 10 4 for β = 2, as observed also in the ±J 2-D EA model2 in [START_REF] Domínguez | Characterizing and improving generalized belief propagation algorithms on the 2d Edwards-Anderson model[END_REF]. This is an artifact of the Kikuchi approximation since the 2-D EA model is thought to be exempt from a spin-glass phase [START_REF] Örg | Strong universality and algebraic scaling in two-dimensional Ising spin glasses[END_REF]. Convergence is still observed in this regime, but huge fluctuations in computational time occur, depending on whether GCBP converges towards a paramagnetic or to a spin-glass fixed point. On the example shown, outliers points with respect to the fitted scaling actually correspond to spin-glass fixed points, while all other points are paramagnetic. This is clearly related to the fact that a long range order has to be found by a GCBP fixed point when converging to a spin-glass state which is not the case for a paramagnetic one. Indeed in the paramagnetic situation, fixed point messages depend on others within distances on the grid of the order of the spatial characteristic scale for the correlations which increases with β. When compared to BP, the computational time for GCBP is larger by a factor of 5 to 25, but in addition to being less precise, BP is by far less robust and actually stops converging around β 1. The same experiments are performed first on small random sparse 20 + 20 regular bipartite graphs, for which exact beliefs can as well be computed by complete enumeration. In these cases the cycle basis is not given in advance and has to be determined. On Figures 6.2 0.002 × x 1.25 0.0015 × x 1. 55 3.10 -7 × x 2 bipartite graphs of size up to N = 10 5 and mean connectivity up to d = 6. We obtain as well good convergence properties, with no convergence failures, thanks again to a damping factor of .7 for d = 3 to .9 for d = 6. Concerning computational time we observe a scaling in N α which deviates from the linear one as expected as the graph becomes denser for a fixed temperature, α ranging from 1.2 at d = 3 to 1.55 at d = 6.

Heterogeneous graphs with larger mean connectivity have a tendency to contains more highly connected nodes for which C v 1. We suspect these nodes to be mainly responsible for a slowing down of convergence. On the same figure we also show the computational time needed by our approximate pre-processing cycle basis stochastic optimization. The scaling is quadratic when the heuristic detailed in Section 4 is used in its complete version, but the very small multiplicative constant allows us to go for relatively large size, before becoming a limiting factor for GCBP around N 10 4 for d = 3 and N 10 5 for d = 4. Since collecting most important small loops has a linear complexity, the way to overcome this issue at large scale is then to limit ourselves to an incomplete set of independent cycles.

Inverse inference

For the inverse Ising problem, we first test the single loop algorithm explained in Section. 6.1 and Section. 6.2 and the results are shown on Figure 6.1. For this we generate loops of increasing sizes S ∈ {3, . . . 8}. Couplings and biases are sampled as before, with an inverse temperature parameter β varied again in the range [0, 5]. The inference is considered successful for a precision threshold, arbitrarily chosen to 10 -5 β, on the max error of the couplings and biases. A comparable computational budget of a maximum of 100 iterations for FP or estimations for LS is given to both methods. Note however that generally when it converges FP does it within 10 or 20 iterations. The fixed point method is always successful for all sizes when β ≤ 1.2, but this rate degrades when β is increased albeit less severely with larger loops. On the contrary, the line search method is not sufficiently precise at small β but sees its success rate increasing with β especially for small loops. Therefore the two methods are very much complementary, and combining them leads to nearly maximal success rates, at least for β ≤ 3.

Our KIC method is then tested and compared with the linear response of the Bethe-Peierls approximation [START_REF] Nguyen | Bethe-Peierls approximation and the inverse Ising model[END_REF] (BA+LR) at infinite sampling and with the pseudo-likelihood method (PLM) [START_REF] Ravikumar | Highdimensional Ising model selection using L 1 -regularized logistic regression[END_REF][START_REF] Decelle | Pseudolikelihood decimation algorithm improving the inference of the interaction network in a general class of Ising models[END_REF] at finite sampling, again on small square grid and on small sparse random bipartite models. Couplings and biases are sampled as before. Comparison with BA+LR indicates a gain in precision between 1 to 2 orders of magnitude for 5 × 5 grids as seen on Figure 7.3 (left). For bipartite models, Figure 7.3 (right) shows a decreasing gain with increasing mean connectivity, BA+LR and KIC returning the same error around d = 3.4. On Figure 7.4 one representative grid and bipartite instances are shown. As expected the error increases with β but stays reasonably close to the order of a few percents in the strong coupling region β > 1, in contrary to BA+LR which is useless in this region. At finite sampling, by comparing with PLM, we see that the precision is either limited by the sampling itself (small β or small sampling N s ≤ 10 5 ) either by the Kikuchi approximation itself for β > 1 and N s = 10 6 on the grid instance and at N s = 10 4 and β > 1 on the bipartite instance.

Conclusion

Our investigations on GBP has led us to propose a systematic way of dealing with cycle regions and a new mean field approach to inverse problems. Our contribution is two-fold: for the direct problem, we propose (i) an original specification of the region graph (MFG) ensuring simple and robust convergence properties (ii) the loop message computation using ordinary BP ensuring fast message exchange between regions. (i)+(ii) characterize GCBP as a new region based algorithm generic to pairwise MRF, which we have made specific in the binary case. For the inverse Ising problem, we propose a new mean-field approach (KIC) general for pairwise MRF models, which is simple and efficient at least for binary models and sparse graphs without necessarily finite tree-width like 2-d grids. In particular the modular aspect of the method, which consists in a decomposition of the problem into small independent inverse problems corresponding to each independent cycle is valid in general, not only for binary MRF. For incomplete data, since it takes as input single and pairwise marginals, it could be a good alternative to PLM which requires instead complete data.

Still, the scalability of GCBP and KIC relies on the scalability of the cycle basis search algorithm for irregular graphs. In [START_REF] Gelfand | Generalized belief propagation on tree robust structured region graphs[END_REF] it is argued that a good choice of basis ensures the algorithm of being tree-robust (TR), namely that GBP converges to an exact fixed point when the underlying graph G is singly connected after eliminating fake links. In our experiments we did not follow this prescription, but instead proposed a simpler one, namely based on the search for a minimal cycle basis, for which a specific heuristic has been developed with reasonable scalability.

Concerning possible applications of this work, it is planned to use both the direct and inverse approach in combination, in order to test some traffic prediction schema based on the Ising model that has been developed in some preceding related work [START_REF] Martin | Latent binary MRF for online reconstruction of large scale systems[END_REF]. In addition, the systematic treatment of the loops that we propose could presumably be extended in a specific way to the Potts model which has been applied in many different contexts like image processing [START_REF] Tanaka | Statistical-mechanical approach to image processing[END_REF] for instance. Yet another perspective of this framework is to be found in the combinatorial optimization context which could help improve approximate heuristics.

A Proof of Proposition 3.1

If G is acyclic, we can build a junction tree using each cycle as a clique, so the form 3.1 is correct except maybe for the specific form chosen for p c . The leaf nodes of G correspond either to dandling trees or to cycle regions of the primal graph G. From the hypothesis on G these components are connected to the rest of the primal graph G either via a single node or via a link. So summing over all variables contained in each of these regions except the contact node or link results in a subgraph of G whose dual is still acyclic, with a modified factor corresponding to the contact link or vertex. By induction, G can be reduced until one single arbitrary loop region remains, which still corresponds to a sub-graph of G. This results therefore in a marginal probability p c having pairwise form with factor graph corresponding to cycle c.

B Dual loop-based instabilities

Let us consider an Ising model on the single dual loop graph of Figure 3.6 with uniform external field h and coupling J. We give the label 0 to the central node with counting number κ 0 = 1 and labels {1, 2, 3} to the peripheral ones, these having κ v = 0. Links with non-vanishing counting numbers (κ = -1) are for ∈ {01, 02, 03}, cycles are labelled {012, 023, 031}. Using the corresponding minimal factor graph, we attach arbitrarily the only v-node, indexed by 0, to = 01. The following exponential parameterization of the messages is adopted: m c→ (s ) = e w c→ +h 1 c→ s 1 +h 2 c→ s 2 +J c→ s 1 s 2 m →0 (s 0 ) = e w →0 +h →0 s0 .

From the update rules (3.8,3.9) we get in particular for (i, j) ∈ {(1, 2), (2, 3), (3, 1)} m 0ij→0i (s 0 ) ←sj exp h 0 0kj→0j s 0 + h j + h j 0kj→0j s j + J 0j + J 0kj→0j s 0 s j , and more specifically

h 0 0ij→0j ←-h 0 0kj→0j + 1 4 log A ++ A -+ A +-A -- with A σ1σ2 def = h 0 0kj→0j + σ 1 (h j + h j 0kj→0j ) + σ 2 (J 0j + J 0kj→0j
). From this we see that these iterative equations are at least marginally unstable, by the presence of an eigenmode of the Jacobian of eigenvalue 1 corresponding to h 0 0kj→0j = cte, ∀kj. One additional dual loop centered on v-node 0 would actually render this mode unstable.

with b i (x i , x j ) b i (x i )b i+1 (x i+1 ) = 1 + b i (x i , x i+1 ) -b i (x i )b i+1 (x i+1 ) b i (x i )b i+1 (x i+1 ) def = 1 + B (i) xixi+1 b i+1 (x i+1 ) ,
and then by expanding the factors when taking averages. Let us call bond ii + 1 the contribution corresponding to the factor

B (i)
x i x i+1 bi+1(xi+1) instead of 1. The point is that one extremity of a bond cannot be left alone in this expansion, if the corresponding variable is summed over, because of the following identities:

xi b i (x i ) B (i) xixi+1 b i+1 (x i+1 ) = xi B (i-1)
xi-1xi = 0.

For the partition function for instance, either all or none of the bound have to be selected, yielding only the two contributions:

Z BP = x n i=1 b i (x i ) + n i=1 B (i) xixi+1 , = 1 + Tr(U ).
For the single variable marginal, say p i (x i ), again either none or either all of the bonds have to be selected, giving For the pairwise marginals p i (x i , x i+1 ) two additional contributions emerge corresponding to selecting only the bond ii + 1 or to selecting all the bonds except this one, yielding the announced expression.

p i (x i ) = 1 Z BP x\xi

E Proof of Proposition 6.2

The problem is to bound in absolute value the largest eigenvalue of the Jacobian. Let λ be an eigenvalue and v be an eigenvector of J Let v = max j v j 47 and i the corresponding index, such that v i = v. We have

|λ| = | j J ij v j v | ≤ j |J ij | ≤ | Q Θ j A i (Q)| + j =i | Q Θ i Θ j | ≤ |Q| Θ (1)
min Θ

(2) min

B(Q) + n -1 ,
with the definition of B(Q) and Θ

(1,2)

min given in the text. Imposing |λ| ≤ 1 leads to the conditions given in the proposition. In particular when magnetizations are absent, i.e. when h i = 0, ∀i, we have

A (Q) = Θi so B(Q) = max i | Θi Θ i | ≤ 1.

Figure 2 . 1 :

 21 Figure 2.1: Factor graph and message propagation.

7 Figure 3 . 1 :

 731 Figure 3.1: Example of cycle bases on 2-D and 3-D lattices and a fundamental cycle basis on an arbitrary graph.

Figure 3 . 2 :

 32 Figure 3.2: Dual graph construction. Dashed link correspond to one virtual added link. we have the following important property, illustrated on Figure 3.3 which justifies the approximation (3.2,3.3). Proposition 3.1. If G is acyclic, the factorization 3.2 is exact.

Figure 3 . 3 :

 33 Figure 3.3: Successive graphical models obtained by deconditioning variables (circled in red) from the leaves, starting from a polygon tree. Factors corresponding to links or vertices in red are modified during the process.

Figure 3 . 4 : 4 Proposition 3 . 3 .

 34433 Figure 3.4: Local dual graph construction. In this case the choice of cycle basis leads to κ v = 2 with d v = 3 and C v = 4. corresponding c-nodes in V v;c they belong to in the primal graph. This construction is illustrated on Figure 3.4 Proposition 3.3. Let d v be the number of components of G v and C v its cyclomatic number. We have κ v = 1 -d v + C v . (3.6)

Figure 3 . 5 :

 35 Figure 3.5: Pairwise MRF (left). Variables and links with non-zero counting number are in bold. Corresponding mixed factor graph (right) with counting numbers.

3 Figure 3 . 6 :

 336 Figure 3.6: One dual loop on top (C v = 1) with corresponding factor-graphs.

Figure 4 . 1 :

 41 Figure 4.1: Example of cycle combinations leading to smaller cycle basis.

Figure 4 . 2 :

 42 Figure 4.2: Example of 7 + 7 regular bipartite graph of mean connectivity 3.4, and corresponding mixed factor graph, with c-nodes, -nodes and v -nodes colored respectively in red, blue and green. v-nodes associated with bridges are absent on this example. 4 auxiliary links (in red on the middle panel) have been inserted in order to ensure single link intersection between cycles as explained in Section 4.3.

Figure 5 . 1 :

 51 Figure 5.1: Message exchange at the cycle level.

Figure 6 . 1 :

 61 Figure 6.1: Success rates for the inverse inference on a single cycle with different sizes (color) for the fixed point (FP), the line search (LS) and the combined methods (LS+FP).

Figure 6 . 2 :

 62 Figure 6.2: Mean error for the direct inference of 2-D random Ising model comparing GCBP with BP as a function of β, on a 5 × 5 square grid (left) and on random 20 + 20 bipartite graphs of mean connectivity 4 (right) with or without local fields of amplitude 0.2β, averaged over 100 instances.

Figure 6 . 3 :

 63 Figure 6.3: Success rates and mean error for the direct inference problem, comparing GCBP with BP on 20+20 random bipartite graphs of mean connectivity 5 with varying β (left) or with fixed β = 1 and increasing the mean connectivity (right), in presence or not of random local fields of max amplitude 0.2β, averaged over 100 instances.

Figure 6 .

 6 Figure6.2 deals with direct inference, GCBP is run on 5 × 5 grids so that the RMSE on the beliefs (single and pairwise) can be computed by exact enumeration. Couplings J ij and local fields h i are i.i.d sampled uniformly respectively in the range [-β, β] and [-0.2β, 0.2β] when local fields are present. β is varied on the range [0, 5], so that weak and strong couplings are tested. 100 instances are generated for each point. With a damping factor up to .5 inserted in the c-node to -node messages needed at low temperature, GCBP always converge on these small grids instances to a fixed point corresponding to a paramagnetic state. At larger scale Figure7.1, thanks again to a damping factor up to .6, the algorithm is also always converging on the considered range of temperature and sizes but two dynamical regimes are observed. At high temperature, for β ≤ 1.5 the computational time grows like N α with a slight departure from linear complexity as β increases, α = 1.05 for β = 0.5 and α = 1.15 at β = 1.5.In that case all the fixed points correspond to paramagnetic states. Instead at β > 1.5 and no external fields, the occurrence of non-paramagnetic states is observed at sufficiently large scale, N ≥ 10 5 for β > 1.5 and N ≥ 10 4 for β = 2, as observed also in

2

 2 GCBP β = 1.5 GCBP β = 1 GCBP β = 0.5BP β = 1 BP β = 0.5 6.10 -4 × N 1.23 6.10 -4 × N 1.15 5.10 -4 × N 1.07 2.510 -4 × N1.05

Figure 7 . 1 :

 71 Figure 7.1: Convergence behaviour of GCBP and BP regarding computational time on 2-D EA models of large sizes. Cases corresponding to β = 0.5, 1 have local random fields in [-0.1β, 0.1β] while other cases are without external fields.

3 GCBP d = 4 GCBP d = 6 MCB d = 3 MCB d = 4

 34634 .b and 6.3.a we again vary the temperature for a fixed mean connectivity d = 4 and d = 5, while on Figure 6.3.b the inverse temperature is kept fixed at β = 1 and the mean connectivity is varied up to d = 9. As seen on Figure 6.3.b. convergence problems are absent below β 2 but occur at small temperatures with increasing frequency above this threshold signaling the presence of a spin glass phase. In addition, up to d = 9 we observe a significant gain factor in the error made by GCBP compared to ordinary BP. On Figure 7.2 are shown results of tests that were performed on random sparse MCB d = 6 5.510 -4 × x 1.24

Figure 7 . 2 :

 72 Figure 7.2: Computational times of GCBP and the MCB search algorithm on random bipartite graphs at β = 1 for different mean connectivity d.

Figure 7 . 3 :

 73 Figure 7.3: Comparison of KIC with BA+LR at infinite sampling on a 5×5 square grid when β is varied (left), on random bipartite graphs at β = 1 with biases of amplitude 0.2β varying the mean connectivity (right).

Figure 7 . 4 :

 74 Figure 7.4: Comparison of KIC with BA+LR at infinite and with PLM at finite sampling on a 5 × 5 square grid (left) and on a bipartite model with connectivity 3(right) when β is varied.

  term in (3.2) having a non-zero counting number is associated with a node in the MFG. There are three families of nodes, c-nodes, -nodes and v-nodes, respectively associated with cycles, links and vertices of the original graph. c-nodes are always factors while v-nodes are always variables. Instead, -nodes associated with links are composite nodes, i.e. can be of both types.• (ii) Edges of the MFG represent Lagrange multipliers and relate variables to factors. A v-node can be linked to -nodes, considered then as factors nodes. -nodes considered as variable nodes can be linked to c-nodes. • (iii) all links of a given cycle c with non-vanishing counting numbers are linked as variables to this c-node. • (iv) to a variable v we associate in general two types of v-nodes depending on d v and C v defined in Section 3.2:

This ensures the independence of these new cycles among each others and with B t+1

Thresholds are comparable after dividing our β by √

to have random models with identical variance of the couplings.

C Proof of Proposition 4.1

By definition of the Lagrange multipliers, when a fixed point is obtained, the corresponding set of beliefs {b i , b , b c } allows one to factorize the joint measure as (3.1), where for all cycles of the basis, b

where the b c i and b c ii+1 are obtained from b c by running BP on the cycle and are in general different from the b i and b computed globally. The relation between the two corresponds to the loop correction. Let us call trivial, an edge (ij) whose factor is trivial

Similarly we say that a cycle has a trivial belief if it is related to variable and pairwise beliefs as

i.e. the b i and b c i coincide. First we remark that a cycle c containing one such trivial edge, not contained in any other cycle, has necessarily a trivial belief, because from the factorization (3.1) for any edge we have in that case

so the pairwise cycle belief has to be of the form b c (x ) = b c i (x i )b c j (x j ). As a result the factorized joint measure actually coincides with the same CVM approximation form (3.2) on a reduced graph, where link has been removed and c is now discarded. From hypothesis (ii) the set of trivial links contained in one single cycle is non empty. As a results all these link can be removed and all corresponding cycles discarded. On the reduced graph, again since all cycles have a trivial belief, there is a non-empty subset of trivial links, that can be removed and so on. The procedure stops after eliminating all trivial links until only the underlying dual tree remains. The definition of the counting numbers ensures that we then end up with the Bethe form of the joint measure associated with this dual tree.

D Proof of Proposition 5.1

The proof is based on the following factorization of the joint measure on a cycle with help of a belief propagation fixed point: