Cyril Furtlehner

An Ising Model for Road Traffic Inference

We review some properties of the "belief propagation" algorithm, a distributed iterative map, used to perform Bayesian inference and present some recent work where this algorithm serves as a starting point to encode observation data into a probabilistic model and to process large scale information in real time. A natural approach is based on the linear response theory and various recent instantiations are presented. We will focus on the particular situation where the data have many different statistical components, representing a variety of independent patterns. As an application, the problem of reconstructing and predicting traffic states based on floating car data is then discussed.

Introduction

The "belief propagation" algorithm BP, originated in the artificial intelligence community for inference problems on Bayesian networks [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Network of Plausible Inference[END_REF]. It is a nonlinear iterative map which propagates information on a dependency graph of variables in the form of messages between variables. It has been recognised to be a generic procedure, instantiated in various domains like error correcting codes, signal processing or constraints satisfaction problems with various names depending on the context [START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF]: the forward-backward algorithm for hidden Markov model selection; the Viterbi algorithm; Gallager's sum-product algorithm in information theory. It has also a nice statistical physics interpretation in the context of mean-field theories, as a minimizer of a Bethe free energy [START_REF] Yedidia | Generalized belief propagation[END_REF], a solver of the cavity equations [START_REF] Mézard | Spin Glass Theory and Beyond[END_REF] and its relation to the TAP equations in the spin-glass context [START_REF] Kabashima | Belief propagation vs. tap for decoding corrupted messages[END_REF]. As a noticeable development in the recent years, related to the connection with statistical physics, is the emergence of a new generation of algorithms for solving difficult combinatorial problems, like the survey propagation algorithm [START_REF] Mézard | The random K-satisfiability problem: from an analytic solution to an efficient algorithm[END_REF] for constraint satisfaction problems or the affinity propagation for clustering [START_REF] Frey | Clustering by passing messages between data points[END_REF].

The subject of the present review is at first a statistical modelling problem. Assuming a set of high dimensional data, in the form of sparse observations covering a finite fraction of segments in a traffic network, we wish to encode the dependencies between the variables in a probabilistic model, which turns out to be a Markov random field (MRF). We proceed in such a way as to insure that inference on this MRF with BP is optimal in some way so that it can be fast and precise at the same time, offering the possibility to address large scale problems like inferring congestion on a macroscopic traffic network. In Section 1.2 we introduce the BP algorithm and review some of its properties. Section 1.3 is devoted to the general problem of encoding observation data, by addressing the inverse Ising problem. In Section 1.4 a traffic application is described along with the construction of an inference model. Section 1.5 is concerned with the problem of multiplicity of BP fixed points and how to turn this into an advantage when the underlying empirical distribution based on observational data is multi-modal. Finally in Section 1.6 we present some preliminary tests of the method.

The Belief Propagation Algorithm

We consider a set of discrete random variables x = {x i , i ∈ V} ∈ {1, . . . , q} |V| obeying a joint probability distribution of the form

P(x) = a∈F ψ a (x a ) i∈V φ i (x i ), (1.1)
where φ i and ψ a are factors associated respectively to a single variable x i and to a subset a ∈ F of variables, F representing a set of cliques and x a def = {x i , i ∈ a}. The ψ a are called the "factors" while the φ i are there by convenience and could be reabsorbed in the definition of the factors. This distribution can be conveniently represented with a bi-bipartite graph, called the factor graph [START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF]; F together with V define the factor graph G, which will be assumed to be connected. The set E of edges contains all the couples (a, i) ∈ F ×V such that i ∈ a. We denote d a (resp. d i ) the degree of the factor node a (resp. to the variable node i). The factor graph on the tributions, called the belief b(x i = x) associated to each variable i can be done efficiently. The BP algorithm does this effectively for all variables in one single procedure, by remarking that the computation of each of these marginals involves intermediates quantities called the messages m a→i (x i ) [resp. n i→a (x i )] "sent" by factor node a to variable node i [resp. variable node i to factor node a], and which are necessary to compute other marginals. The idea of BP is to compute at once all these messages, using the relation among them as a fixed point equation. Iterating the following message update rules sketched on Figure 1.1.b:

       m a→i (x i ) ← xa\xi j∈a\i n j→a (x j )ψ a (x a ), n i→a (x i ) ← φ i (x i ) b∋i m b→i (x i ),
yields, when a fixed point is reached, the following result for the beliefs,

b(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 ).
This turns out to be exact if the factor graph is a tree, but only approximate on multiply connected factor graphs. As mentioned before, this set of beliefs corresponds to a stationary point of a variational problem [START_REF] Yedidia | Generalized belief propagation[END_REF]. Indeed, consider the Kullback-Leibler divergence between a test joint distribution b(x) and the reference p(x). The Bethe approximation leads to the following functional of the beliefs, including the joint beliefs b a (x a ) corresponding to each factor:

D KL (b p) = {x} b({x}) log b({x}) p({x}) ≈ a,xa b a (x a ) log b a (x a ) ψ(x a ) i∈a b i (x i ) + i,xi log b i (x i ) φ i (x i ) def = F Bethe = E -S Bethe .
This is equivalent to say that we look for a minimizer of D KL (b p) in the following class of joint probabilities:

b(x) = a b a (x a ) i∈a b i (x i ) i b i (x i ), (1.2)
under the constraint that

xa\xi b a (x a ) = b i (x i ) ∀a ∈ F, ∀i ∈ a,
and that

x\xa b(x) ≈ b a (x a ), ∀a ∈ F, (1.3)
is valid, at least approximately. For a multi-connected factor graph, the beliefs b i and b a are then interpreted as pseudo-marginal distribution. It is only when G is simply connected that these are genuine marginal probabilities of the reference distribution p.

There are a few properties of BP that are worth mentioning at this point. Firstly, BP is a fast converging algorithm:

• Two sweeps over all edges are needed if the factor-graph is a tree.

• The complexity scales heuristically like KN log(N ) on a sparse factorgraph with connectivity K ≪ N . • It is N 2 for a complete graph. However, when the graph is multiply connected, there is little guarantee on the convergence [START_REF] Mooij | On the properties of the Bethe approximation and loopy belief propagation on binary network[END_REF] even so in practice it works well for sufficiently sparse graphs. Another limit in this case, is that the fixed point may not correspond to a true measure, simply because (1.2) is not normalized and (1.3) is approximate. In this sense, the obtained beliefs, albeit compatible with each other are considered only as pseudo-marginals. Finally, for such graphs, the uniqueness of fixed points is not guaranteed, but it has been shown that:

• stable BP fixed points are local minima of the Bethe free energy [START_REF] Heskes | On the uniqueness of loopy belief propagation fixed points[END_REF];

• the converse is not necessarily true [START_REF] Watanabe | Graph zeta function in the bethe free energy and loopy belief propagation[END_REF].

There are two important special cases, where the BP equations simplify: (i) For binary variables: x i ∈ {0, 1}. Upon normalization, the messages are parametrised as:

m a→i (x i ) = m a→i x i + (1 -m a→i )(1 -x i ),
which is stable w.r.t. the message update rule. The propagation of information reduces then to the scalar quantity m a→i . (ii) For Gaussian variables, the factors are necessarily pairwise, of the form

ψ ij (x i , x j ) = exp -A ij x i x j , φ i (x i ) = exp - 1 2 A ii x 2 i + h i x i .
Since factors are pairwise, messages can be seen as sent directly from one variable node i to another j with a Gaussian form:

m i→j (x j ) = exp - (x j -µ i→j ) 2 2σ i→j .
This expression is also stable w.r.t. the message update rules. Information is then propagated via the 2-component real vector (µ i→j , σ i→j ) with the following update rules:

µ i→j ←- 1 A ij h i + k∈∂i\j µ k→i σ k→i , σ i→j ←-- 1 A 2 ij A ii + k∈∂i\j σ -1 k→i .
At convergence the belief takes the form:

b i (x) = σ i 2π exp - (x -µ i ) 2 2σ i with µ i = σ i h i + j∈∂i µ j→i σ j→i σ -1 i = A ii + j∈∂i σ -1
j→i and the estimated covariance between x i and x j reads

σ ij = 1 A ij (1 -A 2 ij σ i→j σ j→i )
.

In this case, there is only one fixed point even on a loopy graph, not necessarily stable, but if convergence occurs, the single variable beliefs provide the exact marginals [START_REF] Weiss | Correctness of belief propagation in gaussian graphical models of arbitrary topology[END_REF]. In fact, for continuous variables, the Gaussian distribution is the only one compatible with the BP rules. Expectation propagation [START_REF] Minka | Expectation propagation for approximate bayesian inference[END_REF] is a way to address more general distributions in an approximate manner.

The inverse Ising Problem

Once the underlying joint probability measure is given, this algorithm can be very efficient for inferring hidden variables, but in real applications it is often the case that we have first to build the model from historical data.

From now on we assume that we have binary variables (x i ∈ {0, 1}). Let {x j i , i ∈ V ⋆ j , j = 1 . . . M } be a set of observations where M represents the number of distinct, but possibly sparse, observations of the system as a whole and V ⋆ j is the set of nodes observed for the jth observation. We can define an empirical measure based on these historical data as

P(x) = 1 M M j=1 1 2 M -|V ⋆ j | i∈V ⋆ j 1 1 {xi=xi} .
As such this measure is of no use for inference and we have to make some hypothesis to find a suitable inference model. There are of course various possibilities, but a simple one is to consider that the mean and the covariance are given for respectively each variable i and each pair of variable (i, j):

mi def = 1 j 1 1 {i∈V ⋆ j } j,V ⋆ j ∋i (2 * xj i -1) χij def = 1 k 1 1 {(i,j)⊂V ⋆ k } k,V ⋆ k ⊃(i,j) (2 * xk i -1)(2 * xk j -1) -mi mj .
Let us introduce also the notation for the joint expectation of pairs of spins:

mij def = Ê(s i s j ) = χij + mi mj .
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, that is an Ising model for binary variables with s i def = 2x i -1:

P(s) = 1 Z[J, h] exp i h i s i + i,j J ij s i s j
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. They are obtained as minimiser's of the dual optimization problem, namely

(h ⋆ , J ⋆ ) = argmin (h,J) log Z[h, J] - i h i mi - ij J ij mij (1.4)
which corresponds to invert the linear response equations:

∂ log Z ∂h i [h, J] = mi (1.5) ∂ log Z ∂J ij [h, J] = mij , (1.6)
since mi and mij are given as input to the model. 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 based on observation and the Ising model:

D KL [ P P] = log Z[h, J] - i h i mi - i<j J ij mij -S( P).
The set of equations (1.5,1.6) 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].

