Cyril Furtlehner

Yufei Han

Jean-Marc Lasgouttes

Victorin Martin

Pairwise MRF Calibration by Perturbation of the Bethe Reference Point

Keywords: Random Markov Fields, Ising Models, Inference, inverse problem, mean-field, belief propagation) Champs Markovien aléatoires, modèle d'Ising, Inférence, problème inverse, champ moyen, propagation de croyances

We investigate different ways of generating approximate solutions to the inverse problem of pairwise Markov random field (MRF) model learning. We focus mainly on the inverse Ising problem, but discuss also the somewhat related inverse Gaussian problem. In both cases, the belief propagation algorithm can be used in closed form to perform inference tasks. Our approach consists in taking the Bethe mean-field solution as a reference point for further perturbation procedures. We remark indeed that both the natural gradient and the best link to be added to a maximum spanning tree (MST) solution can be computed analytically. These observations open the way to many possible algorithms, able to find approximate sparse solutions compatible with belief propagation inference procedures. Experimental tests on various datasets with refined L 0 or L 1 regularization procedures indicate that this approach may be a competitive and useful alternative to existing ones.

Introduction

The problem at stake is a model selection problem, in the MRF families, where N variables are observed pair by pair. The optimal solution is the MRF with maximal entropy obeying moment constraints. For binary variables, this happens then to be the Ising model with highest log-likelihood. It is a difficult problem, where both the graph structure and the value of the fields and coupling have to be found. In addition, we wish to ensure that the model is compatible with the fast inference algorithm "belief propagation" (BP) to be useful at large scale for real-time inference tasks. This leads us to look at least for a good trade-off between likelihood and sparsity.

Concerning the Inverse Ising Problem (IIP), the existing approaches fall mainly in the following categories:

• Purely computational efficient approaches rely on various optimization schemes of the log likelihood [START_REF] Lee | Efficient structure learning of Markov networks using L 1 -regularization[END_REF] or on pseudo-likelihood [START_REF] Höfling | Estimation of sparse binary pairwise Markov networks using pseudo-likelihood[END_REF] along with sparsity constraints to select the only relevant features.

• Common analytical approaches are based on the Plefka expansion [START_REF] Plefka | Convergence condition of the TAP equation for the infinite-ranged Ising spin glass model[END_REF] of the Gibbs free energy by making the assumption that the coupling constants J ij are small. The picture is then of a weakly correlated uni-modal probability measure. For example, the recent approach proposed in [START_REF] Cocco | Adaptive cluster expansion for the inverse Ising problem: convergence, algorithm and tests[END_REF] is based on this assumption.

• Another possibility is to assume that relevant coupling J ij have locally a tree-like structure. The Bethe approximation [START_REF] Yedidia | Generalized belief propagation[END_REF] can then be used with possibly loop corrections. Again this corresponds to having a weakly correlated uni-modal probability measure and these kind of approaches are referred as pseudo-moment matching methods in the literature for the reason explained in the previous section. For example the approaches proposed in [START_REF] Kappen | Efficient learning in Boltzmann machines using linear response theory[END_REF][START_REF] Welling | Approximate inference in Boltzmann machines[END_REF]28,[START_REF] Yasuda | Approximate learning algorithm in Boltzmann machines[END_REF] are based on this assumption.

• In the case where a multi-modal distribution is expected, then a model with many attraction basins is to be found and Hopfield-like models [START_REF] Hopfield | Neural network and physical systems with emergent collective computational abilities[END_REF][START_REF] Cocco | High-dimensional inference with the generalized Hopfield model: Principal component analysis and corrections[END_REF] are likely to be more relevant. To be mentioned also is a recent mean-field methods [START_REF] Nguyen | Mean-field theory for the inverse Ising problem at low temperatures[END_REF] which allows one to find in some simple cases the Ising couplings of a low temperature model, i.e. displaying multiple probabilistic modes.

