The connection between the statistical physics of disordered systems and optimization problems in computer science dates back from twenty years at least [START_REF] Mézard | Spin glasses and beyond[END_REF]. After all zero temperature statistical physics is simply the search for the state with minimal energy, while the main problem in combinatorial optimization is to look for the configurations of parameters minimizing some cost function (the length of a tour in the traveling salesman problem (TSP), the number of violated constraints in constrained satisfaction problems, ...) [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF]. Yet, despite the beautiful studies of the average properties of the TSP, Graph partitioning, Matching, ..., based on the recently developed mean-field spin glass theory [START_REF] Mézard | Spin glasses and beyond[END_REF], a methodological gap between the fields could not be bridged [30]. In statistical physics statements are usually made on the properties of samples given some quenched disorder distribution such as the typical number of solutions, minimal energy ... In optimization, however, one is interested in solving one (or several) particular instances of a problem, and needs efficient ways to do so, that is, requiring a computational effort growing not too quickly with the number of data defining the instance. Knowing precisely the typical properties for a given, academic distribution of instances does not help much to solve practical cases.

At the beginning of the nineties practitionners in artificial intelligence realized that classes of random constrained satisfaction problems used as artificial benchmarks for search algorithms exhibited abrupt changes of behaviour when some control parameter were finely tuned [START_REF] Mitchell | Proc. of the Tenth Natl. Conf. on Artificial Intelligence (AAAI-92)[END_REF]. The most celebrated example was random K-Satisfiability, where one looks for a solution to a set of random logical constraints over a set of Boolean variables. It appeared that, for large sets of variables, there was a critical value of the number of constraints per variable below which there almost surely existed solutions, and above which solutions were absent. An important feature was that search algorithms performances drastically worsened in the vicinity of this critical ratio.

This phenomenon, strongly reminiscent of phase transitions in condensed matter physics, led to a revival of the interface between statistical physics and computer science, which has not vanished yet. The purpose of the present lecture is to introduce the non specialist reader to the concepts and techniques required to understand the literature in the field. For the sake of simplicity the presentation will be limited to one computational problem, namely, linear systems of Boolean equations. A good reason to do so is that this problem concentrates most of the features encountered in other optimization problems, while being technically simpler to study. In addition it is closely related to error-correcting codes in communication theory, see lectures by A. Montanari and R. Urbanke in the present book. Extension to other problems will be mentioned in the conclusions.

The lecture is divided into three parts. Sections 1 and 2 are devoted to the presentation of the model and of elementary concepts related to phase transitions e.g. finite-size scaling, large deviations, critical exponents, symmetry breaking, ... Sections 3 and 4 expose the specific statistical mechanics techniques and concepts developed in disordered systems to deal with highly interacting and random systems, namely the replica and cavity approaches. Finally Section 5 focuses on dynamics and the study of search algorithms.

B. Linear systems of Boolean equations

Linear systems of Boolean equations look very much like their well known counterparts for integer-valued variables, except that equalities are defined modulo two. Consider a set of N Boolean variables x i with indices i = 1, . . . , N . Any variable shall be False (F) or True (T). The sum of two variables, denoted by +, corresponds to the logical exclusive OR between these variables defined through,

F + T = T + F = T , F + F = T + T = F . (1)
In the following we shall use an alternative representation of the above sum rule. Variables will be equal to 0 or 1, instead of F or T respectively. Then the + operation coincides with the addition between integer numbers modulo two.

The following is a linear equation involving three variables,

x 1 + x 2 + x 3 = 1 . (2)
Four among the 2 3 = 8 assignments of (x 1 , x 2 , x 3 ) satisfy the equation: (1, 0, 0), (0, 1, 0), (0, 0, 1) and (1, 1, 1). A Boolean system of equations is a set of Boolean equations that have to be satisfied together. For instance, the following Boolean system involving four variables

   x 1 + x 2 + x 3 = 1 x 2 + x 4 = 0 x 1 + x 4 = 1 (3)
has two solutions: (x 1 , x 2 , x 3 , x 4 ) = (1, 0, 0, 0) and (0, 1, 0, 1). A system with one or more solutions is called satisfiable.

A trivial example of an unsatisfiable Boolean system is

x 1 + x 2 + x 3 = 1 x 1 + x 2 + x 3 = 0 . ( 4
)
Determining whether a Boolean system admits an assignment of the Boolean variables satisfying all the equations constitutes the XORSAT (exclusive OR Satisfaction) problem. In the following, we shall restrict for some reasons to be clarified in Section II to K-XORSAT, a variant of XORSAT where each Boolean equation include K variables precisely. K-XORSAT belongs to the class P of polynomial problems [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF]. Determining whether a system is satisfiable or not can be achieved by the standard Gaussian elimination algorithm in a time (number of elementary operations) bounded from above by some constant times the cube of the number of bits necessary to store the system[65] [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF].

If the decision version of K-XORSAT is easy its optimization version is not. Assume you are given a system F , run the Gauss procedure and find that it is not satisfiable. Determining the maximal number M S (F ) of satisfiable equations is a very hard problem. Even approximating this number is very hard. It is known that there is no approximation algorithm (unless P=NP) for XORSAT with ratio r > 1 2 , that is, guaranteed to satisfy at least r × M S (F ) equations for any F . But r = 1 2 is achieved, on average, by making a random guess[66]!

C. Models for random systems

There are many different ways of generating random Boolean systems. Perhaps the simplest one is the following, called fixed-size ensemble. To build an equation we pick up uniformly at random K distinct indices among the N ones, say, i 1 , i 2 and i k . Then we consider the equation

x i1 + x i2 + . . . + x i k = v . (5)
The second member, v, is obtained by tossing a coin: v = 0 or v = 1 with equal probabilities (one half) and independently of the indices of the variables in the first member. The process is repeated M times, without correlation between equations to obtain a system with M equations. Another statistical ensemble is the fixed-probability ensemble. One scans the set of all H = 2 N K equations one after the other. Each equation is added to the system with probability p, discarded with probability 1p. Then a system with, on average, p H equations (without repetition) is obtained. In practice one chooses p = M H to have the same (average) number of equations as in the fixed-size ensemble.

The above distributions are not the only possible ones. However they are easy to implement on a computer, are amenable to mathematical studies, and last but not least, lead to a surprisingly rich phenomenology. One of the key quantities which exhibits an interesting behaviour is P SAT (N, α) = Probability that a system of random K-XORSAT with N variables and M = α N equations is satisfiable , which obviously depends on K and the statistical ensemble. Given N P SAT is a decreasing function of α. We will see that, in the infinite size limit (and for K ≥ 2), the decrease is abrupt at some well defined ratio, defining a phase transition between Satisfiable and Unsatisfiable phase [18]. The scope of the lecture is to give some tools to understand this transition and some related phenomena.

II. BASIC CONCEPTS: OVERVIEW OF STATIC PHASE TRANSITIONS IN K-XORSAT

In this Section we introduce the basic concepts necessary to the study of random K-XORSAT. It turns out that even the K = 1 case, trivial from a computer science point of view (each equation contains a single variable!), can be used as an illustration to important concepts such as scaling and self-averageness. Ideas related to the percolation phase transition and random graphs are illustrated on the K = 2 case. Finally the solution space of 3-XORSAT model exemplifies the notion of clusters and glassy states.

A. Finite-size scaling (I): scaling function Figure 1(left) shows the probability P SAT that a randomly extracted 1-XORSAT formula is satisfiable as a function of the ratio α, and for sizes N ranging from 100 to 1000. We see that P SAT is a decreasing function of α and N .

Consider the subformula made of the n i equations with first member equal to x i . This formula is always satisfiable if n i = 0 or n i = 1. If n i ≥ 2 the formula is satisfiable if and only if all second members are equal (to 0, or to 1), an event with probability ( 12 ) ni-1 decreasing exponentially with the number of equations. Hence we have to consider the following variant of the celebrated Birthday problem [67]. Consider a year with a number N of days, how should scale the number M of students in a class to be sure that no two students have the same birthday date?

p = M-1 i=0 1 - i N = exp - M (M -1) 2N + O(M 3 /N 2 ) . (6)
Hence we expect a cross-over from large to small p when M crosses the scaling regime √ N . Going back to the 1-XORSAT model we expect P SAT to have a non zero limit value when the number of equations and variables are both sent to infinity at a fixed ratio y = M/ √ N . In other words, random 1-XORSAT formulas with N variables, M equations or with, say, 100 × N variables, 10 × M equations should have roughly the same probabilities of being satisifiable. To check this hypothesis we replot the data in Figure 1 after multiplication of the abscissa of each point by √ N (to keep y fixed instead of α). The outcome is shown in the right panel of Figure 1. Data obtained for various sizes nicely collapse on a single limit curve function of y.

The calculation of this limit function, usually called scaling function, is done hereafter in the fixed-probability 1-XORSAT model where the number of equations is a Poisson variable of mean value M = y √ N . We will discuss the equivalence between the fixed-probability and the fixed-size ensembles later. In the fixed-probability ensemble the numbers n i of occurence of each variable x i are independent Poisson variables with average value M /N = y/ √ N . Therefore the probability of satisfaction is

P p SAT (N, α = y √ N ) =   e -y/ √ N   1 + n≥1 (y/ √ N ) n n! 1 2 n-1     N = 2e -y/(2 √ N ) -e -y/ √ N N , (7)
where the p subscript denotes the use of the fixed-probability ensemble. We obtain the desired scaling function

Φ 1 (y) ≡ lim N →∞ ln P p SAT (N, α = y √ N ) = - y 2 4 , (8)
in excellent agreement with the rescaled data of Figure 1 (right) [19].

B. Self-averageness of energy and entropy

Let us now consider random 1-XORSAT formulas at a finite ratio α, and ask for the distribution of the minimal fraction of unsatisfied equations, hereafter called ground state (GS) energy e GS . For simplicity we work in the fixedprobability ensemble again. The numbers n 0 i , n 1 i of, respectively, x i = 0, x i = 1 are independent Poisson variables with mean α 2 . The minimal number of unsatisfied equations is clearly min(n 0 i , n 1 i ). The GS energy is the sum (divided by M ) of N such i.i.d. variables; from the law of large number it almost surely converges towards the average value where I ℓ denotes the ℓ th modified Bessel function. In other words almost all formulas have the same GS energy in the infinite N limit, a property called self-averageness in physics, and concentration in probability.

e GS (α) = 1 2 1 -e -α I 0 (α) -e -α I 1 (α) , (9)
How many configurations of variables realize have minimal energy? Obviously a variable is free (to take 0 or 1 value) if n 0 i = n 1 i , and is frozen otherwise. Hence the number of GS configurations is N = 2 N f where N f is the number of free variables. Call

ρ = n≥0 e -α α 2 n 1 (n!) 2 = e -α I 0 (α) (10)
the probability that a variable is free. Then N f is a binomial variable with parameter ρ among N ; it is sharply concentrated around N f = ρ N with typical fluctuations of the order of N 1/2 . As a consequence, the GS entropy per variable, s GS = (log N )/N , is self-averaging and almost surely equal to its average value s GS = ρ log 2. Self-averageness is the very useful property. It allows us to study the average value of a random variable, instead of its full distribution. We shall use it in Section III and also in the analysis of algorithms of Section V E. This property is not restricted to XORSAT but was proven to hold for the GS energy [START_REF] Broder | Proceedings of Symposium of Discrete Algorithms (SODA)[END_REF] and entropy [START_REF] Montanari | Proceedings of Symposium of Discrete Algorithms (SODA)[END_REF] of other optimization problems.

Not all variables are self-averaging of course. A straightforward example is the number N of GS configurations itself. Its q th moment reads N q = (1ρ + ρ 2 q )

N where the overbar denotes the average over the formulas. We see that N q ≫ (N ) q : N exhibits large fluctuations and is not concentrated around its average. Very rare formulas with atypically large number N f of free variables contribute more to the q th moment than the vast majority of formulas, and spoil the output. This is the very reason we will need the introduction of the replica approach in Section III.

C. Large deviations for PSAT (I): 1-XORSAT

As we have seen in the previous sections 1-XORSAT formulas with a finite ratio α are unsatifiable with high probability i.e. equal to unity in the infinite N limit. For finite but large N there is a tiny probability that a randomly extracted formula is actually satisifiable. A natural question is to characterize the 'rate' at which P SAT tends to zero as N increases (at fixed α). Answering to such questions is the very scope of large deviation theory (see A for an elementary introduction). Looking for events with very small probabilities is not only interesting from an academic point of view, but can also be crucial in practical applications. We will see in Section V C that the behaviour of some algorithms is indeed dominated by rare events.

Figure 2 shows minus the logarithm of P SAT , divided by N , as a function of the ratio α and for various sizes N . Once again the data corresponding to different sizes collapse on a single curve, meaning that

P SAT (N, α) = e -N ω1(α)+o(N ) . ( 11
)
Decay exponent ω 1 is called rate function in probability theory. We can derive its value in the fixed-probability ensemble from [START_REF] Altarelli | [END_REF] with y = α × √ N , with the immediate result

ω p 1 (α) = α -ln 2 e α/2 -1 . (12)
The agreement with numerics is very good for small ratios, but deteriorates as α increases. The reason is simple. In the fixed-probability ensemble the number M of equations is not fixed but may fluctuate around the average value ,

Ω(α|α) = α -α -α ln(α/ α) . ( 13
)
However the probability that a random 1-XORSAT formula with M equations is satisfiable is also exponentially small in N , with a rate function ω 1 (α) increasing with α. Thus, in the fixed-probability ensemble, a trade-off is found between ratios α close to α (formulas likely to be generated) and close to 0 (formulas likely to be satisfiable). As a result the fixed-probability rate function is

ω p 1 (α) = min α ω 1 (α) + Ω(α|α) , (14)
and is smaller than ω 1 (α). It is an easy check that the optimal ratio α * = α/(2e -α/2 ) < α as expected. Inverting [START_REF] Chao | [END_REF] we deduce the rate function ω 1 in the fixed-size ensemble, in excellent agreement with numerics (Figure 2). This example underlines that thermodynamically equivalent ensembles have to be considered with care as far as rare events are concerned. Remark that, when α → 0, α = α + O(α 2 ), and

ω p 1 (α) = ω 1 (α) + O(α 3
). This common value coincides with the scaling function -Φ 1 (α) [START_REF] Amit | Modeling Brain Function[END_REF]. This identity is expected on general basis (Section II F) and justifies the agreement between the fixed-probability scaling function and the numerics based on the fixed-size ensemble (Figure 1, right).

D. Percolation in random graphs

Though 1-XORSAT allowed us to understand some general features of random optimization problems it is very limited due to the absence of interactions between variables. A more interesting problem is 2-XORSAT where every equation define a joint constraint on two variables. Formulas of 2-XORSAT can be represented by a graph with N vertices (one for each variable), and αN edges. To each equation of the type x i + x j = e corresponds an edge linking vertices i and j, and carrying 0 or 1 label (the value e of the second member). Depending on the input model chosen (Section I C) multiple edges are present or not.