• A common approach is 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 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. • A second possibility is to assume that relevant coupling J ij have locally a tree-like structure. The Bethe approximation mentioned in the previous section is then used with possibly loop corrections. Again this corresponds to having a weakly correlated uni-modal probability measure and these kind of approached are referred as pseudo-moment matching methods in the literature for the reason explained in the previous section. For example the approach 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][START_REF] Mezard | Constraint satisfaction problems and neural networks: A statistical physics perspective[END_REF][START_REF] Yasuda | Approximate learning algorithm in boltzmann machines[END_REF] are based on this assumptions. • In the case where a multi-modal distribution is expected, then a model with many attraction basin is to be found and Hopfield like model [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 more relevant in this case.

Gibbs free energy:

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 . (1.7)
Note that by duality we have

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

∂ 2 G ∂m i ∂m j = - ∂ 2 F ∂h i ∂h j -1 = χ -1 ij .
(1.9)

i.e. the inverse susceptibility matrix. Finding a set of J ij satisfying this last relation along with (1.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 .
(1.10)

Plefka's expansion:

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] (1.11)
= G 0 [m] + ∞ n=0 α n n! G n [m, J] (1.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 is always fulfill (1.7). This precisely is the Plefka's expansion, and each term of the expansion (1.12) can be obtained by successive derivation of (1.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 (1.7), the first derivative of (1.11) w.r.t α gives dG[m, αJ] dα = -E α H J ,
while the second reads:

d 2 G[m, αJ] dα 2 = -E c α H 2 J - i dh i (m, α) dα E c α H J s i .
In these expressions, it is the connected part of the expectation, noted

E c [XY ] def = E[XY ] -E[X]E[Y ],
which appears when deriving on the free energy. To get successive derivative of h(m, α) one can use (1.8). Another possibility is to express the fact that m is fixed,

dm i dα = 0 = - d dα ∂F [h(α), αJ] ∂h i = i,j h ′ j (α)E c α (s i s j ) + E c α (H J s i ), giving h ′ i (α) = - j [χ -1 ] ij E c α (H J s j ).
To get the first two terms in the Plefka's expansion we need to compute these quantities at α = 0:

E c H 2 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 ), E c 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 reads:

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 high-temperature expansions[END_REF].

Linear response approximate solution:

At this point we are in position to find an approximate solution to the inverse Ising problem, either by inverting equation (1.9) or (1.10). To get a solution at a given order n in the coupling, solving (1.10) requires G at order n + 1, while it is needed at order n in (1.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 (1.10) leads directly for the basic mean-field solution to:

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

At this level of approximation for G, using (1.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 correspond precisely to the TAP equations. Using now (1.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 , (1.13)
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 approximation:

In this case we remark first that when the graph formed by the observed correlations χij is a tree then the form (1.2) of the joint probability corresponding to the Bethe approximation yields actually an exact solution to the inverse problem (1.4)

P(x) = i<j pij (x i , x j ) p(x i )p(x j ) i pi (x i ),
where the p are the single and pair variables empirical marginal given by the observations. Rewriting this expression as an Ising model results in the following expressions for the parameters

h i = 1 -d i 2 log p1 i p0 i + 1 4 j∈i log p11 ij p10 ij p01 ij p00 ij , (1.14)
J ij = 1 4 log p11 ij p00 ij p01 ij p10 ij , (1.15)
while the partition function simply reads

Z Bethe [p] = exp - 1 4 ij log p00 ij p01 ij p10 ij p11 ij - i 1 -d i 2 log(p 0 i p1 i ) , (1.16)
and where the p's are parametrized as

pxi i def = p(x i = 1 + s 2 ) = 1 2 (1 + m i s), (1.17)
pxixj ij def = p(x i = 1 + s 2 , x j = 1 + s ′ 2 ) = 1 4 (1 + m i s + m j s ′ + m ij ss ′ ) (1.18)
are the empirical frequency statistics given by the observations for m ≡ m.

The corresponding Gibbs free energy can then be written explicitly using (1.14,1.15,1.16). With fixed magnetization's 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 (1.15). For the linear response, we get from (1.14):

∂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 (1.15), 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 [2]

[ χ-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 . (1.19)
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:

• [ χ-1 ] ij = 0, then this self-consistency relation has to hold and J ij is given by (1.15) using (1.15) does not hold in that case.

χ ij = χij . • [ χ-1 ] ij = 0 then J ij = 0 while χij can be non-zero, because
Finally complete consistency of the solution is checked on the diagonal elements in (1.19). If full consistency is not verified, these equation can nevertheless be used to find approximate solutions. Remark that if we restrict the set of equations (1.19), 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 [START_REF] Mezard | Constraint satisfaction problems and neural networks: A statistical physics perspective[END_REF][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] Chau 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 (1.15) and (1.19):

χ 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 ,
(1.20) while the corresponding computed of ξ ij , instead of the observed one ξij , has to be inserted in (1.14) to be fully consistent. Note that J Bethe ij and J T AP ij coincide at second order in [ χ-1 ] ij .

Application context

Road Traffic Inference

Once the underlying joint probability measure is given, the BP algorithm can be very efficient for inferring hidden variables, but in real applications it is often the case that we have first to build the model. This is precisely the case for the application that we are considering concerning the reconstruction and prediction of road traffic conditions, typically on the secondary network from sparse observations. Existing solutions for traffic information are classically based on data coming from static sensors (magnetic loops) on main arterial roads. These devices are far too expensive to be installed everywhere on the traffic network and other sources of data have to be found. One recent solution comes from the increasing number of vehicles equipped with GPS and able to exchange data through cellular phone connections for example, in the form of so-called Floating Car Data (FCD). Our objective in this context is to build an inference schema adapted to these FCD, able to run in real time and adapted to large scale road networks, of size ranging from 10 3 to 10 5 road segments. In this respect, the BP algorithm seems very well suited, but the difficulty is to construct a model based on these FCD. To set an inference schema, we assume that a large amount of FCD sent by probe vehicles concerning some area of interest are continuously collected over a reasonable period of time (one year or so) such as to allow a finite fraction (a few percents) of road segments to be covered in real time. Schematically the inference method works as follows:

• Historical FCD are used to compute empirical dependencies between contiguous segments of the road network. • These dependencies are encoded into a graphical model, which vertices are (segment,timestamps) pairs attached with a congestion state, i.e. typically congested/not-congested. • Congestion probabilities of segments that are unvisited or sit in the shortterm future are computed with BP, conditionally to real-time data.

On the factor-graph, the information is propagated both temporally and spatially. In this perspective, reconstruction and prediction are on the same footing, even though prediction is expected to be less precise than reconstruction. For each segments, we may potentially collect a distribution f of travel time and it is not clear how to decide from this distribution, whether a link is congested or not given an newly observed travel time. A straightforward possibility is to consider the mean travel time, or even more robust, the median travel time as a separator of the two state. The way we actually see this encompass this possibility, but is not limited to it. The idea is to define the latent binary state τ (= 1+s 2 ) associated to some travel time x in an abstract way through the mapping:

Λ(x) def = P (τ = 1|x).
This means that an observation x is translated into a conditional probability for the considered segment to be congested. This number Λ(x) ∈ [0, 1], represents our practical definition for the traffic index.

Using Bayes rules and the Boolean notation τ def = 1 -τ , we obtain

P (x|τ ) = Λ(x) p Λ τ + 1 -Λ(x) 1 -p Λ τ f (x). (1.21)
where p Λ def = P (τ = 1). The normalization constraint imposes

p Λ = Λ(x) f (x)dx. (1.22)
A certain amount of information can be stored in this mapping. A special case mentioned before corresponds to having for Λ a step function, i.e.

Λ(x) = 1 1 {x>x * } , (1.23)
with an adjustable parameter corresponding to the threshold x * . Another parameter free possibility is to use the empirical cumulative distribution:

Λ(x) = F (x) def = P (x < x). (1.24)
Now, Given a map Λ, an obvious way to convert back a probability u = P (τ = 1) into a travel time consist then simply in using, when it exits, the inverse map:

x = Λ -1 (u). (1.25)
Actually another legitimate way to proceed is based on the conditional probability (1.21) to yields the following estimator:

x = argmin y E x -y r ,
where the expectation is taken from the probability distribution

P (x) = P (x|τ = 1)u + P (x|τ = 0)(1 -u),
and where x -y r represents the loss function, measuring the error between the prediction x and the actual value y. In this last case, the natural requirement that we seek for Λ is that the mutual information between x and τ be maximal. This reads

I(x, τ ) def = τ ∈{0,1}
dxP (x, τ ) log P (x, τ ) P (x)P (τ ) ,

= dx Λ(x) log Λ(x) + (1 -Λ(x)) log(1 -Λ(x)) f (x) -h(p Λ ), = du h Λ F -1 (u) -h(p Λ ),
after introducing the binary information function h(x)

def = x log x + (1 - x) log(1 -x).
In this form, h being convex, reaching its maximum at x = 0 and x = 1, its minimum at x = 1/2, it is then straightforward to obtain that the step function (1.23) with x * = F -1 (1/2) corresponding to the median observation is the limit function which maximizes I(x, τ ). If instead we use the inverse map Λ -1 , the mutual information between x and τ is not relevant. Without any specific hypothesis on the distribution of beliefs that BP should generate, a simple requirement is then to impose a minimum information i.e. a maximum entropy contained in the variable u = Λ(x), which probability density is given by

dF (u) def = δ u -Λ(x) d F (x) = d F dΛ Λ -1 (u) .
Using this and the change of variable x = Λ -1 (u) yields the entropy

S[u] = -d F (x) d F dΛ (x) = -D KL ( F Λ),
expressed as the opposite of the relative entropy between F and Λ. Without any further constraint, this leads to the fact that Λ = F is the optimal mapping. In both cases, additional constraints comes from the fact that we want a predictor x minimizing a loss function x -x r which depends on the choice of the Euclidean norm L r (see [START_REF] Martin | Encoding dependencies between real-valued observables with a binary latent MRF[END_REF] for details).

Global inference model:

In fact the mapping between real-valued observations and the binary latent states is only one element of the model. The general schema of our Ising based inference model is sketched on Figure 1.2. It can be decomposed into 4 distinct pieces:

• A single variable statistical model translating real-valued observations into binary latent states. • A pairwise statistical model of the dependency between latent states.

• A MRF model to encode the network of dependencies.

• The Belief propagation algorithm to decode a partially observed network.

It is based on a statistical description of traffic data which is obtained by spatial and temporal discretization, in terms of road segments i and discrete time slots t corresponding to time windows of typically a few minutes, leading to consider a set of vertices V = {α = (i, t)}. To each vertex is attached a microscopic degree of freedom x α ∈ E, as a descriptor of the corresponding segment state (e.g. E = {0, 1}, 0 for congested and 1 for fluid). The model itself is based on historical data in form of empirical marginal distributions p(x α ), p(x α , x β ), giving reference states and statistical interactions between degrees of freedom. Finally, reconstruction and prediction are produced in the form of conditional marginal probability distribution p(x α |V * ) of hidden variables in V\V * , conditionally to the actual state of the observed variables in the set V * .

In addition to this microscopic view, it is highly desirable to enrich the description with macroscopic variables, able in particular to capture and encode the temporal dynamics of the global system. These can be obtained by some linear analysis, e.g. PCA or with non-linear methods of clustering providing possibly hierarchical structures. Once some relevant variables are identified, we can expect to have a macroscopic description of the system, which can potentially be easily coupled to the microscopic one, by adding some nodes into the factor-graph. These additional degrees of freedom would be possibly interpreted in terms of global traffic indexes, associated to regions or components.

The binary latent states are used to model the interactions in a simplified way enabling for large scale applications. Trying to model exactly the pairwise dependencies at the observation level is potentially too expensive from the statistical as well as the computational viewpoint. So the pairwise model sketched on Figure 1.2 corresponds to

P (x i , x j ) = τ,τ ′ pij (τ, τ ′ )P (x i |τ )P (x j |τ ′ ),
with P (x|τ ) given in (1.21) and pij to be determined from empirical frequency statistics. Since a probability law of two binary variables requires three independent parameters; two of them are already being given by individuals marginals probabilities p1 i def = P (τ i = 1) according to (1.22). For each pair of variables, one parameter remains therefore to be fixed. By convenience we consider the coefficient

p 11 ij def = P (τ i = 1, τ j = 1),
and write a moment matching constraint in the traffic index space2 . We obtain

p11 ij = p1 i p1 j + cov Λ i (x i ), Λ j (x j ) 2 Ê[Λ i (x)] -1 2 Ê[Λ j (x)] -1 ,
involving the empirical expectation of indexes, Ê[Λ i (x)] and empirical covariance between indexes cov Λ i (x i ), Λ j (x j ) obtained from observation data.

MRF model and pseudo moment matching calibration

At the microscopic level, the next step is to define the MRF i.e. the Ising model, on which to run BP with good inference properties. Recall that we try to answer two related questions:

• Given the set of coefficients p(τ i,t ) and p(τ i,t , τ j,t ), considered now as model input, what is the joint law P ({τ i,t , (i, t) ∈ V})? • Given actual observations {x * i,t , (i, t) ∈ V * }, how to infer {x i,t , (i, t) ∈ V\V * }?

The solution that we have been exploring [START_REF] Furtlehner | A belief propagation approach to traffic prediction using probe vehicles[END_REF] is based on the the Bethe approximation described in Section 1.4.2. It consists to use the Bethe approximation (1.2) for the encoding and the belief-propagation for the decoding, such that the calibration of the model is coherent with the inference algorithm. In particular, when there is no real time observation, the reference point is given by the set of historical beliefs, so we expect that running BP on our MRF delivers precisely these beliefs. Stated differently, we look for the φ and ψ defining the MRF in (1.1) such that the beliefs match the historical marginals:

b i (τ i ) = pi (τ i ), and b ij (τ i , τ j ) = pij (τ i , τ j ).
As explained in Section 1.3 there is an explicit solution to this problem, because BP is coherent with the Bethe approximation [START_REF] Yedidia | Generalized belief propagation[END_REF], and thus any BP fixed point b has to verify

P(τ ) = i∈V φ i (τ i ) (i,j)∈F ψ ij (τ i , τ j ) ∝ i∈V b i (τ i ) (i,j)∈F b ij (τ i , τ j ) b i (τ i )b i (τ j ) . (1.26)
As a result, a canonical choice for the functions φ and ψ is simply

φ i (τ i ) = pi (τ i ), ψ ij (τ i , τ j ) = pij (τ i , τ j ) pi (τ i )p j (τ j ) , (1.27)
along with m i→j (x j ) ≡ 1 as a particular BP fixed point. In addition, from the re-parametrization property of BP [START_REF] Wainwright | Stochastic processes on graphs with cycles: geometric and variational approaches[END_REF], any other choices verifying (1.26) produces the same set of fixed points with the same convergence properties. Note that more advanced methods than the strict use of the Bethe approximation, presented in the preceding section could be used as well, but as we shall see in the next sections, the hypothesis that traffic data could be well represented by one single BP fixed point might not be fulfilled. In that case the linear response, which takes a BP fixed point as a reference starting point might be of limited efficiency. Instead of trying to use more accurate version of the linear response, we have followed a different route, by enriching the Bethe approximation with an adjustable parameter, interpreted as an inverse temperature, in order to better calibrate the model in a multiple BP fixed point context. This will be explained in the next section.

Next, for the decoding part, inserting information in real time in the model is done as follows. In practice, observations are in the form of real numbers like speed or travel time. One possibility is to project such an observation onto the binary state τ i = 0 or τ i = 1, but this proves to be too crude. As explained in Section 1.4.2, since the output of BP is anyway in the form of beliefs, i.e. real numbers in [0, 1], the idea is to exploit the full information by defining a correspondence between observations x i and probabilities p * (τ i = 1). The optimal way of inserting this quantity into the BP equations is obtained variationally by imposing the additional constraint b i (τ i ) = p * (τ i ), which results in modified messages sent from i ∈ V * , now reading [START_REF] Furtlehner | Spatial and temporal analysis of traffic states on large scale networks[END_REF]

n i→j (x i ) = p * i (τ i ) m j→i (τ i ) .
This leads to a new version of BP which convergence properties have been analyzed in [START_REF] Martin | Encoding dependencies between real-valued observables with a binary latent MRF[END_REF]. This works well in practice, in particular when compared to some heuristic method consisting in giving a bias to the local field of the observed variables as shown on the Figure 1.2 discussed in the last Section 1.6.

Multiple BP fixed points for multiple traffic patterns

Some experiments with a preliminary version of this procedure [START_REF] Furtlehner | A belief propagation approach to traffic prediction using probe vehicles[END_REF] indicate that many BP fixed point can exist in absence of information, each one corresponding to some congestion pattern like e.g. congestion/free flow. We have analyzed in [START_REF] Furtlehner | Learning multiple belief propagation fixed points for real time inference[END_REF] the presence of multiple fixed points by looking at a study case, and we outline some of the results in this section. In this study, we considered a generative hidden model of traffic in the form of a probabilistic mixture, with each component having a simple product form:

P hidden (τ ) def = 1 C C c=1 i∈V p c i (τ i ). (1.28)
C represents the number of mixture components. Although (1.28) is quite general, the tests are conducted with C ≪ N , with well separated components of the mixture. The single sites probabilities p c i def = p c i (1), corresponding to each component c, are generated randomly as i.i.d. variables,

p c i = 1 2 (1 + tanh h c i )
with h c i uniformly distributed in some fixed interval [-h max , +h max ]. The mean of p c i is therefore 1/2 and its variance reads

v def = 1 4 E h tanh 2 (h) ∈ [0, 1/4].
This parameter v implicitly fixed by h max represents the average level of "polarizability" of the variables in each cluster: v = 0 corresponds to p c i = 1/2 while v = 1/4 corresponds to p c i ∈ {0, 1} with even probabilities. The interpretation of this model is that traffic congestion is organized in various patterns, which can show up at different times. We then studied the behavior of our inference model on the data generated by this hidden probability by adding a single parameter α into its definition (1.27):

φ i (τ i ) = pi (τ i ), ψ ij (τ i , τ j ) = pij (τ i , τ j ) pi (τ i )p j (τ j ) α (1.29)
where pi and pij are again the 1-and 2-variable frequency statistics that constitute the input of the model, while (1.28) is assumed to be unknown. This parameter α, which can be interpreted as an inverse temperature in the Ising model, is there to compensate for saturation effects (α < 1), when the coupling between variables are too large. This is due to some over-counting of the dependencies between variables which may occur in a multiply connected graph. Still, in complement, some sparsity can be imposed to the factor graph with help of some link selection procedure that reduce the mean connectivity to K.

The typical numerical experiment we perform, given a configuration randomly sampled from (1.28), is to reveal gradually the variables τ V * in a random order and compute conditional predictions for the remaining unknown variables. We then compare the beliefs obtained with the true conditional marginal probabilities P (τ i = τ |τ V * ) computed with (1.28), using an error measure based on the Kullback-Leibler distance:

D KL def = τ ∈{0,1} b i (τ ) log b i (τ ) P (τ i = τ |τ V * ) V * ,
where <> V * means an average taken on the set of hidden variables. A sample test shown on Figure 1.1.b indicates for example that, on a system with 10 3 variables, it is possible with our model to infer with good precision a mixture of 20 components by observing 5% of the variables. To interpret these results, letting s i = 2τ i -1, we first identify the Ising model corresponding to the MRF given by (1.29):

P(s) = 1 Z e -βH[s] ,
with an inverse temperature β and the Hamiltonian

H[s] def = - 1 2 i,j J ij s i s j - i h i s i .
The identification reads:

βJ ij = α 4 log pij (1, 1)p ij (0, 0) pij (0, 1)p ij (1, 0) , βh i = 1 -αK i 2 log pi (1) pi (0) + α 4 j∈i log pij (1, 1)p ij (1, 0) pij (0, 1)p ij (0, 0) ,
Then, in the limit C ≫ 1, N ≫ C and fixed average connectivity K, we get asymptotically a mapping to the Hopfield model [START_REF] Hopfield | Neural network and physical systems with emergent collective computational abilities[END_REF]. The relevant parameters in this limit are η = C/N and the variance v ∈ [0, 1/4] of the variable bias in the components. In this limit, the Hamiltonian is indeed similar to the one governing the dynamics of the Hopfield neural network model:

H[s] = - 1 2N i,j,c ξ c i ξ c j s i s j - i,c h c i ξ c i s i , with ξ c i def = p c i (1) -1 2 √ v and h c i = C 2αK √ v - 2C √ v K j∈i Cov(ξ c i , ξ c j ),
the inverse temperature given by the mapping reads

β = 4αvK C .
Using mean-field methods, the phase diagram of this model has been established [START_REF] Amit | Statistical mechanics of neural networks near saturation[END_REF]. There are 3 phases, separated by the transition lines T g , between the paramagnetic phase and the spin glass phase, and T c , between the spin glass phase and the ferromagnetic phase (see Figure 1.1.a). The latter corresponds to the so-called Mattis states, i.e. to spin configurations correlated with one of the mixture components, of direct relevance w.r.t. inference. Locating the various models obtained in this diagram as on Figure 1.1.a, helps to understand whether inference is possible or not with our MRF model.

We have also tested a multi-parameter version of the model in [START_REF] Furtlehner | Learning multiple belief propagation fixed points for real time inference[END_REF], where the links sorted according to the mutual information they contribute for, and grouped them into a certain number of quartiles: to each quartile q we associated a parameter α q < 1. Using a calibration procedure based on a stochastic optimization algorithm CMAES [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF], we can see on some examples a significant improvement of the model, as seen on the example presented on Figure 1.1.c.

To summarize, the main lessons of this theoretical study are

• the various components of a probabilistic mixture with weak internal correlations maybe correctly accounted for by our inference model. It is able to to associate in an unsupervised way one BP fixed point to each component of the mixture. • the mechanism for that can be understood by some asymptotic analysis which reveals a connection with a Hopfield model, where the main patterns corresponds to the components of the mixture. The phase diagram of the Hopfield model gives then relevant indications on whether inference will be easy, difficult or impossible depending on the ratio N states /N variables and on the mean internal variance of the variables v within each state. • the model can be easily generalized to a multi-parameter version to improve its accuracy with help of a calibration based on a robust optimization strategy like the CMAES algorithm for example.

An example of BP fixed points associated to the mixture's components is given on Figure 1. neighbor (K-NN) predictor, is able to predict correctly real-valued hidden variables.

Experiments with Synthetic and Real Data

Using both synthetic and real data we perform two kind of numerical tests:

(i) reconstruction/prediction experiments (ii) automatic segmentation and BP fixed points identification

In experiments of the type (i), the data set is divided into two parts, one corresponding to the learning set, used to build the model and the other part corresponding to the test set, used to perform the tests. In the reconstruction tests traffic configuration corresponding to a single time layer are extracted from the test set; a fraction ρ of variables are chosen at random to be revealed, and while this fraction is progressively increased, travel time are inferred for the hidden one. In prediction experiments, traffic configuration corresponding to successive time layers are selected, with present time t 0 separating the time layers into 2 equal parts, one corresponding to past and the other to future. Observed variable are necessarily selected in the past window time.

Variables with time stamp t = t0 (present) or t > t0 (future) are inferred, i.e. reconstructed or predicted respectively. In the type (ii) experiments, on one hand an automatic clustering of the data on reduced dimensional space is performed with machine learning techniques [START_REF] Han | Analysis of Network-level Traffic States using Locality Preservative Non-negative Matrix Factorization[END_REF]. On the other hand the BP fixed points obtained at ρ = 0 are listed [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] and compared to segmentation in the reduced dimensional space. A first set of data has been generated with the traffic simulator "METROPO-LIS" [START_REF] De Palma | Real cases applications of the fully dynamic METROPOLIS tool-box: an advocacy for large-scale mesoscopic transportation systems[END_REF], for the benchmark network called Siouxfalls shown on Figure 1.1.a. An example of the automatic clustering of spatial configurations with the corresponding BP fixed points associated to free flow and congestion is shown on Figure 1.1.b. On Figure 1.2.a and b, two ways of inserting information are compared, the variational one mentioned in Section 1.4.3 with an heuristic one based on local fields. In both cases the optimal tuning of α yields two BP fixed points, but the variational method gives better results on reconstruction test and is more consistent with clustering results (same optimal value of α). On Figure 1.2.c a reconstruction experiment on simulated Siouxfalls Metropolis data is shown, using the cumulative distribution for both the encoding (1.24) and inverse decoding (1.25) of traffic indexes. The performance is similar to a K-NN predictor, although much more economical. To set a scale of comparison predictors obtained from historical mean and median values are also shown on the same plot.

To perform tests on real data, we have also considered a dataset consisting of travel times measured every 3 minutes over 2 years of a highway segmented

Conclusion

The work concerning the application of belief propagation and related Boltzmann machine to traffic data is related to some ongoing projects [START_REF]TRAVESTI project[END_REF][START_REF]PUMAS project[END_REF]. It is based on mean-field concepts in physics and basically related to the linear response theory. When combined with machine learning techniques like the automatic segmentation methods it can lead to efficient models able to cope with real-time constraints on large scale networks. We advocate for an Ising model for traffic statistical modelling and propose a proper way for defining traffic indexes which could be also useful for traffic management systems. Still a natural concurrent approach not exposed here can be build analogously using a multivariate model, for which Gaussian belief propagation would apply. More real data will help to decipher from these two possibilities. The main hypothesis underlying our Ising based approach assumes that traffic congestion is well represented by multiple distant pattern superposition. This needs validation with real data on networks. Our reconstruction schema seems to work already with simple underlying binary indexes, but more work is needed for the dynamical part to be able to perform prediction.

Figure 1 .Fig. 1 . 1

 111 Fig. 1.1 Example of factor graph (a) and message propagation rules (b).

1. 4 . 2 Fig. 1 . 1

 4211 Fig. 1.1 Underlying Ising modelling of traffic configurations.

Fig. 1 . 2

 12 Fig. 1.2 Sketch of the Ising based inference schema.

Fig. 1 . 1

 11 Fig. 1.1 (a) Phase diagram of the Hopfield model and optimal points found experimentally.(b) D KL error as a function of observed variables ρ for the single parameter model with N = 1000 and C = 20 and various pruning levels and for the muliparameter model N = 100 C = 5 with various number of calibrated parameters ranging from 1 to 30 (c).

2 .Fig. 1 . 2

 212 Fig. 1.2 Segmentation and BP fixed point identification (a) for synthetic travel time data corresponding to a mixture with five components with internal correlations. Corresponding reconstruction experiment where the L 1 travel time error is plotted against the fraction ρ of observed variables and compared with a K-NN predictor considered here as ground truth.

Fig. 1 . 1 Fig. 1 . 2

 1112 Fig. 1.1 (a) Siouxfalls network. (b) Automatic segmentation of simulated Siouxfalls data and BP fixed point identification projected in the 3-d main PCA space.

Fig. 1 . 3

 13 Fig. 1.3 (a) Automatic segmentation of Highway data projected on 3-d dominant PCA space. (b) Error on travel time for a BP prediction of 3 time layers in future as a function of the fraction of observed variables at t 0 . Comparison is made with a predictor combining recent available observations with historical time dependent mean.

An Ising Model for Road Traffic Inference

Potentially any arbitrary mapping φ(x) could be considered to perform the moment matching

Acknowledgements: This gives me the occasion to express my warm thanks to my colleagues Jean-Marc Lasgouttes and Victorin Martin with whom it is a pleasure to collaborate on the main subjects discussed in this review. I am also grateful to Anne Auger, Yufei Han, Fabrice Marchal and Fabien Moutarde for many aspects mentioned in this work concerning ongoing projects. This work was supported by the grant ANR-08-SYSC-017 from the French National Research Agency.