On the side of inverse Gaussian problem, not surprisingly similar methods have been developed by explicit performing L 0 and L 1 matrix norm penalizations on the inverse covariance matrices, so as to determine sparse non-zero couplings in estimated inverse covariance matrices for large-scale statistical inference applications [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Hsieh | Sparse inverse covariance matrix estimation using quadratic approximation[END_REF] where direct inversion is not amenable. In our context the goal is a bit different. In general cases, the underlying inverse covariance matrix is not necessarily sparse. What we aim to find is a good sparse approximation to the exact inverse covariance matrix. Furthermore, sparsity constraint is not enough for constructing graph structure that is used in conjunction with BP, known sufficient conditions referred as walk-summability [START_REF] Malioutov | Walk-sums and belief propagation in Gaussian graphical models[END_REF] (WS) are likely to be imposed instead of (or in addition to) the sparsity constraint. To the best of our knowledge not much work taking this point into consideration at the noticeable exception of [START_REF] Anandkumar | High-dimensional Gaussian graphical model selection: Walk summability and local separation criterion[END_REF] by restricting the class of learned graph structures. To complete this overview, let us mention also that some authors proposed information based structure learning methods [START_REF] Netrapalli | Greedy learning of Markov network structure[END_REF] quite in line with some approaches to be discussed in the present paper.

In some preceding work dealing with a road traffic inference application, with large scale and real time specifications [START_REF] Furtlehner | A belief propagation approach to traffic prediction using probe vehicles[END_REF][START_REF] Furtlehner | Learning multiple belief propagation fixed points for real time inference[END_REF][START_REF] Furtlehner | Spatial and temporal analysis of traffic states on large scale networks[END_REF], we have noticed that these methods could not be used blindly without drastic adjustment, in particular to be compatible with belief propagation. This led us to develop some heuristic models related to the Bethe approximation. The present work is an attempt to give a theoretical basis and firmer ground to these heuristics.

Preliminaries

Inverse Ising problem

In this section we consider binary variables (x i ∈ {0, 1}), which at our convenience may be also written as spin variables s i = 2x i -1 ∈ {-1, 1}. We assume that from a set of historical observations, the empirical mean mi (resp. covariance χij ) is given for each variable s i (resp. each pair of variable (s i , s j )). In this case, from Jayne's maximum entropy principle [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF], imposing these moments to the joint distribution leads to a model pertaining to the exponential family, the so-called Ising model:

P(s) = 1 Z[J, h] exp i h i s i + i,j J ij s i s j (2.1)
where the local fields h = {h i } and the coupling constants J = {J ij } are the Lagrange multipliers associated respectively to mean and covariance constraints when maximizing the entropy of P.

They are obtained as minimizers of the dual optimization problem, namely

(h ⋆ , J ⋆ ) = argmin (h,J) L[h, J], (2.2)
with

L[h, J] = log Z[h, J] - i h i mi - ij J ij mij (2.

3)

Inria the log likelihood. This leads to invert the linear response equations:

∂ log Z ∂h i [h, J] = mi (2.4) ∂ log Z ∂J ij [h, J] = mij , (2.5
) mij = mi mj + χij being the empirical expectation of s i s j . As noted e.g. in [START_REF] Cocco | Adaptive cluster expansion for the inverse Ising problem: convergence, algorithm and tests[END_REF], the solution is minimizing the cross entropy, a Kullback-Leibler distance between the empirical distribution P based on historical data and the Ising model:

D KL [ P P] = log Z[h, J] - i h i mi - i<j J ij mij -S( P).
(2.6)

The set of equations (2.4,2.5) cannot be solved exactly in general because the computational cost of Z is exponential. Approximations resorting to various mean field methods can be used to evaluate Z[h, J].

Plefka's expansion To simplify the problem, it is customary to make use of the Gibbs free energy, i.e. the Legendre transform of the free energy, to impose the individual expectations m = { mi } for each variable:

G[m, J] = h T (m)m + F [h(m), J] (with F [h, J] def = -log Z[h, J],
h T m is the ordinary scalar product) where h(m) depends implicitly on m through the set of constraints

∂F ∂h i = -m i . (2.7)
Note that by duality we have

∂G ∂m i = h i (m), (2.8)
and

∂ 2 G ∂m i ∂m j = ï dh dm ò ij = ï dm dh ò -1 ij = - ∂ 2 F ∂h i ∂h j -1 = χ -1 ij . (2.9)
i.e. the inverse susceptibility matrix. Finding a set of J ij satisfying this last relation along with (2.8) yields a solution to the inverse Ising problem since the m's and χ's are given. Still a way to connect the couplings directly with the covariance matrix is given by the relation

∂G ∂J ij = -m ij .
(2.10)

The Plefka expansion is used to expand the Gibbs free energy in power of the coupling J ij assumed to be small. Multiplying all coupling J ij by α yields the following cluster expansion:

G[m, αJ] = h T (m, α)m + F [h(m, α), αJ] (2.11) = G 0 [m] + ∞ n=0 α n n! G n [m, J] (2.12)
where each term G n corresponds to cluster contributions of size n in the number of links J ij involved, and h(m, α) depends implicitly on α in order to always fulfill (2.7). This precisely is the Plefka expansion, and each term of the expansion (2.12) can be obtained by successive derivation of (2.11). We have

G 0 [m] = i 1 + m i 2 log 1 + m i 2 + 1 -m i 2 log 1 -m i 2 .
Letting

H J def = i<j J ij s i s j , using (2.7), the two first derivatives of (2.11) w.r.t α read dG[m, αJ] dα = -E α H J , (2.13)
d 2 G[m, αJ] dα 2 = -Var α H J - i dh i (m, α) dα Cov α H J , s i , (2.14)
where subscript α indicates that expectations, variance and covariance are taken at given α. To get successive derivatives of h(m, α) one can use (2.8). Another possibility is to express the fact that m is fixed,

dm i dα = 0 = - d dα ∂F [h(α), αJ] ∂h i = j h ′ j (α) Cov α (s i , s j ) + Cov α (H J , s i ), giving h ′ i (α) = - j [χ -1 α ] ij Cov α (H J , s j ), (2.15)
where χ α is the susceptibility delivered by the model when α = 0. To get the first two terms in the Plefka expansion, we need to compute these quantities at α = 0:

Var H J = i<k,j J ij J jk m i m k (1 -m 2 j ) + i<j J 2 ij (1 -m 2 i m 2 j ), Cov H J , s i = j J ij m j (1 -m 2 i ), h ′ i (0) = - j J ij m j ,
(by convention J ii = 0 in these sums). The first and second orders then finally read:

G 1 [m, J] = - i<j J ij m i m j , G 2 [m, J] = - i<j J 2 ij (1 -m 2 i )(1 -m 2 j ),
and correspond respectively to the mean field and to the TAP approximation. Higher order terms have been computed in [START_REF] Georges | How to expand around mean-field theory using hightemperature expansions[END_REF].

Inria

At this point we are in position to find an approximate solution to the inverse Ising problem, either by inverting equation (2.9) or (2.10). To get a solution at a given order n in the coupling, solving (2.10) requires G at order n + 1, while it is needed at order n in (2.9).

Taking the expression of G up to second order gives

∂G ∂J ij = -m i m j -J ij (1 -m 2 i )(1 -m 2 j ),
and (2.10) leads directly to the basic mean-field solution:

J M F ij = χij (1 -m2 i )(1 -m2 j )
.

(2.16)

At this level of approximation for G, using (2.8) we also have

h i = 1 2 log 1 + m i 1 -m i - j J ij m j + j J 2 ij m i (1 -m 2 j )
which corresponds precisely to the TAP equations. Using now (2.9) gives

∂h i ∂m j = [χ -1 ] ij = δ ij 1 1 -m 2 i + k J 2 ik (1 -m 2 k ) -J ij -2J 2 ij m i m j .
Ignoring the diagonal terms, the TAP solution is conveniently expressed in terms of the inverse empirical susceptibility,

J T AP ij = - 2[ χ-1 ] ij 1 + 1 -8 mi mj [ χ-1 ] ij , (2.17)
where the branch corresponding to a vanishing coupling in the limit of small correlation i.e. small χij and [ χ-1 ] ij for i = j, has been chosen.

Bethe approximate solution

When the graph formed by the pairs (i, j) for which the correlations χij are given by some observations is a tree, the following form of the joint probability corresponding to the Bethe approximation:

P(x) = i<j pij (x i , x j ) p(x i )p(x j ) i pi (x i ), (2.18)
yields actually an exact solution to the inverse problem (2.2), where the p are the single and pair variables empirical marginal given by the observations. Using the following identity

log pij (s i , s j ) pi (s i )p j (s j ) = (1 + s i )(1 + s j ) 2 log p11 ij p1 i p1 j + (1 + s i )(1 -s j ) 2 log p10 ij p1 i p0 j + (1 -s i )(1 + s j ) 2 log p01 ij p0 i p1 j + (1 -s i )(1 -s j ) 2 log p00 ij p0 i p0 j
where the following parametrization of the p's px

i def = p 1 + s i 2 = x = 1 2 (1 + mi (2x -1)), (2.19
)

pxy ij def = p 1 + s i 2 = x, 1 + s j 2 = y = 1 4 (1 + mi (2x -1) + mj (2y -1) + mij (2x -1)(2y -1) (2.20)
relating the empirical frequency statistics to the empirical "magnetizations" m ≡ m, can be used. Summing up the different terms gives us the mapping onto an Ising model (2.1) with

h i = 1 -d i 2 log p1 i p0 i + 1 4 j∈∂i log p11 ij p10 ij p01 ij p00 ij , (2.21)
J ij = 1 4 log p11 ij p00 ij p01 ij p10 ij , ∀ (i, j) ∈ E, (2.22)
where d i is the number of neighbors of i, using the notation j ∈ ∂i for "j neighbor of i". The partition function is then explicitly given by

Z Bethe [p] = exp ï - 1 4
(i,j)∈E log p00 ij p01 ij p10 ij p11 ij - i 1 -d i 2 log(p 0 i p1 i ) ò (2.23)
The corresponding Gibbs free energy can then be written explicitly using (2.21-2.23). With fixed magnetizations m i 's, and given a set of couplings {J ij }, the parameters m ij are implicit function

m ij = m ij (m i , m j , J ij ),
obtained by inverting the relations (2.22). For the linear response, we get from (2.21) a result derived first in [START_REF] Welling | Approximate inference in Boltzmann machines[END_REF]:

∂h i ∂m j = 1 -d i 1 -m 2 i + 1 16 k∈∂i 1 p11 ik + 1 p01 ik 1 + ∂m ik ∂m i + 1 p00 ik + 1 p10 ik 1 - ∂m ik ∂m i δ ij + 1 16
1 p11 ij + 1 p10 ij 1 + ∂m ij ∂m i + 1 p00 ij + 1 p01 ij 1 - ∂m ij ∂m i δ j∈∂i .
Using (2.22), we can also express

∂m ij ∂m i = - 1 p11 ij + 1 p01 ij -1 p10 ij -1 p00 ij 1 p11 ij + 1 p01 ij + 1 p10 ij + 1 p00 ij ,
so that with little assistance of Maple, we may finally reach the expression given in [START_REF] Nguyen | Bethe-Peierls approximation and the inverse Ising model[END_REF] [

χ-1 ] ij = 1 -d i 1 -m 2 i + k∈∂i 1 -m 2 k (1 -m 2 i )(1 -m 2 k ) -χ2 ik δ ij - χij (1 -m 2 i )(1 -m 2 j ) -χ2 ij δ j∈∂i , (2.24)
equivalent to the original one derived in [START_REF] Welling | Approximate inference in Boltzmann machines[END_REF], albeit written in a different form, more suitable to discuss the inverse Ising problem. This expression is quite paradoxical since the inverse of the [χ] ij matrix, which coefficients appear on the right hand side of this equation, should coincide with the left hand side, given as input of the inverse Ising problem. The existence of an exact solution can therefore be checked directly as a self-consistency property of the input data χij : for a given pair (i, j) either: (2.24) has to hold and J ij is given by (2.22) using mij = mi mj + χij .

Inria • [ χ-1 ] ij = 0, then this self-consistency relation
• [ χ-1 ] ij = 0 then J ij = 0 but χij , which can be non-vanishing, is obtained by inverting [ χ-1 ] defined by (2.24).

Finally, complete consistency of the solution is checked on the diagonal elements in (2.24). If full consistency is not verified, this equation can nevertheless be used to find approximate solutions.

Remark that, if we restrict the set of equations (2.24), e.g. by some thresholding procedure, in such a way that the corresponding graph is a spanning tree, then, by construction, χ ij ≡ χij will be solution on this restricted set of edges, simply because the BP equations are exact on a tree. The various methods proposed for example in [28,[START_REF] Yasuda | Approximate learning algorithm in Boltzmann machines[END_REF] actually correspond to different heuristics for finding approximate solutions to this set of constraints. As noted in [START_REF] Nguyen | Bethe-Peierls approximation and the inverse Ising model[END_REF], a direct way to proceed is to eliminate χ ij in the equations obtained from (2.22) and (2.24):

χ 2 ij + 2χ ij (m i m j -coth(2J ij )) + (1 -m 2 i )(1 -m 2 j ) = 0 χ 2 ij - χ ij [χ -1 ] ij -(1 -m 2 i )(1 -m 2 j ) = 0.
This leads directly to

J Bethe ij = - 1 2 atanh 2[ χ-1 ] ij » 1 + 4(1 -m2 i )(1 -m2 j )[ χ-1 ] 2 ij -2 mi mj [ χ-1 ] ij , (2.25)
while the corresponding computed of χ ij , instead of the observed one χij , has to be inserted in (2.21) to be fully consistent. Note that J Bethe ij and J T AP ij coincide at second order in [ χ-1 ] ij .

Hopfield model As mentioned in the introduction when the distribution to be modeled is multi-modal, the situation corresponds to finding an Ising model in the low temperature phase with many modes, referred to as Mattis states in the physics literature. Previous methods assume implicitly a high temperature where only one single mode, "the paramagnetic state" is selected. The Hopfield model, introduced originally to model auto-associative memories, is a special case of an Ising model, where the coupling matrix is of low rank p ≤ N and corresponds to the sum of outer products of p given spin vectors {ξ k , k = 1 . . . p}, each one representing a specific attractive pattern:

J ij = 1 p N k=1 ξ k i ξ k j
In our inference context, these patterns are not given directly, the input of the model being the covariance matrix. In [START_REF] Cocco | High-dimensional inference with the generalized Hopfield model: Principal component analysis and corrections[END_REF] these couplings are interpreted as the contribution stemming from the p largest principle axes of the correlation matrix. This lead in particular the authors to propose an extension of the Hopfield model by introducing repulsive patterns to take as well into account the smallest principal axes. Assuming small patterns coefficients |ξ k | < 1/ √ N , they come up with the following couplings with highest likelihood:

J Hopf ield ij ≡ 1 » (1 -m 2 i )(1 -m 2 j ) p k=1 1 - 1 λ k v k i v k j - 1 λ N -k -1 v N -k i v N -k j
at first order of the perturbation. At this order of approximation the local fields are given by

h i = tanh(m i ) - j J Hopf ield ij m j .
In a previous study [START_REF] Furtlehner | Learning multiple belief propagation fixed points for real time inference[END_REF] we found a connection between the plain direct BP method with the Hopfield model, by considering a 1-parameter deformation of the Bethe inference model (2.18)

P(x) = i<j pij (x i , x j ) p(x i )p(x j ) α i pi (x i ), (2.26)
with α ∈ [0, 1]. We observed indeed that when the data corresponds to some multi-modal measure with well separated components, this measure coincides asymptotically with an Hopfield model made only of attractive pattern, representative of each component of the underlying measure. α represents basically the inverse temperature of the model and is easy to calibrate in practice.

2.2 More on the Bethe susceptibility The explicit relation (2.24) between susceptibility and inverse susceptibility coefficients is not the only one that can be obtained. In fact, it is the specific property of a singly connected factor graph that two variables x i and x j , conditionally to a variable x k are independent if k is on the path between i and j along the tree:

j k i i j k (c) (d) (e) (f) j i l k i j k l i j k l i j k l (a) (b)
p(x i , x j , x k ) = p(x i |x k )p(x j |x k )p(x k ) = p(x i , x k )p(x j , x k ) p(x k )
Using the parametrization (2.19,2.20) with m ij = m i m j + χ ij yields immediately

χ ij = χ ik χ jk 1 -m 2 k , ∀ k ∈ (i, j) along T . (2.27)
By recurrence we get, as noticed in e.g. [START_REF] Mora | Géométrie et inférence dans l'optimisation et en théorie de l'information[END_REF], given the path i 0 = i, i 1 , . . . , i n+1 = j between i and j along the tree T

χ ij = n a=0 χ iaia+1 n a=1 (1 -m 2 ia )
, reflecting the factorization of the joint measure. This expression actually coincides with (2.24) only on a tree. On a loopy graph, this last expression should be possibly replaced by a sum over paths.

Inria

Higher order susceptibility coefficients are built as well in terms of elementary building blocks given by the pairwise susceptibility coefficients χ ij . The notations generalize into the following straightforward manner:

m ijk def = E(s i s j s k ) def = m i m j m k + m i χ jk + m j χ ik + m k χ ij + χ ijk m ijkl def = E(s i s j s k s l ) def = m i m j m k m l + m i m j χ kl + m i m k χ jl + m i m l χ jk + m j m k χ il + m j m l χ ik + m k m l χ ij + m i χ jkl + m j χ ikl + m k χ ijl + m l χ ijk + χ ijkl ,
where χ ijk and χ ijkl are respectively three and four point susceptibilities. The quantities being related to the corresponding marginals similarly to (2.19,2.20):

p(s i , s j , s k ) = 1 8 1 + m i s i + m j s j + m k s k + m ij s i s j + m ik s i s k + m jk s j s k + m ijk s i s j s k p(s i , s j , s k , s l ) = 1 16 1 + m i s i + m j s j + m k s k + m l s l + m ij s i s j + m ik s i s k + m il s i s l + m jk s j s k + m jl s j s l + m kl s k s l + m ijk s i s j s k + m ijl s i s j s l + m ikl s i s k s l + m jkl s j s k s l + m ijkl s i s j s k s l
Using the basic fact that, on the tree

p(s i , s j , s k ) = p(s i , s j )p(s j , s k ) p(s i )
when j is on the path Ù ik given by T , and

p(s i , s j , s k ) = s l p(s i , s l )p(s j , s l )p(s k , s l ) p(s l ) 2
when path Ù ij, Ù ik and ı jk along T intersect on vertex l, we obtain

χ ijk =    -2 m l (1 -m 2 l ) 2 χ il χ jl χ kl with {l} = (i, j) ∩ (i, k) ∩ (j, k) along T , -2m j χ ik if j ∈ (i, k) along T .
For the fourth order, more cases have to be distinguished. When i, j, k and l are aligned as on Figure 2.1.c, in this order on the path Û il along T we have

p(s i , s j , s k , s l ) = p(s i , s j )p(s j , s k , s l ) p(s j ) 2
which leads to

χ ijkl = 4m k m j χ il -χ ik χ jl -χ il χ jk .
When path Ù ij, Ù ik and ı jk along T intersect on vertex l as on Figure 2.1.d, we obtain instead

1 χ ijkl = 2 3m 2 l -1 1 -m 2 l χ ij χ kl .
1 This apparently non-symmetric expression can be symmetrized with help of (2.27).

RR n°8059

For the situation corresponding to Figure 2.1.e, we have

p ( s i , s j , s k , s l ) = sq p(s i , s j , s q )p(s q , s k , s l ) p(s q ) 2
which leads to

χ ijkl = 2 3m 2 q -1 1 -m 2 q χ ij χ kl .
Finally, for the situation corresponding to Figure 2.1.f, we have

p ( s i , s j , s k , s l ) = sq p(s i , s j )p(s j , s k , s l ) p(s j ) 2 leading to χ ijkl = -χ ik χ jl -χ jk χ il + 4 m k m q 1 -m 2 q χ ij χ lq .

Sparse inverse estimation of covariance matrix

Let us leave the Ising modeling issue aside for a while and introduce another related graph selection problem, named sparse inverse covariance estimation, which is defined on real continuous random variables. This method aims at constructing a sparse factor graph structure by identifying conditionally independent pairs of nodes in the graph, given empirical covariances of random variables. Assuming that all nodes in the graph follow a joint multi-variate Gaussian distribution with mean µ and covariance matrix C, the existing correlation between the nodes i and j given all the other nodes in the graph are indicated by the non-zero ijth entry of the precision matrix C -1 , while zero entries correspond to independent pairs of variables. Therefore, under the joint normal distribution assumption, selection of factor graph structures amounts to finding the sparse precision matrix that best describes the underlying data distribution, given the fixed empirical covariance matrix. When the derived inverse estimation is sparse, it becomes easier to compute marginal distribution of each random variable and conduct statistical inference. To achieve that goal, optimizations methods have been developed based on L 0 or L 1 norm penalty for the estimation of C -1 , to enhance its sparse structure constraint on the estimated inverse of covariance matrix and discover underlying conditionally independent parts. Let Ĉ ∈ R n×n be the empirical covariance matrix of n random variables (represented as the nodes in the graph model). The sparsity penalized maximum likelihood estimation A of the precision matrix C -1 can be derived by solving the following positive definite cone program:

A = argmin X≻0 -L(X) + λP (X) (2.28) where L(A) def = log det(A) -Tr(A Ĉ), (2.29)
is the log likelihood of the distribution defined by A, log det being the logarithm of the determinant, and P (A) is a sparsity inducing regularization term [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. λ is the regularization coefficient balancing the data-fitting oriented likelihood and sparsity penalty. Since the precision matrix of joint normal distribution should be positive definite, any feasible solution to this optimization problem is thus required to locate within a positive definite cone. The penalty term P (A) is typically constructed using sparsity inducing matrix norm, also known as sparse learning in the domain of statistical learning. There are two typical configurations of P (A): • The L 0 norm A L0 of the matrix X, which counts the number of non-zero elements in the matrix. It is also known as the cardinality or the non-zero support of the matrix. Given its definition, it is easy to find that L 0 norm based penalty is a strong and intuitive appeal for sparsity structure of the estimated precision matrix. However, it is computationally infeasible to solve exact L 0 -norm minimization directly, due to the fact that exact L 0 norm penalty is discontinuous and non-differentiable. In practice, one either uses a continuous approximation to the form of the L 0 -penalty, or solve it using a greedy method. Due to the non-convexity of the exact L 0 norm penalty, only a local optimum of the feasible solution can be guaranteed. Nevertheless, L 0 norm penalty usually leads to much sparser structure in the estimation, while local optimum is good enough for most practical cases.

Inria
• The L 1 matrix norm A L1 = n i,j |A ij |. L 1 norm penalty was firstly introduced into the standard least square estimation problem by Tibshiharni [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], under the name "Lasso regression". Minimizing the L 1 norm based penalty encourages sparse non-zero entries in the estimated precision matrix A, which achieves a selection of informative variables for regression and reduces complexity of regression model efficiently. Further extension of the basic L 1 -norm penalty function allows one assigning different non-negative weight values

λ ij to different entries A ij , as n i,j λ ij |A ij |.
This weighted combination can constrain the sparse penalties only on the off-diagonal entries, so as to avoid unnecessary sparsity on the diagonal elements. Furthermore, this extension allows us to introduce prior knowledge about the conditional independence structure of the graph into the joint combination problem.

For further understanding of the relation between the exact L 0 and L 1 norm penalty, we illustrate them with respect to one scalar variable in Figure 2.2. As we can see, within [-1, 1], L 1 -norm penalty plays as a convex envelope of the exact L 0 -norm penalty. Due to the convexity property of L 1 norm penalty, the global optimum of the convex programming problem can be achieved with even linear computational complexity [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. However, although L 1 norm based penalty leads to computationally sound solutions to the original issue, it also introduces modeling bias into the penalized maximum likelihood estimation. As illustrated in the figure, when the underlying true values of the matrix entries are sufficiently large, the corresponding L 1 norm based regularization performs linearly increased penalty to those entries, which thus results in a severe bias w.r.t. the maximum likelihood estimation [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF]. In contrast, L 0 norm based penalty avoids such issue by constraining the penalties of non-zero elements to be constant. It has been proved in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] that the L 1 norm penalty discovers the underlined sparse structure when some suitable assumptions are satisfied. However, in general cases, the quality of the solutions is not clear.

3 Link selection at the Bethe point

Bethe approximation and Graph selection

As observed in the previous section, when using the Bethe approximation to find an approximate solution to the IIP, the consistency check should then be that either the factor graph be sparse, nearly a tree, either the coupling are small. There are then two distinct ways of using the Bethe approximation:

• the direct way, where the form of the joint distribution (2.18) is assumed with a complete graph. There is then a belief propagation fixed point for which the beliefs satisfy all the constraints. This solution to be meaningful requires small correlations, so that the belief propagation fixed point be stable and unique, allowing the corresponding log likelihood to be well approximated. Otherwise, this solution is not satisfactory, but a pruning procedure, which amounts to select a sub-graph based on mutual information, can be used. The first step is to find the MST with these weights. Adding new links to this baseline solution in a consistent way is the subject of the next section.

• the indirect way consists in first inverting the potentially non-sparse correlation matrix. If the underlying interaction matrix is actually a tree, this will be visible in the inverse correlation matrix, indicated directly by the non-zero entries. This seems to work better than the previous one when no sparsity but weak coupling is assumed. It corresponds in fact to the equations solved iteratively by the susceptibility propagation algorithm [28].

Let us first justify the intuitive assertion concerning the optimal model with tree like factor graphs, which was actually first proposed in [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF] and which is valid for any type of MRF.

Proposition 3.1. The optimal model with tree like factor graphs to the inverse pairwise MRF is obtained by finding the MST on the graph of weighted links with weights given by mutual information between variables.

Proof. On a tree, the Bethe distribution relating marginal p i of single and pair p ij variables distributions

P(x) = (i,j)∈T p ij (x i , x j ) p i (x i )p j (x j ) i p i (x i )
is exact. Assuming first that the associated factor graph is given by a tree T , the expression (2.6) of the log likelihood leads to the optimal solution given that tree:

p i = pi and p ij = pij , ∀i ∈ V, (i, j) ∈ T .
In turn we have the following expression for the log likelihood:

L = (i,j)∈T I ij - i H i , Inria with H i def = - xi pi (x i ) log pi (x i ) I ij def = xi,xj pij (x i , x j ) log pij (x i , x j ) pi (x i )p j (x j )
.

respectively the single variable entropy and the mutual information between two variables. Since the single variable contributions to the entropy is independent of the chosen graph T , the optimal choice for T correspond to the statement of the proposition.

Optimal 1-link correction to the Bethe point

Adding one link Suppose now that we want to add one link to the max spanning tree.

The question is which link will produce the maximum improvement to the log likelihood. This question is a special case of how to correct a given factor for a general pairwise model. Let P 0 be the reference distribution to which we want to add (or modify) one factor ψ ij to produce the distribution

P 1 (x) = P 0 (x) × ψ ij (x i , x j ) Z ψ (3.1)
with

Z ψ = dx i dx j p 0 ij (x i , x j )ψ ij (x i , x j ).
The log likelihood corresponding to this new distribution now reads

L 1 = L 0 + dx P(x) log ψ ij (x i , x j ) -log Z ψ .
The functional derivative w.r.t. ψ reads

∂L 1 ∂ψ ij (x i , x j ) = p(x i , x j ) ψ ij (x i , x j ) - p 0 ij (x i , x j ) Z ψ , ∀(x i , x j ) ∈ Ω 2 .
which yield immediately the solution

ψ ij (x i , x j ) = pij (x i , x j ) p 0 ij (x i , x j ) with Z ψ = 1, (3.2)
where p 0 (x i , x j ) is the reference marginal distribution obtained from P 0 . The correction to the Log likelihood simply reads

∆L = D KL (p ij p 0 ij ). (3.3)
Sorting all the links w.r.t. this quantity yields the (exact) optimal 1-link correction to be made. As a check, consider the special case of the Ising model. Adding one link amounts to set one J ij to some finite value, but since this will perturb at least the local magnetization of i and j, we have also to modify the local fields h i and h j . The modified measure therefore reads

P 1link = P Bethe exp J ij s i s j + h i s i + h j s j ∆Z(J, h) ,
where ∆Z(J, h) is multiplicative correction factor to the partition function. We have

∆Z(h i , h j , J ij ) = z 0 + z 1 mi + z 2 mj + z 3 m Bethe ij ,
after introducing the following quantities

z 0 def = e Jij +hi+hj + e -Jij -hi+hj + e -Jij +hi-hj + e Jij -hi-hj z 1 def = e Jij +hi+hj -e -Jij -hi+hj + e -Jij +hi-hj -e Jij -hi-hj z 2 def = e Jij +hi+hj + e -Jij -hi+hj -e -Jij +hi-hj -e Jij -hi-hj z 3 def = e Jij +hi+hj -e -Jij -hi+hj -e -Jij +hi-hj + e Jij -hi-hj .
The correction to the log likelihood is then given by

∆L(h i , h j , J ij ) = log ∆Z(h i , h j , J ij ) -h i mi -h j mj -J ij mij (3.4)
This is a concave function of h i , h j and J ij , and the (unique) maximum is obtained when the following constraints are satisfied:

∂ log(∆Z) ∂h i = z 1 + z 0 mi + z 3 mj + z 2 m Bethe ij ∆Z(h i , h j , J ij ) = mi , ∂ log(∆Z) ∂h j = z 2 + z 3 mi + z 0 mj + z 1 m Bethe ij ∆Z(h i , h j , J ij ) = mj , ∂ log(∆Z) ∂J ij = z 3 + z 2 mi + z 1 mj + z 0 m Bethe ij ∆Z(h i , h j , J ij ) = mij .
This constraints can be solved as follows. First let

U def = (1 + mi + mj + m Bethe ij ) e 2(hi+hj ) , V def = (1 -mi + mj -m Bethe ij ) e 2(hj -Jij ) , W def = (1 + mi -mj -m Bethe ij ) e 2(hi-Jij ) .
to obtain the linear system

U (1 -mi ) -V (1 + mi ) + W (1 -mi ) = (1 -mi -mj + m Bethe ij )(1 + mi ), U (1 -mj ) + V (1 -mj ) -W (1 + mj ) = (1 -mi -mj + m Bethe ij )(1 + mj ), U (1 -mij ) -V (1 + mij ) -W (1 + mij ) = -(1 -mi -mj + m Bethe ij )(1 -mij ).
Inverting this system yields the following solution

U = (1 -mi -mj + m Bethe ij ) (1 -mi -mj + mij ) (1 + mi + mj + mij ) V = (1 -mi -mj + m Bethe ij ) (1 -mi -mj + mij ) (1 -mi + mj -mij ) W = (1 -mi -mj + m Bethe ij ) (1 -mi -mj + mij ) (1 + mi -mj -mij ) Inria
Using the parametrization (2.20), we finally arrive at

e hisi+hj sj +Jij sisj = pij (s i , s j ) b ij (s i , s j ) Å b ij (1, 1)b ij (-1, 1)b ij (1, -1)b ij (-1, -1) pij (1, 1)p ij (-1, 1)p ij (1, -1)p ij (-1, -1) ã 1/4 .
where b ij (s i , s j ) is the joint marginal of s i and s j , obtained from the Bethe reference point. From these expressions, we can assess for any new potential link the increase (3.4) in the log likelihood.

After rearranging all terms, it takes indeed as announced the following simple form

∆L(h i , h j , J ij ) = si,sj pij log pij (s i , s j ) b ij (s i , s j ) = D KL (p b). (3.5)
The interpretation is therefore immediate: the best candidate is the one for which the Bethe solution yields the most distant joint marginal b ij to the targeted one pij given by the data. Note that the knowledge of the {b ij , (ij) / ∈ T } requires a sparse matrix inversion through equation (2.24), which renders the method a bit expensive in the Ising case. For Gaussian MRF, the situation is different, because in that case the correction to the log likelihood can be evaluated directly by another means. Indeed, the correction factor (3.2) reads in that case

ψ ij (x i , x j ) = exp - 1 2 (x i , x j ) [ Ĉij ] -1 -[C ij ] -1 (x i , x j ) T ,
where [ Ĉij ] and [C ij ] represent the restricted 2 × 2 covariance matrix corresponding to the pair (x i , x j ) of respectively the reference model and the current model specified by precision matrix A = C -1 . With a small abuse of notation the new model obtained after adding or changing link (ij) reads

A ′ = A + [ Ĉij ] -1 -[C ij ] -1 def = A + V.
with a log likelihoods variation given by:

∆L = C ii Ĉjj + C jj Ĉii -2C ij Ĉij C ii C jj -C 2 ij -2 -log Ĉii Ĉjj -Ĉ2 ij C ii C jj -C 2 ij .
Let us notice the following useful formula (see e.g. [START_REF] Chang | Inversion of a perturbed matrix[END_REF]):

(A + [V ij ]) -1 = A -1 -A -1 [V ij ](1 + A -1 [V ij ]) -1 A -1 = A -1 -A -1 [V ij ](1 + [C ij ][V ij ]) -1 A -1 , (3.6)
valid for a 2 × 2 perturbation matrix [V ij ]. Using this formula, the new covariance matrix reads

C ′ = A ′-1 = A -1 -A -1 [C ij ] -1 1 -[ Ĉij ][C ij ] -1 A -1 . (3.7)
Therefore the number of operations needed to maintain the covariance matrix after each add-on is O(N 2 ).

Let us now examine under which condition adding/modifying links in this way let the covariance matrix remain positive semi-definite. By adding a 2 × 2 matrix, we expect a quadratic correction to the determinant:

det(A ′ ) = det(A) det(1 + A -1 V ) = det(A) 1 + V ii C ii + V ij C ji + V ji C ij + V jj C jj + V ii V jj -V ij V ji C ii C jj -C ij C ji , = det A × det([C ij ]) det([ Ĉij ])
which is obtained directly because A -1 V has non zero entries only on column i and j. Multiplying V by some parameter α ≥ 0, define

P (α) def = det(1 + αA -1 V ) = α 2 det [C ij ][ Ĉij ] -1 - α -1 α .
When increasing α from 0 to 1, P (α) will vary from

1 to det([C ij ])/ det([ Ĉij ]) without canceling at some point iff [C ij ][ Ĉij ] -1 is definite positive. P (α) is proportional to the characteristic polynomial of a 2 × 2 matrix [C ij ]/[ Ĉij ] -1 of argument (α -1)/α, so A ′ remains positive definite if [C ij ][ Ĉij ] -1 is definite positive. Since the product of eigenvalues given by det([C ij ])/ det([ Ĉij ]
) is positive, one has to check for the sum, given by the trace of

[C ij ]/[ Ĉij ] -1 : C ii Ĉjj + C jj Ĉii -2C ij Ĉij > 0. (3.8)
Since [C ij ] and [ Ĉij ] are individually positive definite, we have

C ii C jj -C ij 2 > 0 and Ĉii Ĉjj -Ĉ2 ij > 0
from which we deduce that

C ii Ĉjj + C jj Ĉii 2 2 > C ii Ĉjj C jj Ĉii > C ij 2 Ĉ2 ij ,
giving finally that (3.8) is always fulfilled when both [C ij ] and [ Ĉij ] are non-degenerate. Each time a link is added to the graph, its number of loops increases by one unit, so in a sense (3.3) represent a 1-loop correction to the bare Bethe tree solution.

Removing one link To use this in an algorithm, it would also be desirable to be able to remove links as well, such that with help of a penalty coefficient per link, the model could be optimized with a desired connectivity level.

For the Gaussian model, if A is the coupling matrix, removing the link (i, j) amounts to chose a factor ψ ij in (3.1) of the form:

ψ ij (x i , x j ) = exp A ij x i x j
(x i and x j are assumed centered as in the preceding section). Again, let V denote the perturbation in the precision matrix such that A ′ = A + V is the new one. The corresponding change in the log likelihood then reads

∆L = log det(1 + A -1 V ) -Tr(V Ĉ).

Inria

Arranging this expression leads to

∆L = log 1 -2A ij C ij -A 2 ij det([C ij ]) + 2A ij Ĉij .
Using again (3.6) we get for the new covariance matrix

C ′ = C - A ij 1 -2A ij C ij -A 2 ij det([C ij ]) C ñ A ij C jj 1 -A ij C ij 1 -A ij C ij A ij C ii ô C, (3.9)
with again a slight abuse of notation, the 2 × 2 matrix being to be understood as a N × N matrix with non-zero entries corresponding to (i, i), (i, j), (j, i) and (j, j). To check for the positive-definiteness property of A ′ , let us observe first that

det(A ′ ) = det(A) × P (A ij ), with P (x) = (1 -x(C ij -C ii C jj ))(1 -x(C ij + C ii C jj )).
When x varies from 0 to A ij , P (x) should remains strictly positive to insure that A ′ is definite positive. This results in the following condition:

1 C ij -C ii C jj < A ij < 1 C ii C jj + C ij .
We are now equipped to define algorithms based on addition/deletion of links.

Imposing walk-summability

Having a sparse GMRF gives no guaranty to its compatibility with belief propagation. In order to be able to use the Gaussian Belief Propagation (GaBP) algorithm for performing inference tasks, stricter constraints have to be imposed. The more precise condition known for convergence and validity of the GaBP algorithm is called walk-summability (WS) and is extensively described in [START_REF] Malioutov | Walk-sums and belief propagation in Gaussian graphical models[END_REF]. The two necessary and sufficient conditions for WS that we consider here are:

(i) Diag(A) -|R(A)| is definite positive, with R(A) def = A -Diag(A) the off-diagonal terms of A. (ii) The spectral radius ρ(|R ′ (A)|) < 1, with R ′ (A) ij def = R(A)ij √ AiiAjj .
Adding one link Using the criterion developed in the previous section suppose that in order to increase the likelihood we wish to add a link (i, j) to the graph. The model is modified according to

A ′ = A + [ Ĉij ] -1 -[C ij ] -1 def = A + V.
We assume that A is WS and we want to express conditions under which A ′ is still WS. Using the definition (i) of walk summability and mimicking the reasoning leading to (3.8) we can derive a sufficient condition for WS by replacing A with W (A)

def = Diag(A) -|R(A)|. It yields det (W(A + αV )) = det (W(A)) det 1 + αW(A) -1 W(V ) = det (W ) î 1 + α ¶ W -1 ii V ii + W -1 jj V jj -2W -1 ij |V ij | © + α 2 Ä W -1 ii W -1 jj -(W -1 ji ) 2 ä V ii V jj -|V ij | 2 ó = det(W )Q(α),
since R(A) ij = 0 and shortening W (A) in W . A sufficient condition for WS of A ′ is that Q does not have any root in [0, 1]. Note that checking this sufficient condition imposes to keep track of the matrix (Diag(A) -|R|) -1 = W -1 def = IW and requires O(N 2 ) operations at each step, using (3.6). A more compact expression of Q is

Q(α) = 1 + α Tr IW ij W (V ) + α 2 det (W (V )) det IW ij ,
First let's tackle the special case where det W (V )IW ij = 0, the condition for WS of A ′ is then

Tr(IW ij W (V )) > -1.
Of course if the roots are not real, i.e. Tr(IW ij W (V )) 2 < 4 det W (V )IW ij , A ′ is WS. If none of these conditions is verified we have to check that both roots

-Tr(IW ij W (V ) ± Tr(IW ij W (V )) 2 -4 det (W (V )IW ij ) 2 det (W (V )IW ij ) , are not in [0, 1].
Modifying one link This is equivalent to adding one link in the sense that (3.1) and (3.2) are still valid. If we want to make use of (i) the only difference is that R(A) ij is not zero before adding V so W (A + αV ) = W (A) + αW (V ) does not hold in general. Instead we have

W (A + αV ) = W (A) + φ(αV, A), with φ(V, A) def = ï V ii -|V ij + A ij | + |A ij | -|V ji -A ji | + |A ij | V jj ò
So we can derive a condition for A ′ to be WS using, as for the link addition,

det(W (A + αV )) = det (W ) det 1 + W -1 φ(αV, A) = det (W ) Θ(α)
But now Θ(α), the equivalent of Q(α), is a degree 2 polynomial only by parts. More precisely, if

αV ij -A ij > 0 we have Θ(α) def = Q p (α) and else Θ(α) def = Q m (α)
with both Q p and Q m degree 2 polynomial. So by checking the sign of αV ij -A ij and the roots of Q p and Q m we have sufficient conditions for WS after modifying one link.

Another possible way for both adding or modifying one link is to estimate the spectral radius of |R ′ (A ′ )| through simple power iterations and concludes using (ii). Indeed if a matrix M as a unique eigenvalue of modulus equal to its spectral radius then the sequence

µ k def = b k , M b k b k , b k , with b k+1 def = M b k M b k
converges to this eigenvalue. While the model remains sparse, with connectivity K, a power iteration of R ′ (A ′ ) requires O(KN ) operations, and it is then possible to conclude about the WS of A ′ in O(KN ). Keeping track of W (A) -1 , which requires O(N 2 ) operations, at each step is not needed anymore but we have to test WS for each possible candidate link. Note that computing the spectral radius gives us a more precise information about the WS status of the model. In absence of penalty (ν = 0) the algorithm will simply generate a model for all different mean connectivities, hence delivering an almost continuous Pareto set of solutions, with all possible trade-off between sparsity and likelihood as long as walk summability is satisfied.

Greedy Graph construction Algorithms

Instead, with a fixed penalty the algorithm is converging toward a solution with a connectivity depending implicitly on ν; it corresponding roughly to the point K ⋆ where the slope is

∆L N ∆K (K ⋆ ) = ν.
If we want to use the backtracking mechanism allowed by the penalty term without converging to a specific connectivity, we may also let ν be adapted dynamically. A simple way is to adapt ν with the rate of information gain by letting

ν = η∆L add , with η ∈ [0, 1[,
where ∆L add corresponds to the gain of the last link addition. With such a setting ν is always maintained just below the information gain per link, allowing thus the algorithm to carry on toward higher connectivity. Of course this heuristic assumes a concave Pareto front.

4 Perturbation theory near the Bethe point

Linear response of the Bethe reference point

The approximate Boltzmann machines described in the introduction are obtained either by perturbation around the trivial point corresponding to a model of independent variables, the first order yielding the Mean-field solution and the second order the TAP one, either by using the linear response delivered in the Bethe approximation. We propose to combine in a way the two procedures, by computing the perturbation around the Bethe model associated to the MST with weights given by mutual information. We denote by T ⊂ E, the subset of links corresponding to the MST, considered as given as well as the susceptibility matrix [χ Bethe ] given explicitly by its inverse through (2.24), in term of the empirically observed ones χ. Following the same lines as the one given in Section 2, we consider again the Gibbs free energy to impose the individual expectations m = { mi } given for each variable. Let J Bethe = {K ij , (i, j) ∈ T } the set of Bethe-Ising couplings, i.e. the set of coupling attached to the MST s.t. corresponding susceptibilities are fulfilled and J = {J ij , (i, j) ∈ E} a set of Ising coupling corrections. The Gibbs free energy reads now

G[m, J] = h T (m)m + F h(m), J Bethe + J
where h(m) depends implicitly on m through the same set of constraints (2.7) as before. The only difference resides in the choice of the reference point. We start from the Bethe solution given by the set of coupling J Bethe instead of starting with an independent model. The Plefka expansion is used again to expand the Gibbs free energy in power of the coupling J ij assumed to be small. Following the same lines as in Section 2.1, but with G 0 now replaced by

G Bethe [m] = h T (m)m -log Z Bethe h(m), J Bethe ,
and h i , J Bethe and Z Bethe given respectively by (2.21,2.22,2.23) where E is now replaced by T , letting again

H 1 def = i<j J ij s i s j ,
and following the same steps (2.13,2.14,2.15) leads to the following modification of the local fields

h i = h Bethe i - j [χ -1 Bethe ] ij Cov Bethe (H 1 , s i ) ∀i ∈ V
to get the following Gibbs free energy at second order in α (after replacing H 1 by αH 1 ):

G[m, αJ] = G Bethe (m) -αE Bethe (H 1 ) - α 2 2 Var Bethe (H 1 ) - ij [χ -1 Bethe ] ij Cov Bethe (H 1 , s i ) Cov Bethe (H 1 , s j ) + o(α 2 ).
Inria This is the general expression for the linear response near the Bethe reference point that we now use.

G BLR [J] def = -E Bethe (H 1 ) (4.1) - 1 2 Var Bethe (H 1 ) - i,j [χ -1 Bethe ] ij Cov Bethe (H 1 , s i ) Cov Bethe (H 1 , s j ) . (4.2)
represents the Gibbs free energy at this order of approximation. It it is given explicitly through

E Bethe (H 1 ) = i<j J ij m ij Var Bethe (H 1 ) = i<j,k<l J ij J kl m ijkl -m ij m kl Cov Bethe (H 1 , s k ) = i<j J ij (m ijk -m ij m k )
where

m i def = E Bethe (s i ), m ij def = E Bethe (s i s j ) m ijk def = E Bethe (s i s j s k ), m ijkl def = E Bethe (s i s j s k s l )
are the moments delivered by the Bethe approximation. With the material given in Section 2.2 these are given in closed form in terms of the Bethe susceptibility coefficients χ Bethe . Concerning the log-likelihood, it is given now by:

L[J] = -G Bethe (m) -G BLR [J] - ij (J Bethe ij + J ij ) mij + o(J 2 ). (4.3)
G BLR is at most quadratic in the J's and contains the local projected Hessian of the log likelihood onto the magnetization constraints (2.7) with respect to this set of parameters. This is nothing else than the Fisher information matrix associated to these parameter J which is known to be positive-semidefinite, which means that the log-likelihood associated to this parameter space is convex. Therefore it makes sense to use the quadratic approximation (4.3) to find the optimal point.

Line search along the natural gradient in a reduced space

Finding the corresponding couplings still amounts to solve a linear problem of size N 2 in the number of variables which will hardly scale up for large system sizes. We have to resort to some simplifications which amounts to reduce the size of the problem, i.e. the number of independent couplings. To reduce the problem size we can take a reduced number of link into consideration, i.e. the one associated with a large mutual information or to partition them in a way which remains to decide, into a small number q of group G ν , ν = 1, . . . q. Then, to each group ν is associated a parameter α ν with a global perturbation of the form

H 1 = q ν=1 α ν H ν RR n°8059
where each H ν involves the link only present in G ν .

H ν def = (i,j)∈Gν J ij s i s j ,
and the J ij are fixed in some way to be discussed soon. The corresponding constraints, which ultimately insures a max log-likelihood in this reduced parameter space are then

∂G BLR ∂α ν = -Ê(H ν ).
This leads to the solution:

α µ = q ν=1 I -1 µν Ê(H ν ) -E Bethe (H ν )
where the Fisher information matrix I has been introduced and which reads in the present case

I µν = Cov Bethe (H µ , H ν ) - i =j (i,j)∈T [χ -1 Bethe ] ij Cov Bethe (H µ , s i ) Cov Bethe (H ν , s j ) (4.4)
The interpretation of this solution is to look in the direction of the natural gradient [START_REF] Amari | Natural gradient works efficiently in learning[END_REF][START_REF] Arnold | Information-Geometric Optimization Algorithms: A Unifying Picture via Invariance Principles[END_REF] of the log likelihood. The exact computation of the entries of the Fisher matrix involves up to 4th order moments and can be computed using results of Section 2.2. At this point, the way of choosing the groups of edges and the perturbation couplings J ij of the corresponding links, leads to various possible algorithms. For example, to connect this approach to the one proposed in Section 3.2, the first group of links can be given by the MST, with parameter α 0 and their actual couplings J ij = J Bethe ij at the Bethe approximation; making a short list of the q -1 best links candidates to be added to the graph, according to the information criteria 3.3, defines the other groups as singletons. It is then reasonable to attach them the value

J ij = 1 4 log p11 ij p 11 ij p00 ij p 00 ij p 01 ij p01 ij p 10 ij p10 ij ,
of the coupling according to (3.2), while the modification of the local fields as a consequence of (3.2) can be dropped since the Gibbs free energy take it already into account implicitly, in order to maintain single variable magnetization m i = mi correctly imposed.

Reference point at low temperature

Up to now we have considered the case where the reference model is supposed to be a tree and is represented by a single BP fixed point. From the point of view of the Ising model this corresponds to perturb a high temperature model in the paramagnetic phase. In practice the data encountered in applications are more likely to be generated by a multi-modal distribution and a low temperature model with many fixed points should be more relevant. In such a case we assume that most of the correlations are already captured by the definition of single beliefs fixed points and the residual correlations is contained in the co-beliefs of each fixed point. For a multi-modal distribution with q modes with weight w k , k = 1 . . . q and a pair of variables (s i , s j ) we indeed have

χ ij = q k=1 w k Cov(s i , s j |k) + q k=1 w k (E(s i |k) -E(s i ))(E(s j |k) -E(s j )) def = χ intra ij + χ inter ij , Inria
where the first term is the average intra cluster susceptibility while the second is the inter cluster susceptibility. All the preceding approach can then be followed by replacing the single Bethe susceptibility and higher order moments in equations (4.1,4.4) in the proper way by their multiple BP fixed point counterparts. For the susceptibility coefficients, the inter cluster susceptibility coefficients χ inter are given directly from the single variable belief fixed points. The intra cluster susceptibilities χ k are treated the same way as the former Bethe susceptibility. This means that the co-beliefs of fixed points k ∈ {1, . . . q} are entered in formula (2.24) which by inversion yields the χ k 's, these in turn leading to χ intra by superposition. Higher order moments are obtain by simple superposition. Improved models could be then searched along the direction indicated by this natural gradient.

L 0 norm penalized sparse inverse estimation algorithm

We propose here to use the Doubly Augmented Lagrange (DAL) method [START_REF] Iusem | Augmented Lagrangian methods and proximal point methods for convex optimization[END_REF][START_REF] Eckstein | Nonlinear proximal point algorithm using Bregman functions, with applications to convex programming[END_REF][START_REF] Dong | An efficient algorithm for l 0 minimization in wavelet frame based image restoration[END_REF] to solve the penalized log-determinant programming in (2.28). For a general problem defined as follows:

min x F (x) = f (x) + g(x) (5.1)
where f (x) and g(x) are both convex. DAL splits the combination of f (x) and g(x) by introducing a new auxiliary variable y. Thus, the original convex programming problem can be formulated as : min

x,y

F (x) = f (x) + g(y) s.t. x -y = 0 (5.2)
Then it advocates an augmented Lagrangian method to the extended cost function in (5.2). Given penalty parameters µ and γ, it minimizes the augmented Lagrangian function

L(x, y, ν, x, ỹ) = f (x) + g(y) + ν, x -y + µ 2 x -y 2 2 + γ 2 x -x 2 2 + γ 2 y -ỹ 2 2 (5.3)
where x and ỹ are the prior guesses of x and y that can obtained either from a proper initialization or the estimated result in the last round of iteration in an iterative update procedure. Since optimizing jointly with respect to x and y is usually difficult, DAL optimizes x and y alternatively. That gives the following iterative alternative update algorithm with some simple manipulations:

x k+1 = min x f (x) + µ 2 x -y k + νk 2 2 + γ 2 x -x k 2 2 y k+1 = min y g(y) + µ 2 x k+1 -y + νk 2 2 + γ 2 y -y k 2 2 νk+1 = νk + x k+1 -y k+1 (5.4)
where ν = 1 µ ν. As denoted in [START_REF] Dong | An efficient algorithm for l 0 minimization in wavelet frame based image restoration[END_REF] and [START_REF] Iusem | Augmented Lagrangian methods and proximal point methods for convex optimization[END_REF], DAL improves basic augmented Lagrangian optimization by performing additional smooth regularization on estimations of x and y in successive iteration steps. As a result, it guarantees not only the convergence of the scaled dual variable ν, but also that of the proximal variables x k and y k , which could be divergent in basic augmented Lagrangian method.

We return now to the penalized log-determinant programming in sparse inverse estimation problem, as seen in (2.28). The challenge of optimizing the cost function is twofold. Firstly, the exact L 0 -norm penalty is non-differentiable, making it difficult to find an analytic form of gradient for optimization. Furthermore, due to the log-determinant term in the cost function, it implicitly requires that any feasible solution to the sparse approximation A of the precision matrix should be strictly positive definite. The gradient of the log-determinant term is given by Ĉ -A -1 , which is not continuous in the positive definite domain and makes it impossible to obtain any second-order derivative information to speed up the gradient descent procedure. We hereafter use S ++ as the symmetric positive definite symmetric matrices that form the feasible solution set for this problem. By applying DAL to the cost function (2.28), we can derive the following formulation:

Ĵ(A, Z, Ã, Z, ν) = -log det(A) + Tr( ĈA) + λP (Z) + ν, A -Z + µ 2 A -Z 2 2 + γ 2 A -Ã 2 2 + γ 2 Z -Z 2 2 s.t. A, Z ∈ S ++ (5.5)
where Z is the auxiliary variable that has the same dimension as the sparse inverse estimation A. Ã and Z are the estimated values of A and Z derived in the last iteration step. The penalty parameter γ controls the regularity of A and Z. By optimizing A and Z alternatively, the DAL procedure can be easily formulated as an iterative process as follows, for some δ > 0:

A k+1 = argmin A -log det(A) + Tr( ĈA) + λP (Z k ) + ν k , A -Z k + µ 2 A -Z k 2 2 + γ 2 A -A k 2 2 Z k+1 = argmin Z λP (Z) + ν k , A k+1 -Z + µ 2 A k+1 -Z 2 2 + γ 2 Z -Z k 2 2 ν k+1 = ν k + δ(A k+1 -Z k+1 ) s.t. A k+1 , Z k+1 ∈ S ++ (5.6)
By introducing the auxiliary variable Z, the original penalized maximum likelihood problem is decomposed into two parts. The first one is composed mainly by the convex log-determinant programming term. Non-convex penalty is absorbed into the left part. Separating the likelihood function and the penalty leads to the simpler sub-problems of solving log-determinant programming using eigenvalue decomposition and L 0 norm penalized sparse learning alternatively. Each sub-problem contains only one single variable, making it applicable to call gradient descent operation to search local optimum. Taking ν = 1 µ ν, we can derive the following scaled version of DAL for the penalized log-determinant programming:

A k+1 = argmin A -log det(A) + Tr( ĈA) + µ 2 A -Z k + νk 2 2 + γ 2 A -A k 2 2 Z k+1 = argmin Z µ 2 A k+1 -Z + νk 2 2 + γ 2 Z -Z k 2 2 + λP (Z) νk+1 = νk + A k+1 -Z k+1 s.t. A k+1 , Z k+1 ∈ S ++ (5.7)
To attack the challenge caused by non-differentiability of the exact L 0 norm penalty, we make use of a differentiable approximation to L 0 -norm penalty in the cost function Ĵ, named as "seamless L 0 penalty" (SELO) in [START_REF] Dicker | Variable selection and estimation with the seamless-l 0 penalty[END_REF]. The basic definition of this penalty term is given as:

P SELO (Z) = i,j 1 log(2) log(1 + |Z i,j | |Z i,j | + τ ) (5.8)

Inria

where Z i,j denotes individual entry in the matrix Z and τ > 0 is a tuning parameter. As seen in Figure 2.2, as τ gets smaller, P (Z i,j ) approximates better the L 0 norm I(Z i,j = 0). SELO penalty is differentiable, thus we can calculate the gradient of P (Z) explicitly with respect to each Z i,j and make use of first-order optimality condition to search local optimum solution. Due to its continuous property, it is more stable than the exact L 0 norm penalty in optimization. As proved in [START_REF] Dicker | Variable selection and estimation with the seamless-l 0 penalty[END_REF], the SELO penalty has the oracle property with proper setting of τ . That's to say, the SELO penalty is asymptotically normal with the same asymptotic variance as the unbiased OLS estimator in terms of Least Square Estimation problem.

The first two steps in (5.7) are performed with the positive definite constrains imposed on A and Z. The minimizing with respect to A is accomplished easily by performing Singular Vector Decomposition (SVD). By calculating the gradient of Ĵ with respect to A in (5.7), based on the first-order optimality, we derive:

Ĉ -A -1 + µ(A -Z k + νk ) + γ(A -A k ) = 0
(5.9)

Based on generalized eigenvalue decomposition, it is easy to verify that A k+1 = V Diag(β)V T , where V and {d i } are the eigenvectors and eigenvalues of µ(Z k -νk ) -Ĉ + γA k . β i is defined as:

β i = d i + » d i 2 + 4(ν + γ) 2(ν + γ) (5.10)
Imposing Z ∈ S ++ directly in minimizing the cost function with respect to Y make the optimization difficult to solve. Thus, instead, we can derive a feasible solution to Z by a continuous search on µ. Based on spectral decomposition, it is clear that X k+1 is guaranteed to be positive definite, while it is not necessarily sparse. In contrast, Z is regularized to be sparse while not guaranteed to be positive definite. µ is the regularization parameter controlling the margin between the estimated X k+1 and the sparse Z k+1 . Increasingly larger µ during iterations makes the sequences {X k } and {Z k } converge to the same point gradually by reducing margin between them. Thus, with enough iteration steps, the derived Z k follows the positive definite constraint and sparsity constraint at the same time. We choose here to increase µ geometrically with a positive factor η > 1 after every N µ iterations until its value achieves a predefined upper bound µ max . With this idea, the iterative DAL solution to the L 0 norm penalty is given as:

A k+1 = argmin A -log det(A) + Tr( ĈA) + µ 2 A -Z k + νk 2 2 + γ 2 A -A k 2 2 , Z k+1 = argmin Z µ 2 A k+1 -Z + νk 2 2 + γ 2 Z -Z k 2 2 + λP (Z) νk+1 = νk + A k+1 -Z k+1
µ k+1 = min µ η ⌊k/Nµ⌋ , µ max .

(5.11)

In the second step of (5.11), we calculate the gradient of the cost function with respect to Z and achieve the local minimum by performing the first-order optimum condition on it. Therefore, the updated value of each entry of Z is given by a root of a cubic equation, as defined below: if Z i,j > 0, Z i,j is the positive root of

2Z i,j 3 + (3τ -2θ i,j )Z i,j 2 + (τ 2 -3τ θ i,j )Z i,j -τ 2 θ i,j + λτ µ + γ = 0 if Z i,j < 0, Z i,j is the negative root of 2Z i,j 3 -(3τ + 2θ i,j )Z i,j 2 + (τ 2 + 3τ θ i,j )Z i,j -τ 2 θ i,j - λτ µ + γ = 0 else Z i,j = 0 (5.12)
where Z i,j is one single entry of Z and θ i,j = γZ k i,j + µ(A k+1 i,j + νk ) µ + γ .

Solving the cubic equations can be done rapidly using Cardano's formula within a time cost O(n 2 ). Besides, the spectral decomposition procedure has the general time cost O(n 3 ). Given the total number of iterations K, theoretical computation complexity of DAL is O(Kn 3 ). For our experiments, we initialize µ to 0.06, the multiplier factor η to 1.3 and the regularization penalty parameter γ to 10 -4 . To approximate the L 0 norm penalty, τ is set to be 5 • 10 -4 . In our experiment, to derive the Pareto curve of the optimization result, we traverse different values of λ. Most learning procedures converge with no more than K = 500 iteration steps.

To validate performance of sparse inverse estimation based on the L 0 norm penalty, we involve an alternative sparse inverse matrix learning method using L 1 norm penalization for comparison. Taking P (A) in (2.28) to be the L 1 matrix norm of A, we strengthen conditional dependence structure between random variables by jointly minimizing the negative log likelihood function and the L 1 norm penalty of the inverse matrix. Since L 1 norm penalty is strictly convex, we can use a quadratic approximation to the cost function to search for the global optimum, which avoids singular vector decomposition with complexity of O(p 3 ) and improves the computational efficiency of this solution to O(p), where p is the number of random variables in the GMRF model. This quadratic approximation based sparse inverse matrix learning is given in [START_REF] Hsieh | Sparse inverse covariance matrix estimation using quadratic approximation[END_REF], named as QUIC. We perform it directly on the empirical covariance matrix with different settings of the regularization coefficient λ. According to works in compressed sensing, the equality between L 1 norm penalty and L 0 norm penalty holds if and only if the design matrix satisfies restricted isometry property. However, restricted isometry property is sometimes too strong in practical case. Furthermore, to our best knowledge, there is no similar necessary condition guaranteeing equivalence between L 1 and L 0 norm penalty in sparse inverse estimation problem. Therefore, in our case, L 1 norm penalized log-determinant programming is highly likely to be biased from the underlying sparse correlation structure in the graph, which leads to much denser inverse matrices.

Experiments

In this section, various solutions based on the different methods exposed before are compared. We look first at the intrinsic quality, given either by the exact log likelihood for the Gaussian case, or by the empirical one for the Ising model, and then at its compatibility with belief propagation for inference tasks.

Inverse Ising problem Let us start with the inverse Ising problem. The first set of experiments illustrates how the linear-response approach exposed in Section 2 works when the underlying model to be found is itself an Ising model. The quality of the solution can then be assessed directly by comparing the couplings J ij found with the actual ones. Figure 6.1 are obtained by generating at random 10 3 Ising models of small size N = 10 either with no local fields (h i = 0, ∀i = 1 . . . N ) or with random centered ones h i = U [0, 1] -1/2 and with couplings

J ij = J √ N/3
(2 * U [0, 1] -1), centered with variance J 2 /N , J being the common rescaling factor corresponding to the inverse temperature. A glassy transition is expected at J = 1. The couplings are then determined using (2.16), (2.17), (2.22) and (2.25) respectively for the mean-field, TAP, BP and Bethe (equivalent to susceptibility propagation) solutions. the Bethe approximation yields the most precise results in absence of local fields while it is equivalent to TAP when a local field is present as shown on Figure 6.1.b. Since we want to use these methods in conjunction with BP we have also compared the BP-susceptibilities they deliver. To do that, we simply run BP to get a set of belief and co-beliefs in conjunction with equation (2.24) which after inversion yields a susceptibility matrix to be compared with the exact ones. The comparison shown on Figure 6.1.c indicates that Bethe and TAP yield the best results in absence of local field which are less robust when compared to the more naive BP method when local fields are present as seen on Figure 6.1.d. This is due to the fact that BP delivers exact beliefs when model (2.22,2.21) is used, which is not necessarily the case for other methods when the local fields are non-vanishing. It is actually not a problem of accuracy but of BP compatibility which is raised by this plot.

Sparse inverse models

Let us now test the approach proposed in Section 3.2 to build a model link by link for comparison with more conventional optimization schema based on L 0 and L 1 penalizations. We show the results of tests only for the Gaussian case where the link surgery corresponds to the IAU data of size N = 1000. < K > is the ratio of links to nodes N . The log likelihood on the y axes is unormalized and corresponds to (2.29). The reference is the log likelihood given by the full inverse covariance matrix.

can be treated exactly. In the Ising case the method can be used only marginally to propose little correction on the maximum spanning tree or any other sparse model. In general we cannot expect the method to be able to compete with the inverse susceptibility propagation schema i.e. what we call here the Bethe inverse model (2.25). The reason is that the LL gain given by one link is more costly to assess than in the Gaussian case and it is also only approximate. So the stability of the schema is more difficult to control when many links have to be added because the condition of the validity of the Bethe approximation are not controlled without paying an additional computational price. For the Gaussian case instead the situation is much more favorable because the gain can be computed exactly with low computational cost even when the graph is dense. The test we show on Figure 6.2 are done on simulated data produced by the traffic simulator METROPOLIS [START_REF] De Palma | Real cases applications of the fully dynamic METROPOLIS tool-box: an advocacy for large-scale mesoscopic transportation systems[END_REF], our original motivation for this work being related to traffic inference [START_REF] Furtlehner | A belief propagation approach to traffic prediction using probe vehicles[END_REF][START_REF] Furtlehner | Spatial and temporal analysis of traffic states on large scale networks[END_REF]. The first set corresponds to a small traffic network called Sioux-Fall consisting of 72 links, from which we extract the N = 60 most varying ones (the other one being mostly idle). The second set (IAU) is obtained for a large scale albeit simplified network of the Paris agglomeration of size 13626 links, out of which we extracted a selection of the N = 1000 most varying ones. Each sample data is a N -dimensional vector of observed travel times { ti , i = 1 . . . N }, giving a snapshot of the network at a given time in the day. The total number of samples is S = 3600 for Sioux-Falls and S = 7152 for IAU, obtained by generating many days of various traffic scenarios. Then for each link the travel time distribution is far from being Gaussian, having heavy tails in particular. So to deal with normal variables (if taken individually) we make the following standard transformation:

y i = F -1 Gauss Fi (t i ), ∀i = 1 . . . N (6.1)
which map the travel time t i to a genuine Gaussian variable y i , where Fi and F Gauss are respectively the empirical cdf of t i and of a centered normal variable. The input of the different algorithms under study is then the covariance matrix Cov(y i , y j ). This mapping will actually be important in the next section when using the devised MRF for inference tasks. Figure 6.2 displays the comparison between various methods. Performances of the greedy method are com-Inria parable to the L 0 penalized optimization. To generate one solution both methods are comparable also in term of computational cost, but the greedy is faster in very sparse regime, and since it is incremental, it generate a full Pareto subset for the cost of one solution. On this figure we see also that the L 1 method is simply not adapted to this problem. From the figures, we can see that the estimated inverse matrix derived based on L 1 norm penalty needs distinctively more non-zero entries to achieve similar log-likelihood level as the L 0 penalty, indicating its failure of discovering the underlying sparse structure, the thresholding of small non-zero entries being harmful w.r.t. positive definitness. The reason might be that is adapted to situations where a genuine sparse structure exists, which is not the case in the present data. Inverse models for inference We turn now to experiments related to the original motivation of this work, which is to use calibrated model for some inference tasks. The experiments goes as follows: we have an historical data set consisting of a certain number of samples, each one being a N -dimensional variable vector, say travel times, which serves to build the models2 . Given a sample test data we want to infer the (1 -ρ)N hidden variables when a certain fraction ρ of the variables are revealed. In practice we proceed gradually on each test sample by revealing one by one the variables in a random order and plot as a function of ρ the L 1 error made by the inference model on the hidden variables. Both for Ising and Gaussian MRF, the inference is not performed in the original variable space, but in an associated one obtained through a mapping (a traffic index) using the empirical cumulative distribution of each variable. For the Gaussian model the inference is performed in the index space defined previously by (6.1). For the Ising models we have studied a variety of possible mapping [START_REF] Martin | Encoding dependencies between real-valued observables with a binary latent MRF[END_REF] in order to associate a binary variable to a real one such that a belief associated to a binary state can be converted back into a travel time prediction. Without entering into the details (see [START_REF] Martin | Encoding dependencies between real-valued observables with a binary latent MRF[END_REF] for details), to define in practice this binary state σ i , either we make use of the median value x 1/2 i = F -1 i (1/2) in the distribution of x i for all i = 1 . . . N :

σ i = 1 1 {xi>x 1/2 i } (i).

Furtlehner & Han & others

Either we perform a soft mapping using the cdf:

P (σ i = 1) = Fi (x i ) (ii),
the last one having the advantage of being functionally invertible if F -1 i is defined, while the former one being inverted using Bayes rule. The data we are considering are "low temperature" data in the sense that correlations are too strong for an Ising model with one single fixed point. This is reflected in the fact that none of the basic methods given in the Section 2 is working. To overcome this we use a simple heuristic which consists in to add a parameter α ∈ [0, 1] in the BP model like e.g. in (2.26) or to multiply the J ij by α for the MF, TAP and Bethe models, the local field being consequently modified owing to their dependency on the J ij . Concerning the factor-graph we have considered various graph selection procedures. All are based on the mutual information given empirically between variables. A global/local threshold can be used to construct the graph, the parameter being the mean/local connectivity K; the MST can be used conveniently as a backbone and additional links are obtained through the thresholding selection procedures. These two parameter α and K are calibrated such as to optimize the performance for each type of model so that fair comparisons can be made afterward. One important difference between the Ising model and the Gaussian one is that multiple fixed points may show up in the Ising case while only a single one, stable or not stable, is present in the Gaussian case. This can be an advantage in favor of the Ising model when the data have well separated clusters. Ideally we would like a perfect calibration of the Ising model in order that these fixed points be located at the center of each cluster. The method proposed in Section 4 could help to do this, but has not been implemented yet. On Figure 6.3.a we see also that the knn predictor performs optimally in this case, since the error curve coincides exactly with the one given by the hidden generative Inria model of the data (ground truth). Figure 6.4 shows how the different models compare on the data generated by the traffic simulator. On the Sioux-Falls network, the Gaussian model gives the best results, and a sparse version obtained with the greedy algorithm of section 3.4 reach the same level of performance and outperforms knn. The best Ising model is obtained with the (2.26) with type (ii) encoding. For IAU the full Gaussian model is also competitive w.r.t knn, but the best sparse GaBP model is not quite able to follow. In fact the correlations are quite high in this data, which explain why the best Ising model shows very poor performance. The best Ising model in that case corresponds to the plain BP model with type (ii) encoding and MST graph.

Conclusion

This paper is based on the observation that the Bethe approximation can be in many case a good starting point for building inverse models from data observations. We have developed different ways of perturbing such a mean-field solution valid both for binary and Gaussian variables, and leading to an efficient algorithm in the Gaussian case to generated sparse approximation models compatible with BP. The additional requirement that the model be compatible with BP for large scale application discards dense models and simplifies in a way the search space on model selection. More experimental tests on various data should help to refine and settle the general methods proposed here.

Figure 2 . 1 :

 21 Figure 2.1: Various cumulant topologies of order three (a,b) and four (c-f).

Figure 2 . 2 :

 22 Figure 2.2: Demonstration of L 1 norm, exact L 0 norm and Seamless L 0 norm based penalty function.

Inria

  Removing one linkRemoving one link of the graph will change the matrix A in A ′ such as |R ′ (A ′ )| ≤ |R ′ (A)| where the comparison (≤) between two matrices must be understood as the element-wise comparison. Then dealing with positive matrices elementary results gives us ρ(|R ′ (A ′ )|) ≤ ρ(|R ′ (A)|) and thus removing one link of a WS model provide a new WS model.

Algorithm 1 :Algorithm 2 :

 12 Incremental graph construction by link addition S1 INPUT: the MST graph, and corresponding covariance matrix C. S2 : select the link with highest ∆L compatible with the WS preserving condition of A ′ in the Gaussian case. Update C according to (3.7) for the Gaussian model. For the Ising model, C is updated by first running BP to generate the set of beliefs and co-beliefs supported by the current factor graph, which in turn allows one to use (2.24) to get all the missing entries in C by inverting χ -1 . S3 : repeat S2 until convergence (i) or until a target connectivity is reached (ii) S4 : if (ii) repeat S2 until convergence by restricting the link selection in the set of existing ones. The complexity is respectively O(N 3 ) both for the Gaussian and the Ising model in the sparse domain where O(N ) links are added and respectively O(N 4 ) and O(N 5 ) in the dense domain. Indeed, the inversion of χ -1 in the Ising case costs O(N 2 ) operations as long as the factor graph remains sparse, but O(N 3 ) for dense matrices. Graph surgery by link addition/deletion S1 INPUT: the MST graph, and corresponding covariance matrix C, a link penalty coefficient ν. S2 : select the modification with highest ∆L-sν, with s = +1 for an addition and s = -1 for a deletion, compatible with the WS preserving condition of A ′ in the Gaussian case. Update C according to (3.7) and (3.9) respectively for an addition or a change and a deletion for the Gaussian model. For the Ising model, C is updated by first running BP to generate the set of beliefs and co-beliefs supported by the current factor graph, which in turn allows one to use (2.24) to get all the missing entries in C by inverting χ -1 . S3 : repeat S2 until convergence.

Figure 6 . 1 .Figure 6 . 1 :

 6161 Figure 6.1: Comparison between various approximate solutions to the inverse Ising problem. RMSE errors as a function of the temperature are plotted in (a) and (b) for the couplings J ij , in (c) and (d) for the susceptibility matrix χ ij obtained from the corresponding BP fixed point. Local fields h i are zero in (a) and (c) and finite but zero in average in (b) and (d).

Figure 6 . 2 :

 62 Figure 6.2: Comparison of the greedy information based MRF inference method with L 0 and L 1 norm penalized optimizations. (a) corresponds to the Sioux-Fall data of size N = 60. (b)corresponds to the IAU data of size N = 1000. < K > is the ratio of links to nodes N . The log likelihood on the y axes is unormalized and corresponds to(2.29). The reference is the log likelihood given by the full inverse covariance matrix.

Figure 6 . 3 :

 63 Figure 6.3: Comparison of decimation curves in the case of a multi-modal distribution with five cluster for N = 500 variables (a). Projection of the dataset in the 3-d dominant PCA space along with the corresponding BP fixed points projections obtained with the Ising model.

Figure 6 . 4 :

 64 Figure 6.4: Comparison of decimation curves between various MRF for Sioux-Falls data (a) and IAU data (b)

Figure 6 .

 6 3 illustrates this point. The data are sampled from a distribution containing 5 modes, each one being a product form over N random bimodal distributions attached to each link. On Figure6.3.a which displays the error as a function of the fraction of revealed variables we see that the Ising model obtained with(2.26), encoded with the median value (i), gives better prediction than the exact Gaussian model or the approximated GaBP compatible one. Indeed Figure6.3.b shows a projection of the data in the most relevant 3-d PCA space along with the projected position of BP fixed points (given by their sets of beliefs) delivered by the Ising model. As we see, the model is able to attach one BP fixed point to each component of the distribution.

RR n°8059

In fact the pairwise MRF models exploit only pairwise observations but for sake of comparison with a knn predictor we generate complete historical samples data.RR n°8059

Acknowledgments This work was supported by the French National Research Agency (ANR) grant N°ANR-08-SYSC-017.