As the formula is random so is graph. Figure 3 shows examples of graphs obtained for various values of α. Notice the qualitative change of structure of graphs when the ratio α varies from low values (graphs are mostly made of small isolated trees) to higher ones (a large part of vertices are now connected together). This change is known as the percolation transition in physics, or the appearance of a giant component in mathematics literature.

Before reviewing some of the aspects of the percolation transition let us mention an important fact on the valency of vertices. As a result of the randomness of the graph generation process, each node share edges with a variable number of neighboring vertices. In the large N limit the degree v of a vertex, i.e. the number of its neighbors, is a Poisson variable with mean 2α,

Proba[v] = e -2α (2α) v v! . ( 15
)
For instance the fraction of isolated vertices is e -2α . The average degree of a vertex, c = 2α, is called connectivity.

It is natural to decompose the graphs into its connected subgraphs, called components. Erdös and Rényi were able in 1960 to characterize the distribution of sizes of the largest component [START_REF] Bollobas | Random Graphs[END_REF],

A B C

FIG. 3: Examples of random graphs generated at fixed number M = αN of edges (fixed-size model without repetition). All graph include N = 20 vertices (grey dots). The average degree of valency, 2α, is equal to 0.5 (A), 1 (B), and 2 (C). The labels of the vertices have been permuted to obtain planar graphs, i.e. avoid crossing of edges.

• When c < 1, the largest component includes ∼ ln N/(c-1-ln c) vertices with high probability. Most components include only a finite number of vertices, and are trees i.e. contain no circuit.

• For c = 1 the largest component contain O(N 2/3 ) vertices.

• When c > 1 there is one giant component containing ∼ γ(c)N vertices; the others components are small i.e. look like the components in the c < 1 regime. The fraction of vertices in the giant component is the unique positive solution of

1 -γ = e -c γ . ( 16
)
It is a non analytic function of c, equal to 0 for c ≤ 1, and positive above, tending to unity when c increases.

The phenomenon taking place at c = 1 is an example of (mean-field) percolation transition. We now give a handwaving derivation of (16). Consider a random graph G over N vertices, with connectivity c. Add a new vertex A to the graph to obtain G ′ . If we want G ′ to be drawn from the same distribution as G, a number v of edges must be attached to A, where v an integer-valued random number following the Poisson distribution (15). After addition of A, some connected components of G will merge in G ′ . In particular, with some probability p v , A will be part of the giant component of G ′ . To estimate p v , we note that this event will not happen if and only if none of the v neighbors of A in G ′ belongs to the giant component of G. Thus,

1 -p v = (1 -γ) v , ( 17
)
where γ is the size (fraction of vertices) of the giant component. Summing both sides of (17) over the distribution (15) for v, and asserting that the change in size of the giant component between G and G ′ is o(1) for large N , we obtain (16).

The above derivation illustrates an ubiquitous idea in probability and statistical physics, which could be phrased as follows: 'if a system is very large, its statistical properties should be, in some sense, unaffected by a small increase in size'. This idea will be useful, in a more sophisticated context, in Section IV.

E. Sat/Unsat transition in 2-XORSAT Figure 4 shows the probability P SAT that a randomly extracted 2-XORSAT formula is satisfiable as function of α, and for various sizes N . It appears that P SAT drops quickly to zero for large N when α reaches the percolation threshold α c = 1 2 . For ratios smaller than α c the probability of satisfaction is positive, but smaller than unity. Take α < 1 2 . Then the random graph G associated to a random 2-XORSAT formula is non percolating, and made of many small components. Identical components (differing only by a relabelling of the variables) may appear several times, depending on their topology. For instance consider a connected graph G ′ made of E edges and V vertices. The average number of times G ′ appears in G is a function of E and V only,

N E,V = N V 2α N E 1 - 2α N V (V -1) 2 +V (N -V ) (18)
since any vertex in G ′ can establish edges with other vertices in G ′ , but is not allowed to be connected to any of the N -V outside vertices. When N is very large compared to E, V we have 2-XORSAT FIG. 4: Probability that a random 2-XORSAT formula is satisfiable as a function of the ratio α of equations per variable, and for various sizes N . The full line is the asymptotic analytical formula [START_REF] Dembo | Large deviations techniques and applications[END_REF].

N E,V ≃ N V -E (2α) E V ! e -2α V . ( 19
)
Three cases should distinguished, depending on the value of V -E:

• V -E = 1: this is the largest value compatible with connectedness, and corresponds to the case of trees. From (19) every finite tree has of the order of N copies in G.

• V -E = 0: this correspond to trees with one additional edge, that is, to graphs having one cycle (closed loop).

The average number of unicyclic graphs is, from (19), finite when N → ∞.

• V -E ≤ -1: the average number of components with more than one cycle vanishes in the large N limit; those graphs are unlikely to be found and can be ignored [69].

Obviously a 2-XORSAT formula with tree structure is always satisfiable [70]. Hence dangerous subformulas, as far as satisfiability is concerned, are associated to unicyclic graphs. A simple thought shows that a unicyclic formula is satisfiable if and only if the number of edges carrying label 1 along the cycle is even. Since the values attached to the edges (second members in the formula) are uncorrelated with the topology of the subgraph (first members) each cycle is satisfiable with probability one half. We end up with the simple formula

P SAT (N, α) = 2 -C(G) (20)
where C(G) denotes the number of cycles in G, and . the average over G. For a reason which will become clear below let us classify cycles according to their length L. How many cycles of length L can we construct? We have to choose first L vertices among N , and join them one after the order according to some order. As neither the starting vertex nor the direction along the cycle matter, the average number of L-cycles is

N L = N (N -1) . . . (N -L + 1) 2L × 2α N L → Λ L = (2α) L 2L . ( 21
)
when N → ∞. As the emergence of a cycle between L vertices is a local event (independent of the environment) we expect the number of L-cycles to be Poisson distributed in the large N limit with parameter Λ L . This statement can actually be proven, and extended to any finite collection of cycles of various lengths [START_REF] Bollobas | Random Graphs[END_REF]: in the infinite size limit, the joint distribution of the numbers of cycles of lengths 1, 2, . . . , L is the product of Poisson laws with parameters Λ 1 , Λ 2 , . . . , Λ L calculated in (21). The probability of satisfaction (20) therefore converges to lim

N →∞ P SAT (N, α) = L≥L0    C≥0 e -ΛL (Λ L /2) C C!    = L≥L0 e -ΛL/2 (22)
where L 0 is the minimal cycle length. In normal random graphs L 0 = 3 since triangles are the shortest cycles. However in our 2-XORSAT model any equation, or more precisely, any first member can appear twice or more, hence L 0 = 2. We conclude that [19] lim

N →∞ P SAT (N, α) = e α/2 (1 -2α) 1 4 when α < α c = 1 2 . ( 23
)
The agreement of this result with the large size trend coming out from numerical simulations is visible in Figure 4.

As P SAT is a decreasing function of α it remains null for all ratios larger than α c . The non analyticity of P SAT at α c locates the Sat/Unsat phase transition of 2-XORSAT.

It is an implicit assumption of statistical physics that asymptotic results of the kind of [START_REF] Dembo | Large deviations techniques and applications[END_REF], rigorously valid in the N → ∞ limit, should reflect with good accuracy the finite but large N situation. An inspection of Figure 4 shows this is indeed the case. For instance, for ratio α = .3, (23) cannot be told from the probability of satisfaction measured for formulas with N = 100 variables. This statement does not hold for α = .4, where the agreement between infinite size theory and numerics sets in when N = 1000 at least. It appears that such finite-size effects become bigger and bigger as α gets closer and closer to the Sat/Unsat threshold. This issue, of broad importance in the context of phase transitions and the pratical application of asymptotic results, is studied in Section II H. Consider ratios α > α c . The giant components of the corresponding formulas contain an extensively large number of independent cycles, so we expect from (20) that the probability of satisfaction is exponentially small in N , P SAT = exp(-N ω 2 (α) + o(N )). Lower and upper bounds to the rate function ω 2 can be obtained from, respectively, the first and second moment inequalities described in B. Denoting by N the number of solutions of a formula P SAT is the probability that N ≥ 1, and is bracketed according to (B2).

To calculate the first moment of N remark that an equation is satisfied by one half of the configurations. This result remains true for a restricted set of configurations when we average over the possible choices of (the second member of) the equation. The average number of solutions is thus 2 N /2 M , from which we get

ω 2 (α) ≥ (α -1) ln 2 . ( 24
)
This lower bound is useless for α < 1 since ω 2 is positive by definition. As for the upper bound we need to calculate the second moment N 2 of N . As equations are independently drawn

N 2 = X,Y q(X, Y ) M ( 25
)
where the sum is carried out over the pairs X, Y of configurations of the N variables, and q(X, Y ) is the probability that both X and Y satisfies the same randomly drawn equation. q can be easily expressed in terms of the Hamming distance d between X and Y , defined as the fraction of variables having opposite values in X and Y . The general expression for K-XORSAT is[71]

q(d) = 1 2 (1 -(1 -2d) K ) ( 26
)
and we specialize in this section to K = 2. Going back to [START_REF] Deroulers | [END_REF] we can sum over Y at fixed X, that is, over the distances d taking multiple values of 1 N with the appropriate binomial multiplicity, and then sum over X with the result

N 2 = 2 N d N N d q(d) M = exp(N max d∈[0;1] A(d, α) + o(N )) (27)
in the large N limit, where

A(d, α) = (2α -1) ln 2 -d ln d -(1 -d) ln(1 -d) + α ln q(d) . (28)
For α < 1 2 the maximum of A is located in d * = 1 2 , and equal to A * = 0. When α > 1 2 , A has two global maxima located in d * (α) < 1 2 and 1d * (α), with equal value A * (α) > 0. We plot in Figure 5 the lower [START_REF] Dequen | Proceedings of Theory and Applications of Satisfiability Testing, 6th International Conference[END_REF] and upper bounds to the rate function,

ω 2 (α) ≤ 2(1 -α) ln 2 -max d∈[0;1] A(d, α) (29)
from (B2). At large ratio both bounds asymptotically match, proving that ω 2 (α) = (α -1) ln 2 + O(e -2α ). As the ratio departs from its threshold value by ǫ = αα c the upper bound grows quadratically,

A * (α c + ǫ) ≃ 3 4 ǫ 2 + O(ǫ 3
). Numerics suggest that the increase of the rate function is slower,

ω 2 (α c + ǫ) ≃ Ω ǫ 3 + O(ǫ 4 ) , (30)
for some constant Ω ≃ 1 (Figure 5). We will see in Section III that a sophisticated statistical physics technique, called the replica method, actually predict this scaling with Ω = 32 27 . Actually the rate function can be estimated with the replica approach for any ratio α with the result shown in Figure 5.

G. Order parameter and symmetry breaking

What is the meaning of the Hamming distance d * (α) appearing in the calculation of the second moment of the number of solutions? An easy guess would be the average distance between pairs of solutions

d av (α) = lim N →∞ X,Y solutions of F d(X, Y ) N (F ) 2 F ( 31
)
where the average is taken over the satisfiable formulas F with ratio α, and d(X, Y ) denotes the (intensive) Hamming distance between two solutions X, Y . However an inspection of the calculation of Section II F shows that

d * (α) = lim N →∞ X,Y solutions of F d(X, Y ) F N (F ) 2 F = d av (α) . (32)
Actually, though d * (α) is not the average distance between solutions with the unbiased distribution over formulas, it is the average distance for a biased distribution where each formula is weighted with

w(F ) = N (F ) 2 F ′ N (F ′ ) 2 (33)
as can be readily checked upon insertion of w(F ) in the numerator of (31). We will see in Section (III) how to calculate average properties with the unbiased measure.

Even so definition (32) (and (31) too) is sloppy. If X is a solution so is -X, the configuration where variables values are flipped. Thus the average distance, whatever the weights over formulas, is equal 1 2 for any N ! The difficulty comes from the ambiguity in how the thermodynamic limit is taken, and is the signature of spontaneous symmetry breaking. In the low temperature phase of the Ising model the magnetization is either m * > 0 or -m * < 0 if an external field h with, respectively, positive or negative vanishing amplitude is added prior to taking the infinite size limit. In the present case what plays the role of the field is a coupling between solutions as is well-known in splin-glass theory [START_REF] Parisi | [END_REF].

Inserting exp[-N h d(X, Y )] in the numerator of (32) we obtain, when N → ∞, d * if h → 0 + and 1 -d * if h → 0 -.
The density µ of probability of distances d between solutions, with the biased measure [START_REF] Kaporis | How to prove conditional randomness using the principle of deferred decisions[END_REF], is concentrated below the Sat/Unsat threshold,

µ(d) = δ d - 1 2 for α < α c , (34)
and split into two symmetric peaks above the critical ratio,

µ(d) = 1 2 δ d -d * + 1 2 δ d -(1 -d * ) for α > α c . ( 35
)
The concept of spontaneous symmetry breaking will play a key role in our study of 3-XORSAT (Section IV C). FIG. 6: Probability that a random 3-XORSAT formula is satisfiable as a function of the ratio α of equations per variable, and for various sizes N . The dotted line locates the threshold αc ≃ 0.918.

H. Finite-size scaling (II): critical exponents

Let us summarize what we have found about the probability of satisfying random 2-XORSAT formulas in Section II E and II F. Close to the transition we have from ( 23) and (30),

ln P SAT (N, α c + ǫ) ≃ 1 4 ln(-ǫ) when ǫ < 0, N → ∞ -Ω N ǫ 3 when ǫ > 0, N ≫ 1 .
The lesson of Section II A is that ln P SAT may have a non trivial limit when N → ∞, ǫ → 0 provided we keep y = ǫ N ψ constant. For 1-XORSAT the exponent ψ was found to be equal to 1 2 , and ln P SAT to converge to the scaling function Φ 1 (y) [START_REF] Amit | Modeling Brain Function[END_REF]. The situation is similar but slightly more involved for 2-XORSAT. A natural assumption is to look for the existence of a scaling function such that

ln P SAT (N, ǫ) ≃ N ρ Φ 2 (ǫ N ψ ) . (36)
Let us see if [START_REF] Krzakala | Gibbs States and the Set of Solutions of Random Constraint Satisfaction Problems[END_REF] is compatible with the limiting behaviours [START_REF] Krzakala | Gibbs States and the Set of Solutions of Random Constraint Satisfaction Problems[END_REF]. Fixing ǫ < 0 and sending N → ∞ we obtain, for y = ǫN ψ → -∞, 1 4 ln |y| -ψ 4 ln N for the l.h.s, and N ρ × Φ 2 (y) for the r.h.s. Hence ρ = 0 as in the 1-XORSAT case, but an additive correction is necessary, and we modify scaling Ansatz [START_REF] Krzakala | Gibbs States and the Set of Solutions of Random Constraint Satisfaction Problems[END_REF]

into ln P SAT (N, ǫ) ≃ Φ 2 (y = ǫ N ψ ) - ψ 4 ln N. ( 37
)
The above equation is now compatible with (36) if Φ 2 (y) ∼ 1 4 ln |y| when y → -∞. Fixing now ǫ > 0 and sending N to infinity we see that [START_REF] Krzakala | Gibbs States and the Set of Solutions of Random Constraint Satisfaction Problems[END_REF] is fulfilled if Φ 2 (y) ∼ -Ω y 3 when y → +∞ and

ψ = 1 3 . (38)
The above value for ψ is expected from the study of random graphs [START_REF] Bollobas | Random Graphs[END_REF] and is related to the size N 1-ψ = N 2 3 of the largest components at the percolation threshold (Section II D). ψ is called critical exponent and characterize the width of the critical region of 2-XORSAT. Loosely speaking it means that a formula of with N variables and N 2 + ∆ equations is 'critical' when ∆ ∼ N 2 3 . This information will be useful for the analysis of search algorithms in Section V F.

A consequence of [START_REF] Krzakala | A landscape analysis of constraint satisfaction problems[END_REF][START_REF] Mertens | [END_REF] is that, right at the threshold, the probability of satisfaction decays as[72]

P SAT (N, α c ) ∼ N -1 12 . ( 39
)
This scaling agrees with numerical experiments, though the small value of the decay exponent makes an accurate check delicate (Inset of Figure 5).

I. First and second moments inequalities for the 3-XORSAT threshold

Figure 6 shows the probability that a random 3-XORSAT formula is satisfiable as a function of α for increasing sizes N . It appears that formulas with ratio α < α c ≃ 0.92 are very likely to be satisfiable in the large N limit, while formulas with ratios beyond this critical value are almost surely unsatisfiable. This behaviour is different from the 2-XORSAT case (Figure 4) in that P SAT seems to tend to unity below threshold.

It is important to realize that, contrary to the 2-XORSAT case, the Sat/Unsat transition is not related to connectivity percolation. Consider indeed a variable, say, x 1 . This variable appear, on average, in 3α equations. Each of those equations contain other 2 variables. Hence the 'connectivity' of x 1 is c = 6α, which is larger than unity for α p = 1 6 . In the range [α p , α c ] the formula is percolating but still satisfiable with high probability. The reason is that cycles do not hinder satisfiability as much as in the 2-XORSAT case.

Use of the first and second moment inequalities (B) for the number N of solutions provides us with upper and lower bounds to the Sat/Unsat ratio α c . The calculation follows the same line as the one of the 2-XORSAT case (Section II F). The first moment N = 2 N (1-α) vanishes for ratios larger than unity, showing that

α c ≤ α 1 = 1 . ( 40
)
This upper bound is definitely larger than the true threshold from the numerical findings of Figure 6. We have already encountered this situation in 2-XORSAT: in the 1 2 < α < 1 range formulas are unsatisfiable with probability one (when N → ∞), yet the average number of solutions is exponentially large! The reason is, once more, that the average result is spoiled by rare, satisfiable formulas with many solutions.

As for the second moment expression [START_REF] Engel | Statistical Mechanics of Learning[END_REF][START_REF] Feige | Proceedings of 34th STOC conference[END_REF] still holds with q(d) given by ( 26) with K = 3. The absolute maximum of the corresponding function A(d, α) is located in d * = 1 2 when α < α 2 ≃ 0.889, and d * < 1 2 when α > α 2 . In the latter case N 2 is exponentially larger than N 2 , and the second moment inequality (B2) does not give any information about P SAT . In the former case N 2 and N 2 are equivalent to exponential-in-N order. It is shown in C that their ratio actually tends to one as N → ∞. We conclude that formulas with ratios of equations per variable less than α 2 are satisfiable with high probability in the infinite size limit, or, equivalently [20],

α c ≥ α 2 ≃ 0.889 . ( 41
)
Unfortunately the lower and upper bounds do not match and the precise value of the threshold remains unknown at this stage. We explain in the next section how a simple preprocessing of the formula, before the application of the first and second moment inequalities, can close the gap, and shed light on the structure of the space of solutions.

We start from a simple observation. Assume we have a formula F of 3-XORSAT where a variable, say, x, appears only once, that is, in one equation, say, E : x + y + z = 0. Let us call F ′ the subformula obtained from F after removal of equation E. Then the following statement is true: F is satisfiable if and only if F ′ is satisfiable. The proof is obvious: whatever the values of y, z required to satisfy F ′ equation E can be satisfied by an adequate choice of x, and so can be the whole formula F .

In a random 3-XORSAT formula F with ratio α there are about N × 3α e -3α variables appearing only once in the formula. Removal of those variables (and their equations) produces a shorter formula with O(N ) less equations. Furthermore it may happen that variables with multiple occurrences in the original formula have disappeared from the output formula, or appear only once. Hence the procedure can be iterated until no single-occurrence variables are present. We are left with F 2 , the largest subformula (of the original formula) where every variable appears at least twice.

Many questions can be asked: how many equations are left in F 2 ? how many variables does it involve? how many solutions does it have? Giving the answers requires a thorough analysis of the removal procedure, with the techniques exposed in Section V E [17,[START_REF] Dubois | Proc. of the 43rd annual IEEE symposium on Foundations of Computer Science[END_REF]46]. The outcome depends on the value of the ratio compared to

α d = min b - log(1 -b) 3 b 2 ≃ 0.8184 . . . (42)
hereafter called clustering threshold. With high probability when N → ∞ F 2 is empty if α < α d , and contains an extensive number of equations, variables when α > α d . In the latter case calculation of the first and second moments of the number of solutions of F 2 shows that this number does not fluctuate around the value e N s cluster (α)+o(N ) where

s cluster (α) = (b -3α b 2 + 2α b 3 ) ln 2 ( 43
)
and b is the strictly positive solution of the self-consistent equation Hence F 2 is satisfiable if and only if α < α c defined through s cluster (α c ) = 0, that is,

1 -b = e -3 α b 2 . ( 44
)
α c ≃ 0.9179 . . . . ( 45
)
This value is, by virtue of the equivalence between F and F 2 the Sat/Unsat threshold for 3-XORSAT, in excellent agreement with Figure 6. How can we reconstruct the solutions of F from the ones of F 2 ? The procedure is simple. Start from one solution of F 2 (empty string if α < α d ). Then introduce back the last equation which was removed since it contained n ≥ 1 single-occurrence variable. If n = 1 we fix the value of this variable in a unique way. If n = 2 (respectively n = 3) there are 2 (respectively, 4) ways of assigning the reintroduced variables, defining as many solutions from our initial, partial solution. Reintroduction of equations one after the other according to the Last In -First Out order gives us more and more solutions from the initial one, until we get a bunch of solutions of the original formula F . It turns out that the number of solutions created this way is e N sin(α)+o (N ) where

s in (α) = (1 -α) ln 2 -s cluster (α) . (46)
The above formula is true for α > α d , and should be intended as s in (α) = (1α) ln 2 for α < α d . These two entropies are shown in Figure 7. The total entropy, s * (α) = s in (α) + s cluster (α), is simply (1α) ln 2 for all ratios smaller than the Sat/Unsat threshold. It shows no singularity at the clustering threshold. However a drastic change in the structure of the space of solutions takes place, symbolized in the phase diagram of Figure 8:

• For ratios α < α d the intensive Hamming distance between two solutions is, with high probability, equal to d = 1/2. Solutions thus differ on N/2 + o(N ) variables, as if they were statistically unrelated assignments of the N Boolean variables. In addition the space of solutions enjoys some connectedness property. Any two solutions are connected by a path (in the space of solutions) along which successive solutions differ by a bounded number of variables. Losely speaking one is not forced to cross a big region prived of solutions when going from one solution to another.

• For ratios α > α d the space of solutions is not connected any longer. It is made of an exponentially large (in N ) number N clu = e N s cluster of connected components, called clusters, each containing an exponentially large number N in = e N sin of solutions. Two solutions belonging to different clusters lie apart at a Hamming distance d clu = 1/2 while, inside a cluster, the distance is d in < d clu . b given by ( 44) is the fraction of variables having the same value in all the solutions of a cluster (defined as the backbone).

We present in Sections III and IV statistical physics tools developed to deal with the scenario of Figure 8.

III. ADVANCED METHODS (I): REPLICAS

A. From moments to large deviations for the entropy

The analysis of Section II F has shown that the first, and second moments of the number N of solutions are dominated by rare formulas with a lot of solutions. Let us define the intensive entropy s through N = e N s . As N is random (at fixed α, N ) so is s. We assume that the distribution of s can be described, in the large size limit, by a rate function ω(s) (which depends on α). Hence,

N q = ds e -N ω(s) × e N s q ∼ exp N max s q s -ω(s) (47)
using the Laplace method. If we are able to estimate the leading behaviour of the q th moment of the number of solutions when N gets large at fixed α,

N q ∼ e N g(q) , (48)
then ω can be easily calculated by taking the Legendre transform of g. In particular the typical entropy is obtained by s * = dg dq (q → 0). This is the road we will follow below. We will show how g(q) can be calculated when q takes integer values, and then perform an analytic continuation to non integer q. The continuation leads to substantial mathematical difficulties, but is not uncommon in statistical physics e.g. the q → 1 limit of the q-state Potts model to recover percolation, or the n → 0 limit of the O(n) model to describe self-avoiding walks.

To calculate the q th moment we will have to average over the random components of formulas F , that is, the Kuplets of index variables in the first members and the v = 0, 1 second members. Consider now homogeneous formulas F h whose first members are randomly drawn in the same way as for F , but with all second members v = 0. The number N h of solutions of a homogeneous formula is always larger or equal to one. It is a simple exercise to show that

N q+1 = 2 N (1-α) × N h q , ( 49
)
valid for any positive integer q[73]. Therefore it is sufficient to calculate the moments of N h = e N g h (q) since ( 49) gives a simple identity between g(q + 1) and g h (q). This technical simplification has a deep physical meaning we will comment in Section IV C.

B. Free energy for replicated variables

The q th power of the number of solutions to a homogeneous system reads

N h q = X M ℓ=1 e ℓ (X) q = X 1 ,X 2 ,...,X q M ℓ=1 q a=1 e ℓ (X a ) , (50)
where e ℓ (X) is 1 if equation ℓ is satisfied by assignment X. The last sum runs over q assignments X a , with a = 1, 2, . . . , q of the Boolean variables, called replicas of the original assignment X. It will turn useful to denote by x i = (x 1 i , x 2 i , . . . , x q i ) the q-dimensional vector whose components are the values of variable x i in the q replicas. To simplify notations we consider the case K = 3 only here, but extension to other values of K is straightforward. Averaging over the instance, that is, the triplets of integers labelling the variables involved in each equation ℓ, leads to the following expression for the q th moment,

N h q = X 1 ,X 2 ,...,X q q a=1 e(X a ) M = X 1 ,X 2 ,...,X q   1 N 3 1≤i,j,k≤N δ xi+ xj + x k + O 1 N   M (51)
where δ x = 1 if the compoments of x are all null mod. 2, and 0 otherwise. We now procede to some formal manipulations of the above equation ( 51). First step. Be X = {X 1 , X 2 , . . . , X q } one of the 2 qN replica assignment. Focus on variable i, and its attached assignment vector, x i . The latter may be any of the 2 q possible vectors e.g. x i = (1, 0, 1, 0, 0, . . . , 0) if variable x i is equal to 0 in all but the first and third replicas. The histogram of the assignments vectors given replica assignment X ,

ρ x|X = 1 N N i=1 δ x-xi , (52)
counts the fraction of assignments vectors x i having value x when i scans the whole set of variables from 1 to N . Of course, this histogram is normalised to unity,

x ρ x = 1 , (53)
where the sum runs over all 2 q assignment vectors. An simple but essential observation is that the r.h.s. of ( 51) may be rewritten in terms of the above histogram,

1 N 3 1≤i,j,k≤N δ xi+ xj+ x k = x, x ′ ρ x ρ x ′ ρ x + x ′ . (54)
Keep in mind that ρ in (52,54) depends on the replica assignement X under consideration.

Second step. According to [START_REF] Motwani | Randomized algorithms[END_REF], two replica assignments X 1 and X 2 defining the same histogram ρ will give equal contributions to N h q . The sum over replica assignments X can therefore be replaced over the sum over possible histograms provided the multiplicity M of the latter is taken properly into account. This multiplicity is also equal to the number of combinations of N elements (the x i vectors) into 2 q sets labelled by x and of cardinalities N ρ( x). We obtain

N h q = (norm) {ρ} e N G h {ρ}, α + o(N ) , (55)
where the (norm) subscript indicates that the sum runs over histograms ρ normalized according to [START_REF] Montanari | [END_REF], and

G h {ρ}, α = - x ρ(x) ln ρ(x) + α ln x, x ′ ρ x ρ x ′ ρ x + x ′ . ( 56
)
In the large N limit, the sum in ( 55) is dominated by the histogram ρ * maximizing the functional G h . Third step. Maximisation of function G h over normalized histograms can be done within the Lagrange multiplier formalism. The procedure consists in considering the modified function

G LM h {ρ}, λ, α = G h {ρ}, α + λ 1 - x ρ x , (57)
and first maximise G LM h with respect to histograms ρ without caring about the normalisation constraint, and then optimise the result with respect to λ. We follow this procedure with G h given by (56). Requiring that G LM h be maximal provides us with a set of 2 q coupled equations for ρ * , ln ρ * ( x)

+ 1 + λ -3 α x ′ ρ * x ′ ρ * x + x ′ x ′ , x ′′ ρ * x ′ ρ * x ′′ ρ * x ′ + x ′′ = 0 , (58)
one for each assignment vector x. The optimisation equation over λ implies that λ in ( 58) is such that ρ * is normalised. At this point of the above and rather abstract calculation it may help to understand the interpretation of the optimal histogram ρ * .

C. The order parameter

We have already addressed a similar question at the end of the second moment calculation in Section II G. The parameter d * coming out from the calculation was the (weighted) average Hamming distance [START_REF] Impagliazzo | Proceedings of the IEEE Conference on Computational Complexity[END_REF] between two solutions of the same random instance. The significance of ρ * is identical. Consider q ′ solutions labelled by a = 1, 2, . . . , q ′ of the same random and homogeneous instance and a variable, say, x i . What is the probability, over instances and solutions, that this variable takes, for instance, value 0 in the first and fourth solutions, and 1 in all other solutions? In other words, what is the probability that the assignment vector

x i = (x 1 i , x 2 i , . . . , x q ′ i ) is equal to x ′ = (0, 1, 1, 0, 1, 1, . . . , 1)? The answer is p( x ′ ) = 1 (N h ) q ′ X 1 ,X 2 ,...,X q ′ δ xi-x M l=1 q a=1 e ℓ (X a ) ( 59
)
where the dependence on i is wiped out by the average over the instance. The above probability is an interesting quantity; it provides us information about the 'microscopic' nature of solutions. Setting q ′ = 1 gives us the probabilities p(0), p(1) that a variable is false or true respectively, that is, takes the same value as in the null assignment or not.

For generic q ′ we may think of two extreme situations:

• a flat p over assignment vectors, p( x ′ ) = 1/2 q ′ , corresponds to essentially orthogonal solutions;

• on the opposite, a concentrated probability e.g. p( x ′ ) = δ x ′ implies that variables are extremely constrained, and that the (almost) unique solution is the null assignment.

The careful reader will have already guessed that our calculation of the q th moment gives access to a weighted counterpart of p. The order parameter

ρ * ( x) = 1 (N h ) q X 1 ,X 2 ,...,X q δ xi-x M l=1 q a=1 e ℓ (X a ) , (60)
is not equal to p even when q = q ′ . However, at the price of mathematical rigor, the exact probability p over vector assignments of integer length q ′ can be reconstructed from the optimal histogram ρ * associated to moments of order q when q is real-valued and sent to 0. The underlying idea is the following. Consider (60) and an integer q ′ < q.

From any assignment vector x of length q, we define two assignment vectors x ′ , x ′′ of respective lengths q ′ , qq ′ corresponding to the first q ′ and the last qq ′ components of x respectively. Summing (60) over the 2 q-q ′ assignment vectors x ′′ gives,

x ′′ ρ * ( x ′ , x ′′ ) = 1 (N h ) q {X a } δ x ′ i -x ′ N h q-q ′ l,a e ℓ (X a ) . ( 61
)
As q now appears in the powers of N h in the numerator and denominator only, it can be formally send to zero at fixed q ′ , yielding lim

q→0 x ′′ ρ * ( x ′ , x ′′ ) = p( x ′ ) (62)
from [START_REF] Parisi | [END_REF]. This identity justifies the denomination order parameter given to ρ * .

Having understood the significance of ρ * helps us to find appropriate solutions to [START_REF] Papadimitriou | Proceedings of the 32nd Annual IEEE Symposium on Foundations of Computer Science[END_REF]. Intuitively and from the discussion of the first moment case q = 1, p is expected to reflect both the special role of the null assignment (which is a solution to all homogeneous systems) and the ability of other solutions of a random system to be essentially orthogonal to this special assignment. A possible guess is thus

p( x ′ ) = 1 -b 2 q ′ + b δ x ′ , (63)
where b expresses some degree of 'correlation' of solutions with the null one. Hypothesis (63) interpolates between the fully concentrated (b = 1) and flat (b = 0) probabilities. b measures the fraction of variables (among the N ones) that take the 0 values in all q ′ solution, and coincides with the notion of backbone introduced in Section II J. Hypothesis ( 63) is equivalent, from the connection (62) between p and the annealed histogram ρ * to the following guess for the solution of the maximisation condition [START_REF] Papadimitriou | Proceedings of the 32nd Annual IEEE Symposium on Foundations of Computer Science[END_REF],

ρ * ( x) = 1 -b 2 q + b δ x . (64)
Insertion of Ansatz ( 64) in [START_REF] Papadimitriou | Proceedings of the 32nd Annual IEEE Symposium on Foundations of Computer Science[END_REF] shows that it is indeed a solution provided b is shrewdly chosen as a function of q and α, b = b * (q, α). Its value can be either found from direct resolution of ( 58), or from insertion of histogram (64) in G h (56) and maximisation over b, with the result,

g h (q, α) = max 0≤b≤1 A h (b, q, α) (65)
where

A h (b, q, α) = -1 - 1 2 q (1 -b) ln 1 -b 2 q (66) -b + 1 -b 2 q ln b + 1 -b 2 q + α ln b 3 + 1 -b 3 2 q ,
where the maximum is precisely reached in b * . Notice that, since ρ * in ( 64) is entirely known from the value of b * , we shall indifferently call order parameter ρ * , or b * itself.

D. Results

Numerical investigation of A h (66) shows that: for α < α M (q) the only local maximum of A h is located in b * = 0, and A h (q, α) = q(1α) ln 2; when α M (q) < α < α * (q), there exists another local maximum in b > 0 but the global maximum is still reached in b * = 0; when α > α * (q), the global maximum is located in b * > 0. This scenario extends to generic q the findings of the second moment calculation carried out in Section II F. The α M and α * lines divide the q, α plane as shown in Figure 9. Notice that, while the black dots in Figure 9 correspond to integer-valued q, the continuous lines are the output of the implicit analytic continuation to real q done by the replica calculation.

Taking the derivative of (65) with respect to q and sending q → 0 we obtain the typical entropy of a homogeneous 3-XORSAT system at ratio α,

s * h (α) = ln 2 × max 0≤b≤1 (1 -b) 1 -ln(1 -b) -α (1 -b 3 ) . ( 67
)
The optimal value for b coincides with the solution of [START_REF] Mézard | [END_REF]. The typical entropy is plotted in Figure 10, and is equal to:

• (1α) ln 2 when α < α c ≃ 0.918 (Figure 9); in this range of ratios, homogeneous and full (with random second members) systems have essentially the same properties, with the same cluster organisation of solutions, and identical entropies of solutions.

• a positive but rapidly decreasing function given by (67) when α > α c ; above the critical ratio, a full system has no solution any more, while a homogeneous instance still enjoys a positive entropy. The expression for s * h (α) coincides with the continuation to α > α c of the entropy s in (α) (46) of solutions in a single cluster for a full system. In other words, a single cluster of solutions, the one with the null solution, survive for ratios α > α S in homogeneous systems. FIG. 9: The q, α plane and the critical lines αM (q) (dashed), α * (q) (full tick), and αs(q) (full thin) appearing in the calculation of the q th moment for homogeneous 3-XORSAT systems. Full dots correspond to integer q values, while continuous curves result from the analytic continuation to real q. The fraction of variables in the backbone, b * , vanishes below the line αM (q); the global maximum of A h in (66) is located in b * > 0 for ratios α > α * (q). Ansatz ( 64) is locally unstable in the hardly visible domain q < 0, αM (q) < α < αs(q).

Atypical instances can be studied and the large deviation rate function for the entropy can be derived from (65) for homogeneous systems, and using equivalence (49), for full systems. Minimizing over the entropy we obtain the rate function ω 3 (α) associated to the probability that a random 3-XORSAT system is satisfiable, with the result shown in Figure 10. As expected we find ω 3 = 0 for α < α c and ω 3 > 0 for α > α c , allowing us to locate the Sat/Unsat threshold.

Notice that the emergence of clustering can be guessed from Figure 9. It coincides with the appearance of a local maximum of A h (66) with a non vanishing backbone b. While in the intermediate phase α d < α < α c , the height of the global maximum equals the total entropy s * , the height of the local maximum coincides with the entropy of clusters s cluster [START_REF] Mézard | Spin glasses and beyond[END_REF].

E. Stability of the replica Ansatz

The above results rely on Ansatz (64). A necessary criterion for its validity is that ρ * locates a true local maximum of G h , and not merely a saddle-point. Hence we have to calculate the Hessian matrix of G h in ρ * , and check that the eigenvalues are all negative [START_REF] De Almeida | [END_REF]. Differentiating (56) with respect to ρ( x) and ρ( x ′ ) we obtain the Hessian matrix

H( x, x ′ ) = - δ x+ x ′ ρ * ( x) + 6α ρ * ( x + x ′ ) D -9α N ( x) D N ( x ′ ) D , (68)
where

D = 1-b 3 2 q + b 3 , N ( x) = 1-b 2 2 q + b 2 δ x .
We use b instead of b * to ligthen the notations, but it is intended that b is the backbone value which maximizes A h (66) at fixed q, α. To take into account the global constraint over the histogram [START_REF] Montanari | [END_REF] one can express one fraction, say, ρ( 0), as a function of the other fractions ρ( x), x = 0. G H is now a fonction of 2 q -1 independent variables, with a Hessian matrix H simply related to H,

H( x, x ′ ) = H( x, x ′ ) -H( x, 0) -H( 0, x ′ ) + H( 0, 0) . (69)
Plugging expression (68) into (69) we obtain

H( x, x ′ ) = λ R δ x+ x ′ + 1 2 q -1 λ L -λ R ) where λ R = 6α b D - 2 q 1 -b (70) λ L = 2 q 6α b D - 2 q (1 -b)(1 -b + 2 q b) -9α(1 -2 -q ) b 4 D 2 .
Diagonalization of H is immediate, and we find two eigenvalues:

• λ L (non degenerate). The eigenmode corresponds to a uniform infinitesimal variation of ρ( x) for all x = 0, that is, a change of b in (64). It is an easy check that where A h is defined in (66). As we have chosen b to maximize A h this mode, called longitudinal in replica literature [START_REF] De Almeida | [END_REF], is stable[74].

λ L = 2 q 1 -2 -q ∂ 2 A h ∂b 2 (b, q, α) , (71)
• λ R (2 q -2-fold degenerate): the eigenmodes correspond to fluctuations of the order parameter ρ transverse to the replica subspace described by (64), and are called replicon in spin-glass theory [START_REF] De Almeida | [END_REF]. Inspection of λ R as a function of α, q shows that it is always negative when q > 0. For q < 0 the replicon mode is stable if

α > α s (q) = 1 -b 3 + 2 q b 3 6 b(1 -b) . ( 72
)
which is a function of q only once we have chosen b = b * (q, α s ).

The unstable region q < 0, α M (q) < α < α s (q) is shown in Figure 9 and is hardly visible when q > -3. In this region a continuous symmetry breaking is expected [START_REF] Mézard | Spin glasses and beyond[END_REF]. In particular α s stay below the α * line for small (in absolute value) and negative q. We conclude that our Ansatz (64) defines a maximum of G h . Is it the global maximum of G h ? There is no simple way to answer this question. Local stability does not rule out the possibility for a discontinuous transition to another maximum in the replica order parameter space not described by (64). A final remark is that a similar calculation can be done for any value of K. The outcome for K = 2 is the rate function ω 2 plotted in Figure 5, in good agreement with numerics close to the threshold.

IV. ADVANCED METHODS (II): CAVITY

The cavity method, in the context of disordered systems, was historically developed as an alternative to the the replica method [START_REF] Mézard | Spin glasses and beyond[END_REF]. Its application to spin systems on random graphs is extensively explained in [START_REF] Mézard | [END_REF], and we limit ourselves here to briefly show how it gives back the 3-XORSAT scenario of Section II J [46].

Let us consider a system F involving variables x i , i = 1, . . . , N . In the following we will indifferently use the variable x i = 0, 1 or its spin representation S i = (-1) xi = ±1 when convenient. Let us define the GS energy E F i (S i ) of the system when the i th spin is kept fixed, that is, the minimal number of violated equations in F , taken over the 2 N -1 configurations. We may always write

E F i (S i ) = -w i -h i S i , (73)
where h i is called 'field' acting on spin S i . For a homogeneous system E F i (+1) = 0, and E F i (-1) = n i for some integer n i . Hence h i = ni 2 takes half-integer values. The above definition can be extended to the case of ℓ > 1 fixed spins. Let I ⊂ {1, 2, . . . , N } be a subset of the indices of cardinal |I| ≥ 2, and S I denote one of the 2 |I| configurations of the spins S i , i ∈ I. The GS energy of F for given S I can in general be written as

E F I (S I ) = -w I - i∈I h i S i - I ′ ⊂I:|I ′ |≥2 J I ′ i∈I ′ S i (74)
where the h i s are the fields and the J I ′ s are effective couplings between subsets of spins.

The basic cavity assumption is that effective couplings are vanishingly small: J I ′ = 0 for every subset I ′ . This apparently bold hypothesis critically relies on a general property of random graphs (from which our system is built on). Define the distance between two vertices as the minimal number of edges on pathes linking these two points. Then vertices in a finite subset are, with high probability when N → ∞, typically at infinite distance from each other [75]. When correlations between variables in GS extinguish with the distance i.e. when the correlation length is finite the cavity assumption is correct in the large N limit [START_REF] Mézard | Spin glasses and beyond[END_REF][START_REF] Montanari | [END_REF]. The assumption will break down when correlations subsist on infinite distance, which happens to be the case in the clustered phase.

A. Self-consistent equation for the fields

Under the assumption that couplings between randomly picked up spins are null we are left with the fields only. The goal of this section is to show how to calculate those fields, or more precisely, their probability distribution. The derivation is based on the addition procedure already used in the calculation of the size of the giant component in random graphs (Section II D).

Consider a system F over N variables to which we want to add one equation involving one new variable S, and two variables S 1 , S 2 appearing in F . The energy function associated to this equation is

e(S, S 1 , S 2 ) = 1 2 1 -σ S S 1 S 2 (75)
where σ = +1, respectively -1, when the second member of the equation is 0, resp. 1. Let us calculate the GS energy of the new system F ′ = F + added equation when the new variable S is kept fixed,

E F ′ (S) = min S1,S2 e(S, S 1 , S 2 ) + E F 1,2 (S 1 , S 2 ) = -w -u S . (76)
With the cavity hypothesis the couplings between spins S 1 , S 2 is null and the minimization is straightforward. We deduce the following explicit expression for the field acting on S (called bias in the cavity literature [START_REF] Mézard | [END_REF]),

u = σ 2 sign h 1 h 2 . (77)
Suppose we now add ℓ ≥ 1 (and not only one) equations. The above calculation can be easily repeated. The absence of couplings make the total field acting on S a linear combination of the fields coming from each new equation,

h = ℓ j=1 u j , (78)
where u j is calculated from (77) and each pair of fields (h j 1 , h j 2 ) acting on the spins in the j th equation, j = 1, . . . , ℓ. How many equations should we add for our new system over N + 1 variables to have the same statistical features as old one over N variables? First ℓ should be Poisson distributed with parameter 3α. Then, given ℓ, we randomly chose ℓ pairs of variables; for each pair the corresponding bias u can be calculated from (77). Assume the output is a set of ℓ independent biases, taking values

u =    + 1
2 with probability a + 0 with probability a 0 -1 2 with probability a -

Obviously a + + a 0 + a -= 1. Summing over the equations as in (78) we obtain the distribution of the field h acting on the new spin at fixed ℓ,

p(h|ℓ) = ℓ+,ℓ0,ℓ- ℓ ℓ+, ℓ 0 , ℓ - a ℓ+ + a ℓ0 0 a ℓ- -δ h-1 2 (ℓ+-ℓ-) . (80)
Finally we sum over the Poisson distribution for ℓ to obtain the distribution of fields h,

p(h) = e -3α(1-a0) ℓ+,ℓ- (3α) ℓ++ℓ- ℓ + !ℓ -! a ℓ+ + a ℓ- -δ h-1 2 (ℓ+-ℓ-) . (81)
In turn we calculate the distribution of the biases from the one of the fields through (77). The outcome are the values of the probabilities (79) in terms of p,

a + = h1,h2:h1h2>0 p(h 1 ) p(h 2 ) , a -= h1,h2:h1h2<0 p(h 1 ) p(h 2 ) , a 0 = h1,h2:h1h2=0 p(h 1 ) p(h 2 ) = 2 p(0) -p(0) 2 . ( 82
)
The above equations together with (81) define three self-consistent conditions for a 0 , a + , a -. Notice that the free energy can be calculated along the same lines [START_REF] Mézard | [END_REF].

B. Application to homogeneous and full systems

In the case of homogeneous systems (σ = +1) we expect all the fields to be positive, and look for a solution of (82) with a -= 0. Then p(h) (81) is a Poisson distribution for the integer-valued variable 2h, with parameter 3 α a + . The self-consistent equation (82) reads

a 0 = 1 -a + = 2 e -3αa+ -e -6 αa+ = 1 -1 -e -3 αa+ 2 (83)
which coincides with [START_REF] Mézard | [END_REF] with the definition b = √ a + . As expected

b = h≥ 1 2 p(h) (84)
is the fraction of frozen variables (which cannot be flipped from 0 to 1 in GS assignments), in agreement with the notion of backbone of Section II J. The energy is zero at all ratio α by construction. As for the entropy consider adding a new equation to the system F (but with no new variable). With probability 1b 3 at least one of the three variables in the new equation e was not frozen prior to addition, and the number of solutions of the new system F + e is half the one of F . With probability b 3 all three variables are frozen in F (to the zero value) and the number of solutions of F + e is the same as the one of F . Hence the average decrease in entropy is

N s * h (α + 1 N ) -N s * h (α) ≃ ds * h dα = -(1 -b 3 ) ln 2 . ( 85
)
The same differential equation can be obtained by differentiating (67). With the limit condition s * h (α → ∞) = 0 we obtain back the correct expression for the average entropy of homogeneous systems. The entropy is equal to (1-α) ln 2 at α = α c , and becomes smaller when the ratio decreases. This shows that the solution b = 0 must be preferred in this regime to the metastable b > 0 solution. We conclude that the cavity assumption leads to sensible results for homogeneous systems at all ratios α.

In full systems the sign σ entering (77) takes ±1 values with equal probabilities. We thus expect p(h) to be an even distribution, and a + = a -= 1 2 (1a 0 ). Remark that a solution with a 0 < 1 cannot exist in the satifiable phase. It would allow two added equations to impose opposite non zero biases to the new variable i.e. to constraint this variable to take opposite values at the same time. Given a 0 we calculate from (81) the probability that the field vanishes,

p(0) = e -3 α(1-a0) ∞ ℓ=0 3α 2 (1 -a 0 ) 2ℓ 1 ℓ! 2 (86)
and, in turn, derive from (82) a self-consistent equation for a 0 . Numerical investigations show that a 0 = 1 is the unique solution for α < α T = 1.167. When α > α T there appears another solution with a 0 < 1. The clustering and Sat/Unsat transitions are totally absent. This result, incompatible with the exact picture of random 3-XORSAT exposed in Section II J, shows that the simple cavity hypothesis does not hold for full systems.

C. Spontaneous symmetry breaking between clusters

In the clustered phase variables are known to be strongly correlated and the cavity assumption has to be modified. Actually from what we have done above in the homogeneous case we guess that the independence condition still holds if we can in some way restrict the whole space of solutions to one cluster. To do so we explicitely break the symmetry between clusters as follows [START_REF] Monasson | [END_REF][START_REF] Parisi | [END_REF].

Let S * i , i = 1, . . . , N be a reference solution of a full satisfiable system F , and F h the corresponding homogeneous system. We define the local gauge transform S i → Ŝi = S i × S * i . {S} is a solution of F if and only if { Ŝ} is a solution of F h . As the cavity assumption is correct for the homogeneous system we obtain the distribution of fields ĥi ≥ 0 from (81). Gauging back to the original spin configuration gives us the fields

h i = S * i × ĥi . ( 87
)
It turns out that the above fields depend only on the cluster to which belong the reference solution. Indeed for the fraction 1b of the non frozen spins, ĥi = h i = 0. For the remaining fraction b of spins in the backbone ĥi > 1 2 and S * i has a unique value for all solutions in the cluster (Section II J). Hence the fields h i are a function of cluster (c) containing {S * }, and will be denoted by h i . How do we perform this restriction in practice? A natural procedure is to break the symmetry between clusters in an explicit manner by adding a small coupling to the reference solution [START_REF] Motwani | Randomized algorithms[END_REF][START_REF] Parisi | [END_REF]. Remark that symmetry was broken (naturally but explicitly!) in the case of homogeneous systems when we looked for a distribution p(h) with support on positive fields only. It is a remarkable feature of XORSAT (rather unique among disordered systems) that symmetry between disordered clusters can be broken in a constructive and simple way.

The main outcome of the above discussion is that the field attached to variable i is not unique, but depends on the cluster (c). We define the distribution p i (h) of the fields attached to variable i over the clusters (with uniform weights since all clusters contain the same number of solutions) [50]. The naive cavity assumption corresponds to

p i (h) = δ h-hi . ( 88
)
In presence of many clusters p i (h) is not highly concentred. From (87) and the fact that S * i = ±1 depending on the cluster from which we pick up the reference solution we find that

p i (h) = 1 2 δ h-ĥi + δ h+ ĥi . ( 89
)
As ĥi is itself randomly distributed we are led to introduce the distribution P of the field distributions p i (h). This mathematical object, P(p(h)), is the order parameter of the cavity theory in the clustered phase [START_REF] Mézard | [END_REF]50].

D. Distribution of field distributions

Let us see how P can be obtained within the one-more variable approach of Section IV A. A new equation contains two variables S i , S j from F , with fields h

(c) i , h (c) j
in each cluster (c). The bias u is a deterministic function of those two fields for each cluster (77). We define its distribution over clusters ρ. As u can take three values only and ρ is an even distribution due to the randomness of the second member of the new equation we may write

ρ(u) = (1 -ϕ ij ) δ u + ϕ ij 2 δ u-1 2 + δ u+ 1 2 . ( 90
)
The weight ϕ is a random variable which varies from pair (ij) to pair. What is the probability distribution p(ϕ) of ϕ? Either the two variables in the pair belong to the backbone and they are frozen in all clusters; then u will be non zero and ϕ = 1. Or one (at least) of the two variables is not frozen and u = 0 in all clusters, giving ϕ = 0. We may write

p(ϕ) = (1 -w) δ ϕ + w δ ϕ-1 . ( 91
)
From the above argument we expect w = b 2 . Let us derive this result. Assume we add ℓ ≥ 1 equations to our system. For each one of those equations a bias u j is drawn randomly according to distribution (90). Denote by m(≤ ℓ) the number of those equations with parameter ϕ = 1; m is binomially distributed with probability w among ℓ. Then ℓm biases are null, and m biases are not equal to zero.

For the formula to remain satisfiable the non-zero biases must be all positive or negative [START_REF] Mézard | [END_REF], see Section IV B. Hence the distribution of the field on the new variable is

p m (h) = 1 2 δ h-m 2 + δ h+ m 2 , ( 92
)
in agreement with the expected form (89). The upperscript m underlines that field distributions with non zero probability are can be labelled by an integer m; they define a countable set and the distribution P can be defined as a discrete probability P m over the set of positive integers m. The probability P m of distribution (92) is the convolution of binomial distribution for m at fixed ℓ with the Poisson distribution over ℓ,

P m = e -3 α w (3 α w) m m! . ( 93
)
Identities (92,93) fully determine the distribution of field distributions in term of a single parameter, w.

To close the self-consistency argument consider the two variables in F in, say, the first added equation. Call h, h ′ their fields, distributed according to p m (h), p m ′ (h ′ ) for some m, m ′ . The bias created onto the new variable will be non zero if h and h ′ may both take non zeros value in some clusters, that is, if m and m ′ are not equal to zero. This translates into the mathematical identity

w = m≥ 1 2 ,m ′ ≥ 1 2 P m P m ′ = 1 -P 0 2 = 1 -e -3 α w 2 (94)
from (93). The above equation coincides with [START_REF] Mézard | [END_REF] for w = b 2 . Notice that w is equal to the probability a + that the bias is non zero in the homogeneous case (79), in agreement with the discussion of Section IV C. It is easy to find back the expressions for the entropies of clusters, s cluster , and solutions in a cluster, s in , given in Section II J. As for the latter entropy the argument leading to (85) can be repeated, with the modification that the second member of the added equation is not necessarily zero but the value it should have for the equation to be satisfied when all three variables are frozen. Hence (85) holds with s * h replaced with s in . As for the entropy of clusters the same argument again tells us that, on average, half of the clusters will disappear when the three variables are frozen and the second member of the equation is randomly chosen. Therefore

ds cluster dα = -b 3 ln 2 , (95)
in agreement with equations [START_REF] Mézard | Spin glasses and beyond[END_REF][START_REF] Mézard | [END_REF]. Summing differential equations ( 95) and (85) for s cluster and s in respectively shows that the total entropy of solutions is (1α) ln 2 (Section II J).

V. DYNAMICAL PHASE TRANSITIONS AND SEARCH ALGORITHMS

The deep understanding of the statistical properties of 3-XORSAT makes this problem a valuable benchmark for assessing the performances of various combinatorial search algorithms. At first sight, the idea seems rather odd since 3-XORSAT is a polynomial problem. Interestingly most of the search procedures devised to deal with NP-complete problems e.g. SAT have poor performances i.e. take exponentially long average running times on XORSAT above some algorithmic-dependent critical ratio ... The purpose of this Section is to present two algorithms exhibiting such a dynamical phase transition, and the techniques required for their analysis.

A. Random WalkSAT (RWSAT): definition, worst-case bound

The first algorithm we consider is the Random WalkSAT (RWSAT) algorithm introduced by Papadimitriou [START_REF] Papadimitriou | Proceedings of the 32nd Annual IEEE Symposium on Foundations of Computer Science[END_REF]. RWSAT is based on the observation that a violated equation can be satisfied through negation of one of its variables:

• Start from a randomly chosen configuration of the variables. Call energy the number E of unsatisfied equations.

• While E ≥ 1;

pick up uniformly at random one of the E unsatisfied equations; FIG. 11: Fraction e of unsatisfied equations as a function of time t (number of steps divided by N ) during the operation of RWSAT on a random 3-XORSAT formula at ratio α = 0.2 (left) and α = 0.4 (right) with N = 10 3 (dotted), 10 4 (dashed), 10 5 (full curve) variables. Note the difference of horizontal scales between the two figures. Inset: blow up of the t ∈ [9.5; 10.5] region; the amplitude of fluctuations around the plateau decreases with increasing size N .

pick up uniformly at random one of its 3 variables; -negate the value of this variable, update E;

• Print 'Satisfiable', and Halt.

Notice that, as a result of the negation, some equations that were satisfied may become violated. Therefore the energy is not guaranteed to decrease with the number of steps of the algorithm. RWSAT is able to escape from local minima of the energy landscape, and is a priori capable of better performances. From the other hand, RWSAT may run forever... A major question is how long should the algorithm be running before we stop thinking that the studied system has solutions hard to find and get some confidence that there is really no solution.

This question was addressed by Schöning [61], who showed that RWSAT could easily be used as a one-sided randomized algorithm [54][76]. Consider one instance of 3-XORSAT and run RWSAT for 3N steps from a randomly chosen configuration of variables. Choose again a random initial configuration and run RWSAT another 3N steps, and so on ... The probability that no solution has been found after T repetitions of this procedure though the formula is satisfiable is

p SAT ≤ exp -T × 3 4 N +o(N ) . (96)
Hence we obtain a probabilistic proof that the instance is not satisfiable if the algorithm has run unsuccessfully for more than ( 43 ) N sets of 3N steps. It must be clear that this result holds for any instance, no assumption being made on the distribution of formulas. The probability appearing in (96) is on the random choices done by RWSAT and the choices of the restart configurations for a fixed formula.

The proof of (96) can be sketched as follows. Assume that the formula is satisfiable, and called X * one of its solutions. Consider now the (extensive) Hamming distance between the solution and the configuration X of variables produced by RWSAT at some instant. After each step only one variable is changed so D changes into D + 1 (bad move) or D -1 (good move). Call x, y, z the variables in the equation which was not satisfied by X. One or three of those variables have opposite values in X * . In the latter case the flip is always a good move; in the former case the good move happens with probability 1 3 and a bad move with probability 2 3 . On the overall the probability of a good move is 1 3 at least. Think of D has the position of a random walker on the [0; N ] segment. Initially the position of the walker is a binomial variable, centered in N 2 . At each step the walker moves to the left with probability 1 3 , and to the right with probability 2 3 . We look for the probability ρ that the walker is absorbed by the boundary D = 0 after S steps. A standard calculation shows that ρ is maximal for S = 3N , with the value ρ ≃ ( 3 4 ) N . After T repetitions the probability of not having been absorbed is (1ρ) T < exp(-ρ T ), hence (96). The proof can be easily extended to K-XORSAT with higher values of K. The number of repetitions necessary to prove unsatisifiability scales as ( 2(K-1)

K

) N ; it is essentially equal to 2 N for large K, showing that RWSAT does not beat exhaustive search in this limit.

B. Dynamical transition of RWSAT on random XORSAT instances

Result (96) is true for any instance; what is the typical situation for random systems? Numerical experiments indicate that there is critical value of the ratio of equations per variables, α E ≃ 0.33, hereafter referred to as dynamical threshold, separating two regimes:

• for α < α E , RWSAT generally finds a solution very quickly, namely with a number of flips growing linearly with the number of variables N [77]. Figure 11 shows the plot of the fraction e of unsatisfied clauses as a function of the time (number of steps) T for one randomly drawn system with ratio α = 0.2 and N = 500 variables.

The curve shows a fast decrease from the initial value (e(T = 0) = 1 2 independently of α for large values of N , but deviations can be found at small sizes, see Figure 11) down to zero on a time scale of the order of N [78]. The resolution time T res depends both on the system of equations under consideration and the choices of the algorithm; its average value scales as

T res = N t res + o(N ) . (97)
where t res is an increasing function of α[79].

• for systems with ratios of equations per variable in the α E < α < α c range, the initial relaxation regime taking place on the O(N ) time scale does not allow RWSAT to reach a solution (Figure 11B). The fraction e of unsat equations then fluctuates around some plateau value e plateau for a very long time. Fluctuations are smaller and smaller (and the height of the plateau better and better defined) as the size N increases. As a result of fluctuations, the fraction e of unsatisfied equations may temporarily either increase or decrease. When a fluctuation happens to drive RWSAT to e = 0, a solution is found and the algorithm stops. The corresponding resolution time, T res , is stochastic; numerical experiments for different sizes N indicate that its expectation value scale as

T res = exp(N τ res + o(N )) . (98)
where the coefficient τ res is an increasing function of α. The plateau energy e plateau and the logarithm τ res of the resolution time are shown in Figure 12.

Notice that the dynamical threshold α E above which the plateau energy is positive is strictly smaller than the critical threshold α c ≃ 0.918, where systems go from satisfiable with high probability to unsatisfiable with high probability. In the intermediate range α E < α < α c , systems are almost surely satisfiable but RWSAT needs an exponentially large time to prove so. The reason is that RWSAT remains trapped at a high energy level (plateau of Figure 12) for an exponentially large time. The emergence of metastability can be qualitatively studied with simple tools we now expose.

C. Approximate theory for the metastable plateau and the escape time

Assume that after T steps of the algorithm the energy (number of unsatisfied equations) is E T ≥ 1. Then pick up an unsatisfied equation, say, C, and a variable in C, say, x, and flip it. The energy after the flip is

E T +1 = E T -1 -U + S , (99)
where S (respectively U ) is the number of equations including x which were satisfied (resp. unsatisfied after exclusion of equation C) prior to the flip. S and U are random variables with binomial distributions,

Proba[U ] = E T -1 U 3 N U 1 - 3 N ET -1-U , Proba[S] = M -E T S 3 N S 1 - 3 N M-ET -S . (100)
where the probabilities are intended over the formula content. Taking the average evolution equation (99) we obtain

E T +1 = E T -1 - 3 N E T -1 + 3 N M -E T . (101)
The above equation is exact. It is now tempting to iterate it with time, from the initial condition

E T =0 = M 2 .
This is what we do hereafter but one should realize that this procedure is not correct from a mathematical standpoint. The catch is that one is allowed to average over the formula only once, and certainly not at each time step of the algorithm. Evolution equation (101) amounts to redraw randomly the instance at each time step, conditioned to the energy. This approximation nevertheless allows us to write down a simple equation for E T , which captures much of the true behaviour of RWSAT.

The next step in our analysis is the large size, large time limit. As the energy can typically change by a quantity of the order of unity in one time step we expect the fraction of unsatisfied equations to vary of a time scale of the order of N ,

E T = M e T M = t , (102)
for some smooth function e(t) of the reduced time t. Finite difference equation ( 101) turns into a differential equation after insertion of (102),

de dt = -1 + 3α (1 -2 e) , (103)
with the initial condition e(0) = 1 2 . Clearly (103) makes sense as long as e > 0; if e vanishes the algorithm stops. Resolution of (103) shows the following scenario. If α is smaller than

α E = 1 3 , (104)
the fraction e of unsatisfied equations quickly decreases, and vanishes at some time t res (α). This regime corresponds to a successfull action of RWSAT in a O(N ) number of steps. t res is an increasing function of α which diverges as α → α E . Above this critical ratio e shows a different behaviour: after a decreasing transient regime e saturates to a positive plateau value

e plateau (α) = 1 2 1 - α E α . (105)
The value of the plateau energy is compared to numerics in Figure 12A. The agreement on the location of the dynamical threshold α E as well as the plateau energy are satisfactory. The remaining point is to understand how RWSAT finally finds a solution when α > α E . The above theory, based on taking the N → ∞ limit first, washes out the fluctuations of the energy around its metastable value, of crucial importance for resolution [62]. To take into account these fluctuations let us define the probability Q plateau (E) that the energy takes value E in the plateau regime of Figure 11B. A stationary distribution is well defined if we discard the initial transient regime (choose large t) and collect values for E on exponentially large-in-N time scales. The procedure is standard in the study of long-time metastable states.

Within our draw-instance-at-each-step approximation we may write a self-consistent equation for the stationary distribution of energies,

Q plateau (E) = U,S Proba[U ] Proba[S] Q plateau (E + 1 + U -S) (106)
where the meaning of U, S was explained right after (99). From Section V B we expect fluctuations to decreases sharply with the system size. A reasonable guess for the scaling of the distribution with M is

Q plateau (E) = exp -M ω E M = e + o(M ) ( 107
)
where ω is the rate function associated to the fraction of unsatisfied equations. Plugging the above Ansatz into (106) and taking the large M limit we find that ω fulfills the following differential equation

F ∂ω ∂e , e = 0 , (108)
where F (x, y) = 3αy(e -x -1)+3α(1-y)(e x -1)-x. This equation has to be solved with the condition ω(e plateau ) = 0. An analytical solution can be found for (108) when we restrict to the vicinity of the dynamical transition i.e. to small values of ω. Expanding F to the second order in its first argument and solving (108) we obtain

ω(e) ≃ 2 (e -e plateau ) 2 , ( 109
)
where e plateau is defined in (105). What happens when time increases is now clear. Assume we have run RWSAT up to time t ∼ e Mτ . Then configurations with energy e such that ω(e) < τ have been visited many times and are 'equilibrated' with probability (106), (109). Configurations with energies outside the band e plateau ± τ /2 are not accessible. When the time scales reaches

τ res = ω(0) ≃ 1 2 1 - 1 3α 2 , (110)
zero energy configurations are encountered, and RWSAT comes to a stop. The agreement between the theoretical estimate (110) and the numerical findings (98) visible in Figure 12B is acceptable in regard to the crudeness of the approximation done.

D. Davis-Putnam-Loveland-Logemann (DPLL) algorithm

The second procedure is the Davis-Putnam-Loveland-Logemann (DPLL) algorithm [22]. Contrary to RWSAT DPLL can provide exact proofs for unsatisfiability. The procedure, widely used in practice, is based on the trial-and-error principle. Variables are assigned according to some heuristic rule (split step), and equations involving those variables simplified. If an equation involving a single variable (unit-equation) appears its variable is chosen accordingly prior to any other heuristic assignment (unit-propagation). If a contradiction is found (two opposite unit-equations) DPLL backtracks to the last heuristically assigned variable, flips it, and resumes the search process. The procedure halts either when all equations have been satisfied (a solution is then found), or when all possible values for the variables have been tried in vane and found to be contradictory (a proof of unsatisfiability is then obtained).

DPLL can be described as a recursive function of the variable assignment A. Given a system S DPLL is first called with the empty assignment A = ∅:

Procedure DPLL[A]

• Let S A be what is left from S given variable assignment A;

• if S A is empty, Print 'Satisfiable'; Halt; Rules for assigning variables in the absence of unit-equations are heuristic in that they aim at doing good assumptions i.e. diminishing as much as possible the search process to come from limited information about the current system of equations. Of course, perfect heuristic do exist: trying all possible values for not-yet-assigned variables would ensure that no wrong guess is ever done! But the time required would be exponentially long. In practice, heuristics have to make their decision in polynomial time. Two simple splitting heuristics are: ⋄ UC: choose at random and uniformly any unset variable, and assign it to 0 or 1 with equal probabilities ( 12 ). ⋄ GUC: choose at random and uniformly any equation with minimal length i.e. involving 2 variables if any, or 3 variables otherwise. Pick up at random and uniformly one its variable, and assign it to 0 or 1 with equal probabilities ( 1 2 ). UC, which stands for unit-clause [START_REF] Chao | [END_REF], amounts to make a random guess and is the simplest possible heuristic. GUC (Generalized UC) is more clever: each time a split is done from an equation with 2 variables, this equation is turned into a unit-equation, and eliminated through unit-propagation. In the following, we call DPLL-UC and DPLL-GUC the variants of DPLL based on the UC and GUC heuristics respectively.

A measure of the computational effort required by DPLL is the number T split of variable splittings. This number varies from system to system (at fixed number N of variables and ratio α), and from run to run of DPLL due to the stochasticity introduced by the heuristic rule. The outcome of numerical experiments for the median number of splits [80]. For a given size N T split shows a maximum located around α ≃ α c . If one fixes α T split is an increasing function of the size N ; numerical data support the existence of a dynamical threshold, α E , separating linear and exponential scalings in N ,

T split ∼ N t split + o(N ) if α < α E exp(N τ split + o(N )) if α > α E , (111)
where t split and τ split are functions of the ratio α. The value of the dynamical threshold can be derived from theoretical calculations shown in Section V E and is equal to α E = 2 3 and α E ≃ 0.7507... for UC and GUC heuristics respectively. Three dynamical regimes are therefored identified [2,16]:

• Linear & satisfiable phase (α < α E ): systems with small ratios are solved with essentially no backtracking. A solution is found after O(N ) splits.

• Exponential& satisfiable phase (α E < α < α c ): systems with ratios slightly below threshold have solutions, but DPLL generally requires an exponential number of splits to find one of them. An explanation for this drastic breakdown of performances will be given in Section V E.

• Exponential & unsatisfiable phase (α > α c ): finally, finding a proof of unsatisfiability typically requires an exponentially large number of splits [15]. Note that, as α gets higher and higher, each variable assignment affects more and more equations (of the order of α), and contradictions are detected earlier and earlier. Rigorous calculations show that τ split ∼ 1 α [START_REF] Beame | Proceedings of the ACM Symp. on Theory of Computing[END_REF], and the computational effort decreases with increasing α (Figure 13). The median number of splits is considerably smaller for DPLL-GUC than for DPLL-UC, a result expected from the advantages of GUC against UC discussed above. 3 for the UC heuristic) trajectories remain confined in the sat phase, end in S of coordinates (0, 0), where a solution is found. At αE the trajectory hits tangentially the threshold line in T of coordinates ( 14 , 2 3 ). When α > αE the trajectories intersect the threshold line at some point G (which depends on α), and stops before hitting the αD(p) dotted line (122). After massive backtracking DPLL will find a solution; G corresponds to the highest node in the search tree.

E. Linear phase: resolution trajectories in the 2 + p-XORSAT phase diagram Action of DPLL on an instance of 3-XORSAT causes changes to the numbers of variables and equationses, and thus to the ratio α. Furthermore DPLL turns equations with 3 variables into equation with 2 variables. A mixed 2 + p-XORSAT distribution, where p is the fraction of 3-equations and α the ratio of the total number of 2-and 3equations over the number of variables can be used to model what remains of the input system [81]. Repeating the calculations of Section III for the 2 + p-XORSAT models we derive the phase diagram of Figure 14. The Sat/Unsat critical line α c (p) separates the satisfiable from the unsatisfiable phases. For p ≤ p 0 = 1 4 i.e. to the left of point T, the threshold line coincides with the percolation transition as in the 2-XORSAT model, and is given by α c (p) = 1 2(1-p) . For p > p 0 an intermediate clustered phase is found as in the 3-XORSAT model, and the threshold coincides with the vanishing of the cluster entropy s cluster (Section II J).

The phase diagram of 2+p-XORSAT is the natural space in which DPLL dynamic takes place. An input 3-XORSAT instance with ratio α shows up on the right vertical boundary of Figure 14 as a point of coordinates (p = 1, α). Under the action of DPLL the representative point moves aside from the 3-XORSAT axis and follows a trajectory, very much alike real-space renormalization, which depends on the splitting heuristic. Trajectories enjoy two essential features [1]. First the representative point of the system treated by DPLL does not 'leave' the 2+p-XORSAT phase diagram. In other words, the instance is, at any stage of the search process, uniformly distributed from the 2+p-XORSAT distribution conditioned to its equation per variable ratio α and fraction p of 3-equations. This assumption is not true for all heuristics of split, but holds for UC and GUC [START_REF] Chao | [END_REF][82]. Secondly, the trajectory followed by an instance in the course of resolution is a stochastic object, due to the randomness of the instance and of the assignments done by DPLL. In the large size limit (N → ∞) the trajectory becomes self-averageing i.e. concentrated around its average locus in the 2+p-XORSAT phase diagram [63]. We will come back below on this concentration phenomenon.

Let α 0 denote the equation per variable ratio of the 3-XORSAT instance to be solved. We call E j (T ) the number of j-equations (including j variables) after T variables have been assigned by the solving procedure. T will be called hereafter 'time', not to be confused with the computational effort. At time T = 0 we have E 3 (0) = α 0 N , E 2 (0) = E 1 (0) = 0. Assume that the variable x assigned at time T is chosen through unit-propagation, that is, independently of the j-equation content. Call n j (T ) the number of occurrences of x in j-equations (j = 2, 3). The evolution equations for the populations of 2-,3-equations read

E 3 (T + 1) = E 3 (T ) -n 3 (T ) , E 2 (T + 1) = E 2 (T ) -n 2 (T ) + n 3 (T ) . (112)
Flows n 2 , n 3 are of course random variables that depend on the instance under consideration at time T , and on the choice of variable done by DPLL. What are their distributions? At time T there remain N -T untouched variables;

x appears in any of the E j (T ) j-equation with probability p j = j N -T , independently of the other equations. In the large N limit and at fixed fraction of assigned variables, t = T N , the binomial distribution converges to a Poisson law with mean

n j T = j e j 1 -t where e j = E j (T ) N (113)
is the density of j-equations at time T . The key remark is that, when N → ∞, e j is a slowly varying and non stochastic quantity and is a function of the fraction t = T N rather than T itself. Let us iterate (112) between times T 0 = t N and T 0 + ∆T where 1 ≪ ∆T ≪ N e.g. ∆T = O( √ N ). Then the change ∆E 3 in the number of 3-equations is (minus) the sum of the stochastic variables n j (T ) for T = T 0 , T 0 + 1, . . . , T 0 + ∆T . As these variables are uncorrelated Poisson variables with O(1) mean (113) ∆E 3 will be of the order of ∆T , and the change in the density e 3 will be of order of ∆T /N → 0. Applying central limit theorem ∆E 3 /∆T will be almost surely equal ton 3 t given by ( 113) and with the equation density measured at reduced time t. The argument can be extended to 2-equations, and we conclude that e 2 , e 3 are deterministic (self-averaging) quantities obeying the two coupled differential equations [START_REF] Chao | [END_REF]

de 3 dt (t) = - 3 e 3 1 -t , de 3 dt (t) = 3 e 3 1 -t - 2 e 2 1 -t . (114)
Those equations, together with the initial condition e 3 (0) = α 0 , e 2 (0) = 0 can be easily solved,

e 3 (t) = α 0 (1 -t) 3 , e 2 (t) = 3 α 0 t (1 -t) 2 . ( 115
)
To sum up, the dynamical evolution of the equation populations may be seen as a slow and deterministic evolution of the equation densities to which are superimposed fast, small fluctuations. The distribution of the fluctuations adiabatically follows the slow trajectory. This scenario is pictured in Figure 15. Expressions (115) for the equation densities allow us to draw the resolution trajectories corresponding to the action of DPLL on a 3-XORSAT instance. Initially the instance is represented by a point with coordinates (p = 1, α = α 0 ) in Figure 14. As more and more variables are assigned the representative point moves away from the rightmost vertical axis. After a fraction t of variables have been assigned the coordinates of the point are

p(t) = e 3 e 2 + e 3 = 1 -t 1 + 2t , α(t) = e 2 + e 3 1 -t = α 0 (1 -t)(1 + 2t) . (116)
Trajectories corresponding to various initial ratios are shown in Figure 14. For small ratios α 0 < α E trajectories remain confined in the sat phase, end in S of coordinates (0, 0), where a solution is found. At α E (= 2 3 for the UC heuristic), the single branch trajectory hits tangentially the threshold line in T of coordinates ( 14 , 2 3 ). When α 0 > α E the trajectories enter the Unsat phase, meaning that DPLL has turned a satisfiable instance (if α 0 < α c ) into an unsatisfiable one as a result of poor assignments. It is natural to expect that α E is the highest ratio at which DPLL succeeds in finding a solution without resorting to much backtracking.

F. Dynamics of unit-equations and universality

The trajectories we have derived in the previous Section are correct provided no contradiction emerges. But contradictions may happen as soon as there are E 1 = 2 unit-equations, and are all the more likely than E 1 is large. 16: Evolution of the number E1 of 1-equations as one more variable is assigned. n2 denotes the number of 2-equations reduced to 1-equations, s1 the number of 1-equations satisfied. If E1 ≥ 1 a variable is fixed through unit-propagation: E1 decreases by one plus s1, and increases by n2. In the absence of unit-equation (E1 = 0) the number of 1-equations after the assignment is simply

E' 1 E' 0 n 2 E 1 1 s 1 2 n FIG.
E ′ 1 = n2.
Actually the set of 1-equations form a 1-XORSAT instance which is unsatisfiable with a finite probability as soon as E 1 is of the order of √ N from the results of Section II A. Assume now that E 1 (T ) ≪ N after T variables have been assigned, what is the probability ρ T that no contradiction emerges when the T th variable is assigned by DPLL? This probability is clearly one when E 1 = 0. When E 1 ≥ 1 we pick up a 1-equation, say, x 6 = 1, and wonder whether the opposite 1-equation, x 6 = 0, is present among the (E 1 -1) 1-equations left. As equations are uniformly distributed over the set of N -T untouched variables

ρ T = 1 - 1 2(N -T ) max(E1(T )-1,0) . ( 117
)
The presence of the max in the above equation ensures it remains correct even in the absence of unit-equations (E 1 = 0). E 1 (T ) is a stochastic variable. However from the decoupling between fast and slow time scales sketched in Figure 15 the probability distribution of E 1 (T ) depends only on the slow time scale t. Let us call µ(E 1 ; t) this probability. Multiplying (117) over the times T = 0 to T = N -1 we deduce the probability that DPLL has successfully found a solution without ever backtracking,

ρ success = exp   - 1 0 dt 2(1 -t) E1≥1 µ(E 1 ; t) (E 1 -1)   (118)
in the large N limit. We are left with the calculation of µ [START_REF] Frieze | [END_REF]. Figure 16 sketches the stochastic evolution of the number E 1 during one step. The number of 1-equations produced from 2-equations, n 2 , is a Poisson variable with average value, from (115),

d(t) = 2 e 2 (t) 1 -t = 6 α 0 t(1 -t) (119)
when N → ∞. The number of satisfied 1-equations, s 1 , is negligible as long as E 1 remains bounded. The probability that the number of 1-equations goes from E 1 to E ′ 1 when T → T + 1 defines the entry of the transition matrix

M (E ′ 1 , E 1 ; t) = n2≥0 e -d(t) d(t) n2 n 2 ! δ E ′ 1 -(E1+n2-δE 1 ) . (120)
from which a master equation for the probability of E 1 at time T may be written. On time scales 1 ≪ ∆T ≪ N this master equation converges to the equilibrium distribution µ [16,[START_REF] Frieze | [END_REF], conveniently expressed in terms of the generating function

G(x; t) = E1≥0 µ(E 1 ; t) x E1 = (1 -d(t))(x -1) x e d(t) (1-x) -1 . ( 121
)
The above is a sensible result for d(t) ≤ 1 but does not make sense when d(t) > 1 since a probability cannot be negative! The reason is that we have derived (121) under the implicit condition that no contradiction was encountered. This assumption cannot hold when the average rate of 1-equation production, d(t), is larger that one, the rate at which 1-equations are satisfed by unit-propagation. From (119) we see, when α > α E = 2 3 , the trajectory would cross the generated from a branch with population E once the T th variable is assigned. Transition matrix M is an extension of (120) to the whole population vector E and not only E 1 . We have 0 ≤ M ≤ 2, the extreme values corresponding to a contradiction and to a split respectively. We claim that

α D (p) = 1 2(1 -p) (122)
L( E ′ ; T + 1) = E M ( E ′ , E; T ) L( E; T ) . (126)
Evolution equation ( 126) is somewhat suspicious since it looks like the approximation (106) we have done in the analysis of RWSAT. Yet a major difference exists which makes (126) exact [51]. Drawing randomly many times the same instance, as we are doing, is in principle forbidden but not along one branch for the very reason the analysis of Section V E was correct. Actually what we have done in Section V E is to draw randomly at time T the equations containing the T th variable. But this is correct since those equations are immediately simplified into shorter equations and their remaining content remains unknown [START_REF] Kaporis | How to prove conditional randomness using the principle of deferred decisions[END_REF]. The situation seems more complicated in the case of the whole tree since the same equation can appear at different depth along distinct branches. Indeed the number of branches produced from two distinct branches after assignment of one variable are correlated variables. But thanks to the linearity of expectation those correlations do not matter and (126) is correct.

Transition matrix M can be explicitely written down. It is more convenient to write (126) for the generating function of the number of branches, B( x; T ) = E L( E; T )x E1

where f is the vector with components

f 1 = x 1 + 1 2 -x 1 N -T , f 2 = x 2 + 2(x 2 -x 1 ) N -T , f 3 = x 3 + 3(x 2 -x 3 ) N -T . (128)
The three terms on the r.h.s. of (127) correspond, from left to right: unit-propagation (the branch keeps growing), variable splitting (2 branches are created from the previous one), branches carrying empty instances (satisfied instance). Equation (127) together with the initial condition B( x; 0) = x α0N 3 completely defines the average dynamics of the search tree. We sketch the main steps of its resolution below [16]:

1. To count the number of branches irrespectively of the number of unit-equations we should consider the value x 1 = 1. However, as long as branches grow the number E 1 of unit-equations cannot be large, and remains bounded. We can therefore choose x 1 = 1 2 which simplifies (127) without affecting the large size scaling of L and B. This technical trick is reminiscent of Knuth's kernel method [START_REF] Knuth | The Art of Computer Programming[END_REF].

2. For large N it is reasonable to expect that the number of branches grows exponentially with the depth, or, equivalently,

E1

L(E with the initial condition γ(x 2 , x 3 ; t) = α 0 ln x 3 .

3. The first order PDE can be solved exactly with the characteristic method. The output, after Legendre inversion through (130), is the entropy λ(e 2 , e 3 ; t) of branches at reduced depth t. Let us call λ * (t) the maximum value of λ over the equation densities for a fixed fraction t of assigned variables.

4. λ * (t) is a function growing from λ * = 0 at t = 0, reaching a maximum value λ * M in t M , and decreasing for larger times t ≤ 1. t M is the depth in the tree of Figure 17C where most contradictions are found; the number of C leaves is, to exponential order, e N λ * M . We conclude that the size of the tree we were looking for is

τ split = λ * M , (132)
compare with (111). For large α ≫ α c one finds τ split ∼ ln 2/(6α) in agreement with [START_REF] Beame | Proceedings of the ACM Symp. on Theory of Computing[END_REF]. The calculation can be extended to highers values of K.

The above calculation holds for the unsatisfiable, exponential phase. How can we understand the satisfiable but exponential regime α E < α 0 < α c ? The resolution trajectory crosses the Sat/Unsat critical line at some point G shown in Figure 14. Immediately after G the instance left by DPLL is unsatisfiable. A subtree with all its leaves carrying contradictions will develop below G (Figure 17B). The size τ G split of this subtree can be easily calculated from the above theory. The only change is the initial condition over γ: γ(x 2 , x 3 ; 0) = α G (p G ln x 3 + (1p G ) ln x 2 ) where (p G , α G ) are the coordinates of G which can be calculated from α 0 and the knowledge of the critical Sat/Unsat line. Once this subtree has been built DPLL backtracks to G, flips the attached variable and will finally end up with a solution. Hence the (log of the) number of splits necessary will be typically equal to τ split = (1t G ) τ G split [16].

VI. CONCLUSIONS

Previous Sections have allowed us to illustrate rather general techniques and ideas to deal with random systems. It does not come as a surprise that other problems than XORSAT e.g. the satisfaction of Boolean constraints, graph coloring, the covering of vertices, ... have been successfully studied with these tools. Many of those problems, when given an input distribution based on random graphs, actually share a lot of common features with XORSAT. The reader is referred to [3,[START_REF] Biroli | [END_REF][START_REF] Krzakala | Gibbs States and the Set of Solutions of Random Constraint Satisfaction Problems[END_REF]40,45,49] (satisfiability), [START_REF] Krzakala | Gibbs States and the Set of Solutions of Random Constraint Satisfaction Problems[END_REF][START_REF] Mulet | [END_REF] (coloring), [31,64] (vertex cover), ... for entry points to the literature. Let us also mention that many other interesting optimization problems, not directly related to random graphs, have been studied with the techniques of Sections 4 and 5, and the results sometimes rigorously proven e.g. matching [5,41,56], traveling salesman [42], number partitioning [START_REF] Mertens | [END_REF]39], graph partitioning [30], ... Finally, from a historical point of view, one should not forget that statistical mechanics tools have found numerous and beautiful applications in the study of the learning and storage properties of neural networks [START_REF] Amit | Modeling Brain Function[END_REF][START_REF] Engel | Statistical Mechanics of Learning[END_REF], all the more so the random satisfiability problem can be recast as an Ising perceptron problem [START_REF] Krauth | [END_REF].

The study of random optimization problems is obviously interesting from a probabilistic point of view. As far as computer science is concerned they can be seen as useful benchmarks for testing and improving resolution procedures. A successful example is the traduction of the cavity equations of Section 5 into an algorithm for solving given instances of the satisfiability problem [45]. This algorithm, called Survey Propagation, extends to the clustered phase the Belief Propagation procedure of wide-spread use in statistical inference, and is a very efficient procedure to find solutions to 3-Satisfiability slightly below threshold. Another application of statistical physics ideas is the conception of new heuristics for DPLL capable of proving the unsatisfiability of formulas with 700 hundreds variables at threshold [START_REF] Dequen | Proceedings of Theory and Applications of Satisfiability Testing, 6th International Conference[END_REF].

Despite those successes important question remain open. First is there a relationship between clustering and hardness of resolution? This question is reminiscent of a very general issue in statistical physics, namely the relationship between dynamical and static properties of disordered or glassy systems [21]. The onset of clustering, or more precisely of strong correlations between variables over the space of solutions drastically worsens the performances of sampling algorithms e.g. Monte Carlo procedures [START_REF] Krzakala | Gibbs States and the Set of Solutions of Random Constraint Satisfaction Problems[END_REF][START_REF] Montanari | [END_REF]. However, in practical applications, one looks for a solution rather than for the sampling of the solution space... From this point of view knowing whether solutions are clustered or not does not seem to be of crucial relevance. Actually a local and polynomial search strategy capable of finding solutions well above the clustering threshold has been explicitely found for various optimizations problems [START_REF] Krzakala | A landscape analysis of constraint satisfaction problems[END_REF].

Another open question is what happens at large K, that is, when constraints involve more and more variables. The performances of all known algorithms, be they local search procedures or DPLL solvers, seem to deteriorate. Worstcase bound indicate that the large K case is very difficult [START_REF] Impagliazzo | Proceedings of the IEEE Conference on Computational Complexity[END_REF]. From statistical mechanics point of view problems look like more and more the random energy model [START_REF] Mézard | Spin glasses and beyond[END_REF] as K increases, but can we beat the worst-case bounds on average? Finally let us mention a recent work by Feige [START_REF] Feige | Proceedings of 34th STOC conference[END_REF] which, for the first time, showed that the complexity of solving random SAT (or XORSAT) model had a fundamental interest in worst-case approximation theory. Consider 3-SAT instances with ratio α ≫ α c . Most of them have GS energy close to αN/2, but a very tiny fraction of those instances have energy smaller than, say, ǫN where ǫ ≪ α is fixed. Is there a polynomial algorithm capable of recognizing all such atypical formulas from the vast majority of typical instances? Insights from statistical physics suggest that, the answer is positive for SAT (if we want most satisfiable instances to be detected and not all of them) while XORSAT seems to be much harder [START_REF] Altarelli | [END_REF]! Actually, to the knowledge of the author, no local search algorithm (based on random walk, variable assigment, Monte Carlo, message-passing, cooling procedure, ...) is efficient for solving XORSAT. This makes the study of this problem even more valuable from a computer science point of view.

FIG. 1 :

 1 FIG.1: Left: Probability that a random 1-XORSAT formula is satisfiable as a function of the ratio α of equations per variable, and for various sizes N . Right: same data as in the left panel after the horizontal rescaling α → α × √ N = y; note the use of a log scale for the vertical axis. The dashed line shows the scaling function Φ1(y)[START_REF] Amit | Modeling Brain Function[END_REF].

FIG. 2 :

 2 FIG.2: Same data as Figure1 (left) with: logarithmic scale on the vertical axis, and rescaling by -1/N . The scaling functions ω1[START_REF] Chao | [END_REF] and ω p 1[START_REF] Bollobas | Random Graphs[END_REF] for, respectively, the fixed-size and fixed-probability ensembles are shown.

F.

  Large deviations for PSAT (II): bounds in the Unsat phase of 2-XORSAT.

FIG. 5 :

 5 FIG.5: Rate function ω2(α) associated to the probability of satisfaction of 2-XORSAT formulas with ratio α. The dotted line is the upper bound[START_REF] Feige | Proceedings of 34th STOC conference[END_REF] and the dot-dashed line the lower bound[START_REF] Dequen | Proceedings of Theory and Applications of Satisfiability Testing, 6th International Conference[END_REF]. The full line is the output of the replica calculation of Section III D, squares and circles represent numerical results for N = 200, 100 respectively from 10 6 formulas. Inset: PSAT as a function of the size N at the Sat/Unsat ratio. The slope -1 12 (39) is shown for comparison.

FIG. 7 :

 7 FIG.7: Entropies (base 2 logarithms divided by size N ) of the numbers of solutions and clusters as a function of the ratio α. The entropy of solutions equals 1α for α < αc ≃ 0.918. For α < α d ≃ 0.818, solutions are uniformly scattered on the N -dimensional hypercube. At α d the solution space discontinuously breaks into disjoint clusters. The entropies of clusters, s cluster , and of solutions in each cluster, sin, are such that s cluster + sin = s. At αc the number of clusters stops being exponentially large (s cluster = 0). Above αc there is almost surely no solution.

FIG. 8 :

 8 FIG.8: Phase diagram of 3-XORSAT. A 'geometrical' phase transition takes place in the satisfiable phase at α d ≃ 0.818. At small ratios α < α d solutions are uniformely scattered on the N -dimensional hypercube, with a typical normalized Hamming distance d = 1 2 . At α d the solution space discontinuously breaks into disjoint clusters: the Hamming distance din ≃ 0.14 between solutions inside a cluster is much smaller than the typical distance d clu = 1 2 between two clusters.

FIG. 10 :

 10 FIG.10: Rate function ω3 for the probability of satisfaction of full (bottom curve) and entropy s * h of solutions for homogeneous (top curve) 3-XORSAT systems vs. α. The vertical dotted lines indicate the critical Sat/Unsat threshold, αc ≃ 0.918. For α < αc ω3 = 0, and s * h = (1α) ln 2 is the same as for full systems. Above the threshold ω * < 0. Homogeneous systems are, of course, always satisfiable: the entropy s * h is a positive but quickly decreasing function of α.

  modification has to be brought to the cavity assumption of Section IV A is now clear. Given a subset I of the spins with configuration S I we define E (c) F (S I ) as the GS energy over configurations in the cluster (c). Then the cavity assumption is correct (spins in I are uncorrelated) and E (c) F define the fields h (c)

FIG. 12 :

 12 FIG.12: Fraction e plateau of unsatisfied equations on the plateau (A) and logarithm τres of the average resolution time divided by N (B) as a function of the ratio α of equations per variable. Diamonds are the output of numerical experiments, and have been obtained through average of data from simulations over 1,000 systems and runs of RWSAT for various sizes N , and extrapolation to N → ∞[62]. Full lines are theoretical approximations (105),(110).

•FIG. 13 :

 13 FIG.13: Median number of splits required by DPLL with the UC (left) and GUC (right) heuristics as a function of the ratio α, and for N = 20, 40, 60 variables (from bottom to top). Data have been extracted from the resolution of 10,000 randomly drawn systems; continuous lines are guidelines for the eye. Note the difference of (logarithmic) scale between UC and GUC curves showing that DPLL-GUC is much more efficient than DPLL-UC. The polynomial/exponential transition is located at ratios αE =2 3 and αE = 0.7507... for UC and GUC respectively.

FIG. 14 :

 14 FIG.14: Phase diagram of 2 + p-XORSAT and dynamical trajectories of DPLL. The threshold line αc(p) (bold full line) separates sat from unsat phases. Departure points for DPLL trajectories are located on the 3-XORSAT vertical axis with ratios .4,2 3 , .8, 1. from bottom to top. The arrow indicates the direction of motion along trajectories parametrized by the fraction t of variables set by DPLL. For small ratios α < αE (=2 3 for the UC heuristic) trajectories remain confined in the sat phase, end in S of coordinates (0, 0), where a solution is found. At αE the trajectory hits tangentially the threshold line in T of coordinates (14 , 2 3 ). When α > αE the trajectories intersect the threshold line at some point G (which depends on α), and stops before hitting the αD(p) dotted line (122). After massive backtracking DPLL will find a solution; G corresponds to the highest node in the search tree.

FIG. 15 :

 15 FIG.15: Deterministic versus stochastic dynamics of the equation population E as a function of the number of steps T of the algorithm. On the slow time scale (fraction t = T /N ) the density e = E/N of (2-or 3-) equations varies smoothly according to a deterministic law. Blowing up of the dynamics around some point t ′ , e ′ shows the existence of small and fast fluctuations around this trajectory. Fluctuations are stochastic but their probability distribution depends upon the slow variables t ′ , e ′ only.

1

 1

FIG. 17 :

 17 FIG. 17: Search trees in three regimes of Section V D: A. linear, satisfiable (α < αE); B. exponential, satisfiable (αE < α < αc); C. exponential, unsatisfiable (α > αc). Leaves are marked with S (solutions) or C (contradictions). G is the highest node to which DPLL backtracks, see Figure14.

  f 2 , f 3 ; T ) -2B( 0; T ) ,

  Introduction to Phase Transitions in Random Optimization Problems R. Monasson Laboratoire de Physique Théorique de l'ENS 24 rue Lhomond, 75005 Paris, France Notes of the lectures delivered in Les Houches during the Summer School on Complex Systems (July 2006).

	I. INTRODUCTION
	A. Preamble

  2 , e 3 are the densities of equations as usual. From point 1 the Legendre transform of λ γ(x 2 , x 3 ; t) = max e2,e3 λ(e 2 , e 2 ; t) + e 2 ln x 2 + e 3 ln x 3

									(130)
	fulfills the partial differential equation (PDE)
	∂γ ∂t	= ln 2 +	1 -2x 2 1 -t	∂γ ∂x 2	+	3(x 2 -x 3 ) 1 -t	∂γ ∂x 3	.	(131)

1 , E 2 , E 3 ; T ) ∼ e N λ(e2,e2;t)+o(N ) (129)

where e

on which d = 1 for some time t D < 1. A contradiction is very likely to emerge before the crossing.

When α < α E d remains smaller than unity at any time. In this regime the probability of success reads, using (118) and (121),

ρ success is a decreasing function of the ratio α, down from unity for α = 0 to zero for α = α E . The present analysis of the UC heuristic can be easily transposed to the GUC heuristic. Details are not given here but can be found in [2,16]. The result is an expression for ρ success larger than its UC counterpart (123), and vanishing in α E ≃ 0.7507. Interestingly the way ρ success vanishes when α reaches α E ,

is the same for both heuristics. This similarity extends to a whole class of heuristics which can be described by the flow of equation densities only and based on unit-propagation [START_REF] Deroulers | [END_REF]. The probability that DPLL finds a solution without backtracking to a 3-XORSAT instance of size N satisfies finite-size scaling at the dynamical critical point,

where the scaling function Φ is independent of the heuristics and can be calculated exactly [START_REF] Deroulers | [END_REF]. The exponent characterizing the width of the critical region is the one associated to percolation in random graphs [START_REF] Mertens | [END_REF]. A consequence of (125) is that, right at α E , ρ success ∼ exp(-Cst × N 6 ) decreases as a stretched exponential of the size. The value of the exponent, and its robustness against the splitting heuristics can be understood from the following argument [START_REF] Deroulers | [END_REF].

Let us represent 1-and 2-equations by a graph G over the set of N -T vertices (one for each variable x i ) with E 1 marked vertices (one for each unit-equation x i = 0, 1), and E 2 signed edges (x i + x j = 0, 1), see Section II D. d is simply the average degree of vertices in G. Unit-propagation corresponds to removing a marked vertex (and its attached edges), after having marked its neighbours; the process is iterated until the connected component is entirely removed (no vertex is marked). Meanwhile, new edges have been created from the reduction of 3-equations into 2-equations. Then a vertex is picked up according to the heuristic and marked, and unit-propagation resumes. The success/failure transition coincides with the percolation transition on G: d = 1 as expected. From random graph theory [START_REF] Bollobas | Random Graphs[END_REF] the percolation critical window is of width |d -1| ∼ N -1/3 . As d is proportional to the ratio α 0 (119) we find back ψ = 1 3 . The time spent by resolution trajectories in the critical window is ∆t ∼ |d -1| ∼ N -1/6 , corresponding to ∆T = N ∆t ∼ N 5/6 eliminated variables. As the largest components have size S ∼ N 2/3 the number of such components eliminated is C = ∆T /S ∼ N 1/6 . What is the probability q that a large component is removed without encountering a contradiction? During the removal of the component the number of marked vertices 'freely' diffuses, and reaches E 1 ∼ √ S ∼ N 1/3 . The probability that no contradiction occurs is, from (117), q ∼ (1-Cst N ) E1×S , a finite quantity. Thus ρ success ∼ q C ∼ exp(-N 1/6 ). The presence of numerous, smaller components does not affect this scaling.

G. Exponential phase: massive backtracking

For ratios α 0 > α E DPLL is very likely to find a contradiction. Backtracking enters into play, and is responsible for the drastic slowing down of the algorithm (Figure 13).

The history of the search process can be represented by a search tree, where the nodes represent the variables assigned, and the descending edges their values (Figure 17). The leaves of the tree correspond to solutions (S), or to contradictions (C). The analysis of the α < α E regime leads us to the conclusion that search trees look like Figure 17A at small ratios [83]. Consider now the case of unsatisfiable formulas (α 0 > α c ) where all leaves carry contradictions after DPLL halts (Figure 17C). DPLL builds the tree in a sequential manner, adding nodes and edges one after the other, and completing branches through backtracking steps. We can think of the same search tree built in a parallel way [16]. At time (depth T ) our tree is composed of L(T ) ≤ 2 T branches, each carrying a partial assignment over T variables. Step T consists in assigning one more variable to each branch, according to DPLL rules, that is, through unit-propagation or split. Possible consequences are: emergence of a contradiction and end of the branch, simplification of the attached formulas and the branch keeps growing.

The number of branches L(T ) is a stochastic variable. Its average value can be calculated as follows [51]. Let us define the average number L( E; T ) of branches with equation populations E = (E 1 , E 2 , E 3 ) at depth T . Initially L( E; 0) = 1 for E = (0, 0, α 0 N ), 0 otherwise. Call M ( E ′ , E; T ) the average number of branches with population E ′ APPENDIX A: A PRIMER ON LARGE DEVIATIONS Large deviation theory is the field of probability which deals with very unlikely events [START_REF] Dembo | Large deviations techniques and applications[END_REF]. You are given a fair (unbiased) coin and toss it N times. The number H of head draws has probability

When N gets large H is highly concentrated around H * = N/2 with small relative fluctuations of the order of O( √ N ). Yet we can ask for the probability of observing a fraction h = H/N equal to say, 25%, of heads, far away from the likely value h * = 50%. To calculate this probability we use Stirling's asymptotic expression for the binomial coefficient in (A1) to obtain

where

is called rate function. The meaning of (A2) is that events with value of h = h * are exponentially rare in N , and ω(h) give the decay (rate) exponent. The answer to our question is e -N ω(.25) ∼ e -0.13 N when N is large. Some comments are:

• ω(h) is strictly positive, except in h = h * = 1 2 where it vanishes. This is the only value for the fraction of head draws with non exponentially small-in-N probability.

• Let h = h * + δh where δh is small. Using ω(h * ) = ω ′ (h * ) = 0 we have

that is, δh is Gaussianly distributed with zero mean and variance (N ω ′′ (h * )) -1 = (4N ) -1 . Hence central limit theorem is found back from the parabolic behaviour of the rate function around its minimum[84].

• ω is here a convex function of its argument. This property is true rate functions describing independent events. Indeed, suppose we have H positive (according to some criterion e.g. being a head for a coin) events among a set of N events, then another set of N ′ events among which H ′ are positive. If the two sets are uncorrelated

since the same total number H + H ′ of positive events could be observed in another combination of N + N ′ events. Taking the logarithm and defining

for any u ∈ [0; 1]. Hence the representative curve of ω lies below the chord joining any two points on this curve, and ω is convex. Non-convex rate functions are found in presence of strong correlations [85].

APPENDIX B: INEQUALITIES OF FIRST AND SECOND MOMENTS

Let N be a random variable taking values on the positive integers, and call p N its probability. We denote by N and N 2 the first and second moments of N (assumed to be finite), and write

the probability that N is not equal to zero. Our aim is to show the inequalities

The right inequality, call 'first moment inequality', is straightforward:

Consider now the linear space made of vectors v = (v 0 , v 1 , v 2 , . . .} whose components are labelled by positive integers, with the scalar product

The left inequality in (B2) is simply the Cauchy-Schwarz inequality for v, v ′ : (v

APPENDIX C: CORRECTIONS TO THE SADDLE-POINT CALCULATION OF N 2

In this Appendix we show that N 2 is asymptotically equivalent to N 2 , where N is the number of solutions of a 3-XORSAT formula with ratio α < α 2 ≃ 0.889. This requires to take care of the finite-size corrections around the saddle-point calculations of Section II F. Let Z = (z 1 , z 2 , . . . , z N ) denotes a configuration of variables at distance d from the zero configuration i.e. dN variables z i are equal to 1, the other (1d)N variables are null. Let q(d, N ) be the probability that Z satisifies the equation z i + z j + z k = 0 where (i, j, k) is a random triplet of distinct integers (unbiased distribution):

and q(d) is defined in [START_REF] Dubois | Proc. of the 43rd annual IEEE symposium on Foundations of Computer Science[END_REF] with K = 3. Terms of the order of N -2 have been discarded.

Using formula [START_REF] Engel | Statistical Mechanics of Learning[END_REF] with q(d) substituted with q(d, N ) and the Stirling formula for the asymptotic behaviour of combinatorial coefficients we have

where A(d, α) is defined in [START_REF] Feige | Proceedings of 34th STOC conference[END_REF], and ∼ indicates a true asymptotic equivalence (no multiplicative factor omitted). The r.h.s. of (C2) is the Riemann sum associated to the integral

We now estimate the integral through the saddle-point method. For α < α 2 ≃ 0.889 the dominant contribution to the integral comes from the vicinity of d * = 1 2 . There are quadratic fluctuations around this saddle-point, with a variance equal to N times the inverse of (the modulus of) the second derivative A dd of A with respect to d. Carrying out the Gaussian integral over those fluctuations we obtain

Therefore, from the second moment inequality, P SAT → 1 when N → ∞ at ratios smaller than α 2 .

equations.

[67] The Birthday problem is a classical elementary probability problem: given a class with M students, what is the probability that at least two of them have the same birthday date? The answer for M = 25 is p ≃ 57%, while a much lower value is expected on intuitive grounds when M is much smaller than the number N = 365 of days in a year.

[68] M obeys a Poisson law with parameter M . Using Stirling formula,

where Ω is defined in [START_REF] Broder | Proceedings of Symposium of Discrete Algorithms (SODA)[END_REF].

[69] The probability that such a graph exists is bounded from above by the average number, see B.

[70] Start from one leaf, assign the attached variable to 0, propagate to the next variable according to the edge value, and so on, up to the completion of the tree.

[71] The equation is satisfied if the number of its variables taking opposite values in Y as in X is even. definition of d the probability (over its index i) that a variable takes different value in X and Y is d. Hence expression [START_REF] Dubois | Proc. of the 43rd annual IEEE symposium on Foundations of Computer Science[END_REF] for q(d). Beware of the O( 1 N ) corrections to this expression e.g. if variable x1 = y1 (which happens with probability d) then the probability that x2 = y2 is (dN -1)/(N -1) = d + (1d)/ -N -1). Those corrections are relevant for the calculation of Gaussian fluctuations around the saddle-point (C).

[72] This scaling is correct provided there is no diverging e.g. O(ln ln N ) corrections to [START_REF] Krzakala | A landscape analysis of constraint satisfaction problems[END_REF].

[73] Actually the identity holds for q = 0 too, and is known under the name of harmonic mean formula [4].

[74] Actually b * is chosen to minimize A h when q < 0, thus λL has always the right negative sign.

[75] An alternative formulation is, for finite size N , that the shortest loops (in extensive number) have lengths of the order of log N [START_REF] Bollobas | Random Graphs[END_REF].

[76] Schöning's original work was devoted to the analysis of RWSAT on K-SAT, but his result holds for K-XORSAT too.

[77] A proof of this statement was obtained by [START_REF] Alekhnovich | Analysis of the Random Walk Algorithm on Random 3-CNFs[END_REF] for the random SAT model.

[78] This decrease characterises the overall operation of RWSAT. A precise look at the e(T ) curve reveals that the energy e may occasionally increase.

[79] On intuitive grounds, as a step of the algorithm can satisfy θ(1) equations at a time, we expect the average value of Tres to be of the order of the number M of equations at least. Thus tres should grow at least linearly with α. Experiments shows that the growth is in fact more than linear.

[80] The median is more representative of the typical value of the number of splits than the expectation value, since the latter may be dominated by huge and unlikely samples, see discussion of Section II B.

[81] Equations with a single variable are created too, but are eliminated through unit-propagation. When a heuristic assignment has to be made the system is a mixture of equations with 2 and 3 variables only.

[82] Analysis of more sophisticated heuristics e.g. based on the number of occurences of variables require to handle more complex instance distributions [START_REF] Kaporis | How to prove conditional randomness using the principle of deferred decisions[END_REF].

[83] A small amount of backtracking may be necessary to find the solution since ρsuccess < 1 [START_REF] Frieze | [END_REF], but the overall picture of a single branch is not qualitatively affected.

[84] Non standard behaviour e.g. fluctuations of the order of N ν with ν = 1 2 as found in Levy flights correspond to nonanalyticies of ω in h * or the vanishing of the second derivative.

[85] Consider the following experiment. You are given three coins: the first one is fair (coin A), the second and third coins, respectively denoted by B and C, are biased and give head with probabilities, respectively,