text:
# Approximate analysis of search algorithms with “physical” methods

## I Introduction

The computational effort needed to deal with large combinatorial structures considerably varies with the task to be performed and the resolution procedure usedpapadimi . The worst case complexity of a task, more precisely an optimization or decision problem, is defined as the time required by the best algorithm to treat any possible inputs to the problem. For instance, the sorting problem of a list of $`N`$ numbers has worst-case complexity $`N\mathrm{log}N`$: there exists several algorithms that can order any list in at most $`N\mathrm{log}N`$ elementary operations, and none with asymptotically less operations. Unfortunately, the worst-case complexities of many important computational problems, called NP-Complete, is not known. Partitioning a list of $`N`$ numbers in two sets with equal partial sums is one among hundreds of such NP-complete problems. It is a fundamental conjecture of theoretical computer science that there exists no algorithm capable of partitioning any list of length $`N`$, or of solving any other NP-Complete problem with inputs of size $`N`$, in a time bounded by a polynomial of $`N`$. Therefore, when dealing with such a problem, one necessarily uses algorithms which may takes exponential times on some inputs. Quantifying how ‘frequent’ these hard inputs are for a given algorithm is the question answered by the analysis of algorithms. In this paper, we will present an overview of recent works done by physicists to address this point, and more precisely to characterize the average performances, called hereafter complexity, of a given algorithm over a distribution of inputs to an optimization problem.

The history of algorithm analysis by physical methods/ideas is at least as old as the use of computers by physicists. One well-established chapter in this history is, for instance, the analysis of Monte Carlo sampling algorithms for statistical mechanics models. In this context, it is well known that phase transitions, i.e. abrupt changes in the physical properties of the model, can imply a dramatic increase in the time necessary to the sampling procedure. This phenomenon is commonly known as critical slowing down. The physicists’ insight in this problem comes mainly from the analogy between the dynamics of algorithms and the physical dynamics of the system. This analogy is quite natural: in fact many algorithms mimick the physical dynamics itself.

A quite new idea is instead to abstract from physically motivated problems and use statistical mechanics ideas for analyzing the dynamics of algorithms. In effect there are many reasons which suggest that analysis of algorithms and statistical physics should be considered close relatives. In both cases one would like to understand the asymptotic behavior of dynamical processes acting on exponentially large (in the size of the problem) configuration spaces. The differences between the two disciplines mainly lie in the methods (and, we are tempted to say, the style) of investigation. Theoretical computer science derives rigorous results based on probability theory. However these results are sometimes too weak for a complete characterization of the algorithm. Physicists provide instead heuristic results based on intuitively sensible approximations. These approximations are eventually validated by a comparison with numerical experiments. In some lucky cases, approximations are asymptotically irrelevant: estimates are turned into conjectures left for future rigorous derivations.

Perhaps more interesting than stylistic differences is the point of view which physics brings with itself. Let us highlight two consequences of this point of view.

First, a particular importance is attributed to “complexity phase transitions” i.e. abrupt changes in the resolution complexity as some parameter defining the input distribution is variedAI ; Friedgut . We shall consider two examples in the next Sections:

* Random Satisfiability of Boolean constraints (SAT). In $`K`$-SAT one is given an instance, that is, a set of $`M`$ logical constraints (clauses) among $`N`$ boolean variables, and wants to find a truth assignment for the variables which fulfill all the constraints. Each clause is the logical OR of $`K`$ literals, a literal being one of the $`N`$ variables or its negation e.g. $`(x_1x_{17}\overline{x_{31}})`$ for 3-SAT. Random $`K`$-SAT is the $`K`$-SAT problem supplied with a distribution of inputs uniform over all instances having fixed values of $`N`$ and $`M`$. The limit of interest is $`N,M\mathrm{}`$ at fixed ratio $`\alpha =M/N`$ of clauses per variableMit ; Hans .
* Vertex cover of random graphs (VC). An input instance of the VC decision problem consists in a graph $`G`$ and an integer number $`X`$. The problem consists in finding a way to distribute $`X`$ covering marks over the vertices in such a way that every edge of the graph is covered, that is, has at least one of its ending vertices marked. A possible distribution of inputs is provided by drawing random graphs $`G`$ à la Erdös-Renyì i.e. with uniform probability among all the graphs having $`N`$ vertices and $`E`$ edges. The limit of interest is $`N,E\mathrm{}`$ at fixed ratio $`c=2E/N`$ of edges per vertex.

The algorithms for random SAT and VC we shall consider in the next Sections undergo a complexity phase transition as the input parameter $`\pi `$ ($`=\alpha `$ for SAT, $`c`$ for VC) crosses some critical threshold $`\pi _{\mathrm{alg}}`$. Typically resolution of a randomly drawn instance requires linear time below the threshold $`\pi <\pi _{\mathrm{alg}}`$ and exponential time above $`\pi >\pi _{\mathrm{alg}}`$. The observation that most difficult instances are located near the phase boundary confirms the relevance of the phase-transition phenomenon.

Secondly, a key role is played by the intrinsic (algorithm independent) properties of the instance under study. The intuition is that, underlying the dramatic slowing down of a particular algorithm, there can be some qualitative change in some structural property of the problem e.g. the geometry of the space of solutions. While there is no general understanding of this question, we can further specify the above statements case-by-case. Let us consider, for instance, a local search algorithm for a combinatorial optimization problem. If the algorithm never increases the value of the cost function $`F(C)`$ where $`C`$ is the configuration (assignment) of variables to be optimized over, the number and geometry of the local minima of $`F(C)`$ will be crucial for the understanding of the dynamics of the algorithm. This example is illustrated in Sec. III.3. While the “dynamical” behavior of a particular algorithm is not necessarily related to any “static” property of the instance, this approach is nevertheless of great interest because it could provide us with some ‘universal’ results. Some properties of the instance, for example, may imply the ineffectiveness of an entire class of algorithms.

While we shall mainly study in this paper the performances of search algorithms applied to hard combinatorial problems as SAT, VC, we will also consider easy, that is, polynomial problems as benchmarks for these algorithms. The reason is that we want to understand if the average hardness of resolution of solving NP-complete problems with a given distribution of instances and a given algorithm truly reflects the intrinsic hardness of these combinatorial problems or is simply due to some lack of efficiency of the algorithm under study. The benchmark problem we shall consider is random XORSAT. It is a version of a satisfiability problem, much simpler than SAT from a computational complexity point of viewCrei . The only but essential difference with SAT is that a clause is said to be satisfied if the exclusive, and not inclusive, disjunction of its literals is true. XORSAT may be recast as a linear algebra problem, where a set of $`M`$ equations involving $`N`$ Boolean variables must be satisfied modulo 2, and is therefore solvable in polynomial time by various methods e.g. Gaussian elimination. Nevertheless, it is legitimate to ask what are the performances of general search algorithms for this kind of polynomial computational problem. In particular, we shall see that some algorithms requiring exponential times to solve random SAT instances behave badly on random XORSAT instances too. A related question we shall focus on in Sec. III.2 is decoding, which may also, in some cases, be expressed as the resolution of a set of Boolean equations.

The paper is organized as follows. In Sec. II.1 we shall review backtracking search algorithms, which, roughly speaking, work in the space of instances. We explain the general ideas and then illustrate them on random SAT (Sec. II.2) and VC (Sec. II.3). In Sec. II.4 we consider the fluctuations in running times of these algorithms and analyze the possibility of exploiting these fluctuations in random restart strategies. In Sec. III we turn to local search algorithms, which work in the space of configurations. We review the analysis of such algorithms for decoding problems (Sec. III.2), random XORSAT (Sec. III.3), and SAT (Sec. III.4). Finally in the Conclusion we suggest some possible future developments in the field.

## II Analysis of the Davis-Putnam-Loveland-Logeman search procedure

### II.1 Overview of the algorithm and physical concepts

In this section, we briefly review the Davis-Putnam-Loveland-Logemann (DPLL) procedure DP ; survey . A decision problem can be formulated as a constrained satisfaction problem, where a set of variables must be sought for to fulfill some given constraints. For simplicity, we suppose here that variables may take a finite set of values with cardinality $`v`$ e.g. $`v=2`$ for SAT or VC. DPLL is an exhaustive search procedure operating by trials and errors, the sequence of which can be graphically represented by a search tree (Fig. 1). The tree is defined as follows: (1) A node in the tree corresponds to a choice of a variable. (2) An outgoing branch (edge) codes for the value of the variable and the logical implications of this choice upon not yet assigned variables and clauses. Obviously a node gives birth to $`v`$ branches at most. (3) Implications can lead to: (3.1) a violated constraint, then the branch ends with $`C`$ (contradiction), the last choice is modified (backtracking of the tree) and the procedure goes on along a new branch (point 2 above); (3.2) a solution when all constraints are satisfied, then the search process is over; (3.3) otherwise, some constraints remain and further assumptions on the variables have to be done (loop back to point 1).

A computer independent measure of computational complexity, that is, the amount of operations necessary to solve the instance, is given by the size $`Q`$ of the search tree i.e. the number of nodes it contains. Performances can be improved by designing sophisticated heuristic rules for choosing variables (point 1). The resolution time (or complexity) is a stochastic variable depending on the instance under consideration and on the choices done by the variable assignment procedure. Its average value, $`\overline{Q}`$, is a function of the input distribution parameters $`\pi `$ e.g. the ratio $`\alpha `$ of clauses per variable for SAT, or the average degree $`c`$ for the VC of random graphs, which can be measured experimentally and that we want to calculate theoretically. More precisely, our aim is to determine the values of the input parameters for which the complexity is linear, $`\overline{Q}=\gamma N`$ or exponential, $`\overline{Q}=2^{N\omega }`$, in the size $`N`$ of the instance and to calculate the coefficients $`\gamma ,\omega `$ as functions of $`\pi `$.

The DPLL algorithm gives rise to a dynamical process. Indeed, the initial instance is modified during the search through the assignment of some variables and the simplification of the constraints that contain these variables. Therefore, the parameters of the input distribution are modified as the algorithm runs. This dynamical process has been rigorously studied and understood in the case of a search tree reducing to one branch (tree A in Figure 1)fra2 ; Fra ; Achl ; Fri ; Kir ; kir2 . Study of trees with massive backtracking e.g. trees B and C in Fig. 1 is much more difficult. Backtracking introduces strong correlations between nodes visited by DPLL at very different times, but close in the tree. In addition, the process is non Markovian since instances attached to each node are memorized to allow the search to resume after a backtracking step.

The study of the operation of DPLL is based on the following, elementary observation. Since instances are modified when treated by DPLL, description of their statistical properties generally requires additional parameters with respects to the defining parameters $`\pi `$ of the input distribution. Our task therefore consists in

1. identifying these extra parameters $`\pi ^{}`$kir2 ;
2. deriving the phase diagram of this new, extended distribution $`\pi ,\pi ^{}`$ to identify, in the $`\pi ,\pi ^{}`$ space, the critical surface separating instances having solution with high probability (satisfiable phase) from instances having generally no solution (unsatisfiable phase), see Fig. 2.
3. tracking the evolution of an instance under resolution with time $`t`$ (number of steps of the algorithm), that is, the trajectory of its characteristic parameters $`\pi (t),\pi ^{}(t)`$ in the phase diagram.

Whether this trajectory remains confined to one of the two phases or crosses the boundary inbetween has dramatic consequences on the resolution complexity. We find three average behaviours, schematized on Fig. 2:

* if the initial instance has a solution and the trajectory remains in the sat phase, resolution is typically linear, and almost no backtracking is present (Fig. 1A). The coordinates of the trajectory $`\pi (t),\pi ^{}(t)`$ of the instance in the course of the resolution obey a set of coupled ordinary differential equations accounting for the changes in the distribution parameters done by DPLL.
* if the initial instance has no solution, solving the instance, that is, finding a proof of unsatisfiability, takes exponentially large time and makes use of massive backtracking (Fig. 1B). Analysis of the search tree is much more complicated than in the linear regime, and requires a partial differential equation that gives information on the population of branches with parameters $`\pi ,\pi ^{}`$ throughout the growth of the search tree.
* in some intermediary regime, instances have solutions but finding one requires an exponentially large time (Fig. 1C). This may be related to the crossing of the boundary between sat and unsat phases of the instance trajectory. We have therefore a mixed behaviour which can be understood through combination of the two above cases.

We now explain how to apply concretly this approach to the cases of random SAT and VC.

### II.2 Average analysis of the Random SAT problem

The input distribution of 3-SAT is characterized by a single parameter $`\pi `$, the ratio $`\alpha `$ of clauses per variable. The action of DPLL on an instance of 3-SAT, illustrated in Fig. 3, causes the changes of the overall numbers of variables and clauses, and thus of $`\alpha `$. Furthermore, DPLL reduces some 3-clauses to 2-clauses. We use a mixed 2+p-SAT distributionSta , where $`p(=\pi ^{})`$ is the fraction of 3-clauses, to model what remains of the input instance at a node of the search tree. Using experiments and methods from statistical mechanicsSta  and rigorous calculationsAchl1 , the threshold line $`\alpha _C(p)`$, separating sat from unsat phases, may be estimated with the results shown in Fig. 4. For $`pp_0=2/5`$, i.e. left to point T, the threshold line is given by $`\alpha _C(p)=1/(1p)`$, and saturates the upper bound for the satisfaction of 2-clauses. Above $`p_0`$, no exact value for $`\alpha _C(p)`$ is known. The phase diagram of 2+p-SAT is the natural space in which the DPLL dynamics takes place. An input 3-SAT instance with ratio $`\alpha `$ shows up on the right vertical boundary of Fig. 4 as a point of coordinates $`(p=1,\alpha )`$. Under the action of DPLL, the representative point moves aside from the 3-SAT axis and follows a trajectory in the $`(p,\alpha )`$ plane.

In this section, we show that the location of this trajectory in the phase diagram allows a precise understanding of the search tree structure and of complexity as a function of the ratio $`\alpha `$ of the instance to be solved (Inset of Fig. 4). In addition, we shall present an approximate calculation of trajectories accounting for the case of massive backtracking, that is for unsat instances, and slightly below the threshold in the sat phase. Our approach is based on a non rigorous extension of works by Chao and Franco who first studied the action of DPLL (without backtracking) on easy, sat instancesfra2 ; Fra  as a way to obtain lower bounds to the threshold $`\alpha _C`$, see Achl  for a recent review.

Let us emphasize that the idea of trajectory is made possible thanks to an important statistical property of the heuristics of split we consider fra2 ; Fra ,

* Unit-Clause (UC) heuristic: pick up randomly a literal among a unit clause if any, or any unset variable otherwise.
* Generalized Unit-Clause (GUC) heuristic: pick up randomly a literal among the shortest avalaible clauses.
* Short Clause With Majority (SC1) heuristic: pick up randomly a literal among unit clauses if any, or pick up randomly an unset variable $`v`$, count the numbers of occurences $`\mathrm{},\overline{\mathrm{}}`$ of $`v`$, $`\overline{v}`$ in 3-clauses, and choose $`v`$ (respectively $`\overline{v}`$) if $`\mathrm{}>\overline{\mathrm{}}`$ (resp. $`\mathrm{}<\overline{\mathrm{}}`$). When $`\mathrm{}=\overline{\mathrm{}}`$, $`v`$ and $`\overline{v}`$ are equally likely to be chosen.

These heuristics do not induce any bias nor correlation in the instances distributionfra2 ; kir2 . Such a statistical “invariance” is required to ensure that the dynamical evolution generated by DPLL remains confined to the phase diagram of Fig. 4. In the following, the initial ratio of clauses per variable of the instance to be solved will be denoted by $`\alpha _0`$.

#### II.2.1 Lower sat phase and branch trajectories.

Let us consider the first descent of the algorithm, that is the action of DPLL in the absence of backtracking. The search tree is a single branch (Fig. 1A). The numbers of 2 and 3-clauses are initially equal to $`C_2=0,C_3=\alpha _0N`$ respectively. Under the action of DPLL, $`C_2`$ and $`C_3`$ follow a Markovian stochastic evolution process, as the depth $`T`$ along the branch (number of assigned variables) increases. Both $`C_2`$ and $`C_3`$ are concentrated around their average values, the densities $`c_j(t)=E[C_j(tN)/N]`$ ($`j=2,3`$) of which obey a set of coupled ordinary differential equations (ODE)fra2 ; Fra ; Achl ,

$$\frac{dc_3}{dt}=\frac{3c_3}{1t},\frac{dc_2}{dt}=\frac{3c_3}{2(1t)}\frac{2c_2}{1t}\rho _1(t)h(t),$$
(1)
where $`\rho _1(t)=1c_2(t)/(1t)`$ is the probability that DPLL fixes a variable at depth $`t`$ through unit-propagation. Function $`h`$ depends upon the heuristic: $`h_{UC}(t)=0`$, $`h_{GUC}(t)=1`$ (if $`\alpha _0>2/3`$), $`h_{SC_1}(t)=ae^a(I_0(a)+I_1(a))/2`$ with $`a3c_3(t)/(1t)`$ and $`I_{\mathrm{}}`$ denotes the $`\mathrm{}^{th}`$ modified Bessel function. To obtain the single branch trajectory in the phase diagram of Fig. 4, we solve the ODEs (1) with initial conditions $`c_2(0)=0,c_3(0)=\alpha _0`$, and perform the change of variables

$$p(t)=\frac{c_3(t)}{c_2(t)+c_3(t)},\alpha (t)=\frac{c_2(t)+c_3(t)}{1t}.$$
(2)
Results are shown for the GUC heuristics and starting ratios $`\alpha _0=2`$ and 2.8 in Fig. 4. Trajectories, indicated by light dashed lines, first head to the left and then reverse to the right until reaching a point on the 3-SAT axis at a small ratio. Further action of DPLL leads to a rapid elimination of the remaining clauses and the trajectory ends up at the right lower corner S, where a solution is found.

Frieze and Suen Fri  have shown that, for ratios $`\alpha _0<\alpha _L3.003`$ (for the GUC heuristics), the full search tree essentially reduces to a single branch, and is thus entirely described by the ODEs (1). The number of backtrackings necessary to reach a solution is bounded from above by a power of $`\mathrm{log}N`$. The average size $`\overline{Q}`$ of the branch then scales linearly with $`N`$ with a multiplicative factor $`\gamma (\alpha _0)=Q/N`$ that can be analytically computed Coc .

The boundary $`\alpha _L`$ of this easy sat region can be defined as the largest initial ratio $`\alpha _0`$ such that the branch trajectory $`p(t),\alpha (t)`$ issued from $`\alpha _0`$ never leaves the sat phase in the course of DPLL resolution.

#### II.2.2 Unsat phase and tree trajectories.

For ratios above threshold ($`\alpha _0>\alpha _C4.3`$), instances almost never have a solution but a considerable amount of backtracking is necessary before proving that clauses are incompatible. As shown in Fig. 1B, a generic unsat tree includes many branches. The number of branches (leaves), $`B`$, or the number of nodes, $`Q=B1`$, grow exponentially with $`N`$Chv . It is convenient to define its logarithm $`\omega `$ through $`B=2^{N\omega }`$. Contrary to the previous section, the sequence of points $`(p,\alpha )`$ characterizing the evolution of the 2+p-SAT instance solved by DPLL does not define a line any longer, but rather a patch, or cloud of points with a finite extension in the phase diagram of Fig. 2.

We have analytically computed the logarithm $`\omega `$ of the size of these patches, as a function of $`\alpha _0`$, extending to the unsat region the probabilistic analysis of DPLL. This is, a priori, a very difficult task since the search tree of Fig. 1B is the output of a complex, sequential process: nodes and edges are added by DPLL through successive descents and backtrackings. We have imagined a different building up, that results in the same complete tree but can be mathematically analyzed: the tree grows in parallel, layer after layer (Fig. 5).

A new layer is added by assigning, according to DPLL heuristic, one more variable along each living branch. As a result, a branch may split (case 1), keep growing (case 2) or carry a contradiction and die out (case 3). Cases 1,2 and 3 are stochastic events, the probabilities of which depend on the characteristic parameters $`c_2,c_3,t`$ defining the 2+p-SAT instance carried by the branch, and on the depth (fraction of assigned variables) $`t`$ in the tree. We have taken into account the correlations between the parameters $`c_2,c_3`$ on each of the two branches issued from splitting (case 1), but have neglected any further correlation which appear between different branches at different levels in the treeCoc . This Markovian approximation permits to write an evolution equation for the logarithm $`\omega (c_2,c_3,t)`$ of the average number of branches with parameters $`c_2,c_3`$ as the depth $`t`$ increases,

$$\frac{\omega }{t}(c_2,c_3,t)=H[c_2,c_3,\frac{\omega }{c_2},\frac{\omega }{c_3},t].$$
(3)
$`H`$ incorporates the details of the splitting heuristics. In terms of the partial derivatives $`y_2=\omega /c_2`$, $`y_3=\omega /c_3`$, we have for the UC and GUC heuristics

$`H_{UC}`$ $`=`$ $`1+{\displaystyle \frac{1}{\mathrm{ln}2}}\left[{\displaystyle \frac{3c_3}{1t}}\left(e^{y_3}{\displaystyle \frac{1+e^{y_2}}{2}}1\right)+{\displaystyle \frac{c_2}{1t}}\left({\displaystyle \frac{3}{2}}e^{y_2}2\right)\right]`$
$`H_{GUC}`$ $`=`$ $`\mathrm{log}_2\nu (y_2)+{\displaystyle \frac{1}{\mathrm{ln}2}}\left[{\displaystyle \frac{3c_3}{1t}}\left(e^{y_3}{\displaystyle \frac{1+e^{y_2}}{2}}1\right)+{\displaystyle \frac{c_2}{1t}}\left(\nu (y_2)2\right)\right]`$
where $`\nu (y_2)={\displaystyle \frac{1}{2}}e^{y_2}\left(1+\sqrt{1+4e^{y_2}}\right).`$ (4)
Partial differential equation (PDE) (3) is analogous to growth processes encountered in statistical physics Gro . The surface $`\omega `$, growing with “time” $`t`$ above the plane $`(c_2,c_3)`$, or equivalently from (2), above the plane $`(p,\alpha )`$ (Fig. 6), describes the whole distribution of branches. The average number of branches at depth $`t`$ in the tree equals $`B(t)=𝑑p𝑑\alpha \mathrm{\hspace{0.33em}2}^{N\omega (p,\alpha ,t)}2^{N\omega ^{}(t)}`$, where $`\omega ^{}(t)`$ is the maximum over $`p,\alpha `$ of $`\omega (p,\alpha ,t)`$ reached in $`p^{}(t),\alpha ^{}(t)`$. In other words, the exponentially dominant contribution to $`B(t)`$ comes from branches carrying 2+p-SAT instances with parameters $`p^{}(t),\alpha ^{}(t)`$, which define the tree trajectories on Fig. 4.

The hyperbolic line in Fig. 4 indicates the halt points, where contradictions prevent dominant branches from further growing. Each time DPLL assigns a variable through unit-propagation, an average number $`u(p,\alpha )`$ of new 1-clauses is produced, resulting in a net rate of $`u1`$ additional 1-clauses. As long as $`u<1`$, 1-clauses are quickly eliminated and do not accumulate. Conversely, if $`u>1`$, 1-clauses tend to accumulate. Opposite 1-clauses $`x`$ and $`\overline{x}`$ are likely to appear, leading to a contradiction Fra ; Fri . The halt line is defined through $`u(p,\alpha )=1`$. As far as dominant branches are concerned, the equation of the halt line reads

$$\alpha =\left(\frac{3+\sqrt{5}}{2}\right)\mathrm{ln}\left[\frac{1+\sqrt{5}}{2}\right]\frac{1}{1p}\frac{1.256}{1p}.$$
(5)
Along the tree trajectory, $`\omega ^{}(t)`$ grows from 0, on the right vertical axis, up to some final positive value, $`\widehat{\omega }`$, on the halt line. $`\widehat{\omega }`$ is our theoretical prediction for the logarithm of the complexity (divided by $`N`$). Values of $`\widehat{\omega }`$ obtained for $`4.3<\alpha _0<20`$ by solving equation (3) compare very well with numerical results Coc .

We have plotted the surface $`\omega `$ above the $`(p,\alpha )`$ plane, with the results shown in Fig. 6. It must be stressed that, though our calculation is not rigorous, it provides a very good quantitative estimate of the complexity. Furthermore, complexity is found to scale asymptotically as

$$\widehat{\omega }(\alpha _0)\frac{3+\sqrt{5}}{(6\mathrm{ln}2)\alpha _0}\left[\mathrm{ln}\left(\frac{1+\sqrt{5}}{2}\right)\right]^2\frac{0.292}{\alpha _0}(\alpha _0\alpha _C).$$
(6)
This result exhibits the expected scalingBea , and could indeed be exact. As $`\alpha _0`$ increases, search trees become smaller and smaller, and correlations between branches, weaker and weaker.

#### II.2.3 Upper sat phase and mixed branch–tree trajectories.

The interest of the trajectory approach proposed in this paper is best seen in the upper sat phase, that is ratios $`\alpha _0`$ ranging from $`\alpha _L`$ to $`\alpha _C`$. This intermediate region juxtaposes branch and tree behaviors, see Fig. 1C. The branch trajectory starts from the point $`(p=1,\alpha _0)`$ corresponding to the initial 3-SAT instance and hits the critical line $`\alpha _c(p)`$ at some point G with coordinates ($`p_G,\alpha _G`$) after $`Nt_G`$ variables have been assigned by DPLL (Fig. 4). The algorithm then enters the unsat phase and generates 2+p-SAT instances with no solution. A dense subtree, that DPLL has to go through entirely, forms beyond G till the halt line (Fig. 4). The size of this subtree, $`2^{N(1t_G)\widehat{\omega }_G}`$, can be analytically predicted from our theory. G is the highest backtracking node in the tree (Fig. 1C) reached back by DPLL, since nodes above G are located in the sat phase and carry 2+p-SAT instances with solutions. DPLL will eventually reach a solution. The corresponding branch (rightmost path in Fig. 1C) is highly non typical and does not contribute to the complexity, since almost all branches in the search tree are described by the tree trajectory issued from G (Fig. 4). We have checked experimentally this scenario for $`\alpha _0=3.5`$. The coordinates of the average highest backtracking node, $`(p_G0.78,\alpha _G3.02`$), coincide with the analytically computed intersection of the single branch trajectory and the critical line $`\alpha _c(p)`$Coc . As for complexity, experimental measures of $`\omega `$ from 3-SAT instances at $`\alpha _0=3.5`$, and of $`\omega _G`$ from 2+0.78-SAT instances at $`\alpha _G=3.02`$, obey the expected identity $`\omega =\omega _G(1t_G)`$ and are in very good agreement with theoryCoc . Therefore, the structure of search trees for 3-SAT reflects the existence of a critical line for 2+p-SAT instances.

### II.3 Average analysis of the vertex cover of random graphs

We now consider the VC problem, where inputs are random graphs drawn from the $`G(N,p=c/N)`$ ensembleBollo . In other words, graphs have $`N`$ vertices and the probability that a pair of vertices are linked through an edge is $`c/N`$, independently of other edges. When the number $`X=xN`$ of covering marks is lowered, the model undergoes a COV/UNCOV transition at some critical density of covers $`x_\mathrm{c}(c)`$ for $`N\mathrm{}`$. For $`x>x_\mathrm{c}(c)`$, vertex covers of size $`Nx`$ exist with probability one, for $`x<x_\mathrm{c}(c)`$ the available covering marks are not sufficient. The statistical mechanics analysis of Ref. WeHa1  gave the result

$`x_\mathrm{c}(c)=1{\displaystyle \frac{2W(c)+W(c)^2}{2c}},\text{for}c<e,`$ (7)
where $`W(c)`$ solves the equation $`We^W=c`$. This result is compatible with the bounds of Refs. Ga\_VC ; Fr\_VC , and was later shown to be exact Bauer . For $`c>e`$, Eq. (7) only gives an approximate estimate of $`x_\mathrm{c}(c)`$. More sophisticated calculations can be found in Ref. WeHaLong .

Let us consider a simple implementation of the DPLL procedure for the present problem. During the computation, vertices can be covered, uncovered or just free, meaning that the algorithm has not yet assigned any value to that vertex. At the beginning all the vertices are set free. At each step the algorithm chooses a vertex $`i`$ at random among those which are free. If $`i`$ has neighboring vertices which are either free or uncovered, then the vertex $`i`$ is declared covered first. In case $`i`$ has only covered neighbors, the vertex is declared uncovered. The process continues unless the number of covered vertices exceeds $`X`$. In this case the algorithm backtracks and the opposite choice is taken for the vertex $`i`$ unless this corresponds to declaring uncovered a vertex that has one or more uncovered neighbors. The algorithm halts if it finds a solution (and declares the graph to be COV) or after exploring all the search tree (in this case it declares the graph to be UNCOV).

Of course one can improve over this algorithm by using smarter heuristics WeHeur . One remarkable example is the “leaf-removal” algorithm defined in Ref. Bauer . Instead of picking any vertex randomly, one chooses a connectivity-one vertex, declare it uncovered, and declare covered its neighbor. This procedure is repeated iteratively on the subgraph of free nodes, until no connectivity-one nodes are left. In the low-connectivity, COV region $`\{c<e,x>x_\mathrm{c}(c)\}`$, it stops only when the graph is completely covered. As a consequence, this algorithm can solve VC in linear time with high probability in all this region. No equally good heuristics exists for higher connectivity, $`c>e`$.

#### II.3.1 Branch trajectories

Under the action of one of the above algorithms, the instance is progressively modified and the number of variables is reduced. In fact, at each step a vertex is selected and can be eliminated from the graph regardless whether it is declared covered or uncovered. The analysis of the first algorithm is greatly simplified by the remark that, as long as backtracking has not begun, the new vertex is selected randomly. This implies that the modified instance produced by the algorithm is still a random graph. Its evolution can be effectively described by a trajectory in the $`(c,x)`$ space. If one starts from the parameters $`c_0`$, $`x_0`$, after $`Nt`$ steps of the algorithm, he will end up with a new instance of size $`N(1t)`$ and parameters WeHa2

$`c(t)=c_0(1t),x(t)={\displaystyle \frac{x_0t}{1t}}+{\displaystyle \frac{e^{c_0(1t)}e^{c_0}}{c_0(1t)}}.`$ (8)
Some examples of the two types of trajectories (the ones leading to a solution and the ones which eventually enter the UNCOV region) are shown in Fig. 7. The separatrix is given by

$`x_\mathrm{s}(c)=1{\displaystyle \frac{1e^c}{c}},`$ (9)
and corresponds to the dotted line in Fig. 7. Above this line the algorithm solves the problem in linear time.

For more general heuristics the analysis becomes less straightforward because the graph produced by the algorithm does not belong to the standard random-graph ensemble. It may be necessary to augment the number of parameters which describe the evolution of the instance. As an example, the leaf-removal algorithm mentioned in the previous Section is conveniently described by keeping track of three numbers which parametrize the degree profile (i.e. the fraction of vertices $`p_d(t)`$ having a given degree $`d`$) of the graph WeHeur .

#### II.3.2 Tree trajectories

Below the critical line $`x_\mathrm{c}(c)`$, cf. Eq. (7), no solution exists to the typical random instance of VC. Our algorithm must explore a large backtracking tree to prove it and this takes an exponential time. The size of the backtracking tree could be computed along the lines of Sec. II.B.2. However a good result can be obtained with a simple “static” calculation WeHa1 .

As explained in Sec. II.B.2, we imagine the evolution of the backtracking tree as proceeding “in parallel”. At the level $`M`$ of the tree a set of $`M`$ vertices has been visited. Call $`𝒢_M`$ the subgraph induced by these vertices. Since we always put a covering mark on a vertex which is surrounded by vertices declared uncovered, each node of the backtracking tree will carry a vertex cover of the associated subgraph $`𝒢_M`$ . Therefore the number of backtracking nodes is given by

$`Q={\displaystyle \underset{M=1}{\overset{N}{}}}𝒩_{\mathrm{VC}}(𝒢_M;X),`$ (10)
where $`𝒩_{\mathrm{VC}}(𝒢_M;X)`$ is the number of VC’s of $`𝒢_M`$ using at most $`X`$ marks. A very crude estimate of the right-hand side of the above equation is:

$`Q{\displaystyle \underset{M=1}{\overset{N}{}}}{\displaystyle \underset{X^{}=0}{\overset{\mathrm{min}(X,M)}{}}}\left(\begin{array}{c}M\\ X^{}\end{array}\right),`$ (13)
where we bounded the number of VC’s of size $`X^{}`$ on $`𝒢_M`$ with the number of ways of placing $`X^{}`$ marks on $`M`$ vertices. The authors of WeHa2  provided a refined estimate based on the annealed approximation of statistical mechanics. The results of this calculation are compared in Fig. 8 with the numerics.

#### II.3.3 Mixed trajectories

If the parameters which characterize an instance of VC lie in the region between $`x_\mathrm{c}(c)`$, cf. Eq. (7), and $`x_\mathrm{s}(c)`$, cf. Eq. (9), the problem is still soluble but our algorithm takes an exponential time to solve it. In practice, after a certain number of vertices has been visited and declared either covered or uncovered, the remaining subgraph $`𝒢_{free}`$ cannot be any longer covered with the leftover marks. This happens typically when the first descent trajectory (8) crosses the critical line (7).

It takes some time for the algorithm to realize this fact. More precisely, it takes exactly the time necessary to prove that $`𝒢_{free}`$ is uncoverable. This time dominates the computational complexity in this region and can be calculated along the lines sketched in the previous Section. The result is, once again, reported in Fig. 8, which clearly shows a computational peak at the phase boundary.

Finally, let us notice that this mixed behavior disappears in the entire $`c<e`$ region if the leaf-removal heuristics is adopted for the first descent.

### II.4 Distribution of resolution times

Up to now we have studied the typical resolution complexity. The study of fluctuations of resolution times is interesting too, particularly in the upper sat phase where solutions exist but are found at a price of a large computational effort. We may expect that there exist lucky but rare resolutions able to find a solution in a time much smaller than the typical one. Due to the stochastic character of DPLL complexity indeed fluctuates from run to run of the algorithm on the same instance. In Fig. 9 we show this run-to-run distribution of the logarithm $`\omega `$ of the resolution complexity for four instances of random 3-SAT with the same ratio $`\alpha =3.5`$. The run to run distribution are qualitatively independent of the particular instances, and exhibit two bumps. The wide right one, located in $`\omega 0.035`$, correspond to the major part of resolutions. It acquires more and more weight as $`N`$ increases and corresponds to the typical behavior analysed in Section II.B.3. The left peak corresponds to much faster resolutions, taking place in linear time. The weight of this peak (fraction of runs with complexities falling in the peak) decreases exponentially fast with $`N`$, and can be numerically estimated to $`W_{lin}=2^{N\zeta }`$ with $`\zeta 0.011`$. Therefore, instances at $`\alpha =3.5`$ are typically solved in exponential time but a tiny (exponentially small) fraction of runs are able to find a solution in linear time only.

A systematic stop-and-restart procedure may be introduced to take advantage of this fluctuation phenomenon and speed up resolution. If a solution is not found before $`N`$ splits, DPLL is stopped and rerun after some random permutations of the variables and clauses. The expected number $`N_{rest}`$ of restarts necessary to find a solution being equal to the inverse probability $`1/W_{lin}`$ of linear resolutions, the resulting complexity scales as $`NW_{lin}^12^{N\zeta }`$.

To calculate $`\zeta `$ we have analyzed, along the lines of the study of the growth of the search tree in the unsat phase, the whole distribution of the complexity for a given ratio $`\alpha `$ in the upper sat phase. Calculations can be found in Cocrs . Linear resolutions are found to correspond to branch trajectories that cross the unsat phase without being hit by a contradiction, see Fig. 4. Results are reported in Fig. 10 and compare very well with the experimentally measured number $`N_{rest}`$ of restarts necessary to find a solution. In the whole upper sat phase, the use of restarts offers an exponential gain with respect to usual DPLL resolution (see Fig. 10 for comparison between $`\zeta `$ and $`\omega `$), but the completeness of DPLL is lost.

A slightly more general restart strategy consists in stopping the backtracking procedure after a fixed number of nodes $`Q_R=e^{N\omega _R^{}}`$ has been visited. A new (and statistically independent) DPLL procedure is then started from the beginning. In this case one exploits lucky, but still exponential, stochastic runs. The tradeoff between the exponential gain of time and the exponential number of restarts, can be optimized by tuning the parameter $`\omega _R^{}`$. This approach has been analyzed in Ref. VCRestart  taking VC as a working example. In Fig. 11 we show the computatonal complexity of such a strategy as a function of the restart parameter $`\omega _R^{}`$. We compare the numerics with an approximate calculation VCRestart . The instances were random graphs with average connectivity $`c=3.2`$, and $`x=0.6`$ covering marks per vertex. The optimal choice of the parameter seems to be (in this case) $`\omega _R^{}0`$, corresponding to polynomial runs.

The analytical prediction reported in Fig. 11 requires, as for 3-SAT, an estimate of the execution-time fluctuations of the DPLL procedure (without restart). It turns out that one major source of fluctuations is, in the present case, the location in the $`(c,x)`$ plane of the highest node in the backtracking tree. In the typical run this coincides with the intersection $`(c_G,x_G)`$ between the first descent trajectory (8) and the critical line (7). One can estimate the probability $`P(c,x)\mathrm{exp}\{N\psi (c,x)\}`$ for this node to have coordinates $`(c,x)`$ (obviously $`\psi (c_G,x_G)=0`$).

When an upper bound $`\omega _R^{}`$ on the backtracking time is fixed, the problem is solved in those lucky runs which are characterized by an atypical highest backtracking node. Roughly speaking, this means that the algorithm has made some very good (random) choices in its first steps. In Fig. 12 we plot the position of the highest backtracking point in the (last) successful runs for several values of $`\omega _R^{}`$. Once again the numerics compare favourably with an approximate calculation.

## III Analysis of local search algorithms

We now turn to the description and study of algorithms of another type, namely local search algorithms. As a common feature, these algorithms start from a configuration (assignment) of the variables, and then make successive improvements by changing at each step few of the variables in the configuration (local move). For instance, in the SAT problem, one variable is flipped from being true to false, or vice versa, at each step. Whereas complete algorithms of the DPLL type give a definitive answer to any instance of a decision problem, exhibiting either a solution or a proof of unsatisfiability, local search algorithms give a sure answer when a solution is found but cannot prove unsatisfiability. However, these algorithms can sometimes be turned into one-sided probabilistic algorithms, with an upper bound on the probability that a solution exists and has not been found after $`T`$ steps of the algorithm, decreasing to zero when $`T\mathrm{}`$rando .

### III.1 Landscape and search dynamics

Local search algorithms perform repeated changes of a configuration $`C`$ of variables (values of the Boolean variables for SAT, status –marked or unmarked– of vertices for VC) according to some criterion, usually based on the comparison of the cost function $`F`$ (number of unsatisfied clauses for SAT, of uncovered edges for VC) evaluated at $`C`$ and over its neighborhood. It is therefore clear that the shape of the multidimensional surface $`CF(C)`$, called cost function landscape, is of high importance. On intuitive grounds, if this landscape is relatively smooth with a unique minimum, local procedures as gradient descent should be very efficient, while the presence of many local minima could hinder the search process (Fig. 13). The fundamental underlying question is whether the performances of the dynamical process (ability to find the global minimum, time needed to reach it) can be understood in terms of an analysis of the cost function landscape only.

This question was intensively studied and answered for a limited class of cost functions, called mean field spin glass models, some years agoleti . The characterization of landscapes is indeed of huge importance in physical systems. There, the cost function is simply the physical energy, and local dynamics are usually low or zero temperature Monte Carlo dynamics, essentially equivalent to gradient descent. Depending on the parameters of the input distribution, the minima of the cost functions may undergo structural changes, a phenomenon called clustering in physics.

Clustering has been rigorously shown to take place in the random 3-XORSAT problemCrei ; xorsat1 ; xorsat2 , and is likely to exist in many other random combinatorial problems as 3-SATvaria ; cavi . Instances of the 3-XORSAT problem with $`M=\alpha N`$ clauses and $`N`$ variables have almost surely solutions as long as $`\alpha <\alpha _c0.918`$xorsat1 ; xorsat2 . The clustering phase transition takes place at $`\alpha _s0.818`$ and is related to a change in the geometric structure of the space of solutions, see Fig. 13:

* when $`\alpha <\alpha _s`$, the space of solutions is connected. Given a pair of solutions $`C,C^{}`$, i.e. two assignments of the $`N`$ Boolean variables that satisfy the clauses, there almost surely exists a sequence of solutions, $`C_j,j=0,1,2,\mathrm{},J`$, with $`C_0C`$, $`C_JC^{}`$, $`J=O(N)`$, connecting the two solutions such that the Hamming distance (number of different variables) between $`C_j`$ and $`C_{j+1}`$ is bounded from above by some finite constant when $`N\mathrm{}`$.
* when $`\alpha _s<\alpha <\alpha _c`$, the space of solutions is not connected any longer. It is made of an exponential (in $`N`$) number of connected components, called clusters, each containing an exponentially large number of solutions. Clusters are separated by large voids: the Hamming distance between two clusters, that is, the smallest Hammming distance between pairs of solutions belonging to these clusters, is of the order of $`N`$.

From intuitive grounds, changes of the statistical properties of the cost function landscape e.g. of the structure of the solutions space may potentially affect the search dynamics. This connection between dynamics and static properties was established in numerous works in the context of mean field models of spin glasses leti , and subsequently also put forward in some studies of local search algorithms in combinatorial optimization problemsvaria ; suedois ; cavi . So far, there is no satisfying explanation to when and why features of a priori algorithm dependent dynamical phenomena should be related to, or predictable from some statistical properties of the cost function landscape. We shall see some examples in the following where such a connection indeed exist (Sec. III.B) and other ones where its presence is far less obvious (Sec. IIIC,D).

### III.2 Algorithms for error correcting codes

Coding theory is a rich source of computational problems (and algorithms) for which the average case analysis is relevant Barg ; Spielman . Let us focus, for sake of concreteness, on the decoding problem. Codewords are sequences of symbols with some built-in redundancy. If we consider the case of linear codes on a binary alphabet, this redundancy can be implemented as a set of linear constraints. In practice, a codeword is a vector $`\underset{¯}{𝚡}\{0,1\}^N`$ (with $`N1`$) which satisfies the equation

$`\underset{¯}{𝚡}=\underset{¯}{\mathrm{𝟶}}(\mathrm{mod}\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}2}),`$ (14)
where $``$ is an $`M\times N`$ binary matrix (parity check matrix). Each one of the $`M`$ linear equations involved in Eq. (14) is called a parity check. This set of equation can be represented graphically by a Tanner graph, cf. Fig. 14. This is a bipartite graph highlighting the relations between the variables $`𝚡_i`$ and the constraints (parity checks) acting on them. The decoding problem consists in finding, among the solutions of Eq. (14), the “closest” one $`\underset{¯}{𝚡}_\mathrm{d}`$ to the output $`\underset{¯}{𝚡}_{\mathrm{out}}`$ of some communication channel. This problem is, in general, NP-hard Berlekamp .

The precise meaning of “closest” depends upon the nature of the communication channel. Let us make two examples:

* The binary symmetric channel (BSC). In this case the output of the communication channel $`\underset{¯}{𝚡}_{\mathrm{out}}`$ is a codeword, i.e. a solution of (14), in which a fraction $`p`$ of the entries has been flipped. “Closest” has to be understood in the Hamming-distance sense. $`\underset{¯}{𝚡}_\mathrm{d}`$ is the solution of Eq. (14) which minimizes the Hamming distance from $`\underset{¯}{𝚡}_{\mathrm{out}}`$.
* The binary erasure channel (BEC). The output $`\underset{¯}{𝚡}_{\mathrm{out}}`$ is a codeword in which a fraction $`p`$ of the entries has been erased. One has to find a solution $`\underset{¯}{𝚡}_\mathrm{d}`$ of Eq. (14) which is compatible with the remaining entries. Such a problem has a unique solution for small enough erasure probability $`p`$.

There are two sources of randomness in the decoding problem: $`(i)`$ the matrix $``$ which defines the code is usually drawn from some random ensemble; $`(ii)`$ the received message which is distributed according to some probabilistic model of the communication channel (in the two examples above, the bits to be flipped/erased were chosen randomly). Unlike many other combinatorial problems, there is therefore a “natural” probability distribution defined on the instances. Average case analysis with respect to this distribution is of great practical relevance.

Recently, amazingly good performances have been obtained by using low-density parity check (LDPC) codes Chung . LDPC codes are defined by parity check matrices $``$ which are large and sparse. As an example we can consider Gallager regular codes GallagerThesis . In this case $``$ is chosen with flat probability distribution within the family of matrices having $`l`$ ones per column, and $`k`$ ones per row. These are decoded using a suboptimal linear-time algorithm known as “belief-propagation” or “sum-product” algorithm GallagerThesis ; Pearl . This is an iterative algorithm which takes advantage of the locally tree-like structure of the Tanner graph, see Fig. 14, for LDPC codes. After $`n`$ iterations it incorporates the information conveyed by the variables up to distance $`n`$ from the one to be decoded. This can be done in a recursive fashion allowing for linear-time decoding.

Belief-propagation decoding shows a striking threshold phenomenon as the noise level $`p`$ crosses some critical (code-dependent) value $`p_d`$. While for $`p<p_d`$ the transmitted codeword is recovered with high probability, for $`p>p_d`$ decoding will fail almost always. The threshold noise $`p_d`$ is, in general, smaller than the threshold $`p_c`$ for optimal decoding (with unbounded computational resources).

The rigorous analysis of Ref. RichardsonUrbankeIntroduction  allows a precise determination of the critical noise $`p_d`$ under quite general circumstances. Nevertheless some important theoretical questions remain open: Can we find some smarter linear-time algorithm whose threshold is greater than $`p_d`$? Is there any “intrinsic” (i.e. algorithm independent) characterization of the threshold phenomenon taking place at $`p_d`$? As a first step towards the answer to these questions, Ref. DynamicCodes  explored the dynamics of local optimization algorithms by using statistical mechanics techniques. The interesting point is that “belief propagation” is by no means a local search algorithm.

For sake of concreteness, we shall focus on the binary erasure channel. In this case we can treat decoding as a combinatorial optimization problem within the space of bit sequences of length $`Np`$ (the number of erased bits, the others being fixed by the received message). The function to be minimized is the energy density

$`ϵ(\underset{¯}{𝚡})={\displaystyle \frac{2}{N}}d_H(\underset{¯}{𝚡},\underset{¯}{\mathrm{𝟶}}),`$ (15)
where we denote as $`d_H(\underset{¯}{𝚡}_1,\underset{¯}{𝚡}_2)`$ the Hamming distance between two vectors $`\underset{¯}{𝚡}_1`$ and $`\underset{¯}{𝚡}_2`$, and we introduced the normalizing factor for future convenience. Notice that both arguments of $`d_H(,)`$ in Eq. (15) are vectors in $`\{0,1\}^M`$.

We can define the $`R`$-neighborhood of a given sequence $`\underset{¯}{𝚡}`$ as the set of sequences $`\underset{¯}{𝚣}`$ such that $`d_\mathrm{H}(\underset{¯}{𝚡},\underset{¯}{𝚣})R`$, and we call $`R`$-stable states the bit sequences which are optima of the decoding problem within their $`R`$-neighborhood.

One can easily invent local search algorithms papadimi  for the decoding problem which use the $`R`$-neighborhoods. The algorithm start from a random sequence and, at each step, optimize it within its $`R`$-neighborhood. This algorithm is clearly suboptimal and halts on $`R`$-stable states. Let us consider, for instance, a $`(k=6,l=3)`$ regular code and decode it by local search in $`1`$-neighborhoods. In Fig. 15 we report the resulting energy density $`ϵ`$ after the local search algorithm halts, as a function of the erasure probability $`p`$. We averaged over 100 different realizations of the noise and of the matrix $``$. For sake of comparison we recall that the threshold for belief-propagation decoding is $`p_d0.429440`$ RichardsonUrbankeIntroduction , while the threshold for optimal decoding is at $`p_c0.488151`$ DynamicCodes . It is evident that local search by $`1`$-neighborhoods performs quite poorly.

A natural question is whether (and how much), these performances are improved by increasing $`R`$. It is therefore quite natural to study metastable states. These are $`R`$-stable states for any $`R=o(N)`$111We use the standard notation $`f_N=o(N)`$ if $`lim_N\mathrm{}f_N/N=0`$.. There exists no completely satisfying definition of such states: here we shall just suggest a possibility among others. The tricky point is that we do not know how to compare $`R`$-stable states for different values of $`N`$. This forbids us to make use of the above asymptotic statement. One possibility is to count without really defining them. This can be done, at least in principle, by counting $`R`$-stable states, take the $`N\mathrm{}`$ limit and, at the end, the $`R\mathrm{}`$ limitgiu . On physical grounds, we expect $`R`$-stable states to be exponentially numerous. In particular, if we call $`𝒩_R(ϵ)`$ the number of $`R`$-stable states taking a value $`ϵ`$ of the cost function (15), we have

$`𝒩_R(ϵ)\mathrm{exp}\{NS_R(ϵ)\}.`$ (16)
We can therefore define the so called (physical) complexity $`\mathrm{\Sigma }(ϵ)`$ as follows,

$`\mathrm{\Sigma }(ϵ)\underset{R\mathrm{}}{lim}S_R(ϵ).`$ (17)
Roughly speaking we can say that the number of metastable states is $`\mathrm{exp}\{N\mathrm{\Sigma }(ϵ)\}`$. Of course there are several alternative ways of taking the limits $`R\mathrm{}`$, $`N\mathrm{}`$, and we do not yet have a proof that these procedures give the same result for $`\mathrm{\Sigma }(ϵ)`$. Nevertheless it is quite clear that the existence of an exponential number of metastable states should affect dramatically the behavior of local search algorithms.

Statistical mechanics methods DynamicCodes  allows to determine the complexity $`\mathrm{\Sigma }(ϵ)`$ Remi . In “difficult” cases (such as for error-correcting codes), the actual computation may involve some approximation, e.g. the use of a variational Ansatz. Nevertheless the outcome is usually quite accurate. In Fig. 16 we consider a $`(6,3)`$ regular code on the binary erasure channel. We report the resulting complexity for three different values of the erasure probability $`p`$. The general picture is as follows. Below $`p_d`$ there is no metastable state, except the one corresponding to the correct codeword. Between $`p_d`$ and $`p_c`$ there is an exponential number of metastable states with energy density belonging to the interval $`ϵ_{GS}<ϵ<ϵ_D`$ ($`\mathrm{\Sigma }(ϵ)`$ is strictly positive in this interval). Above $`p_c`$, $`ϵ_{GS}=0`$. The maximum of $`\mathrm{\Sigma }(ϵ)`$ is always at $`ϵ_D`$.

The above picture tell us that any local algorithm will run into difficulties above $`p_d`$. In order to confirm this picture, the authors of Ref. DynamicCodes  made some numerical computations using simulated annealing as decoding algorithm for quite large codes ($`N=10^4`$ bits). For each value of $`p`$, we start the simulation fixing a fraction $`(1p)`$ of spins to $`\sigma _i=+1`$ (this part will be kept fixed all along the run). The remaining $`pN`$ spins are the dynamical variables we change during the annealing in order to try to satisfy all the parity checks. The energy of the system counts the number of unsatisfied parity checks.

The cooling schedule has been chosen in the following way: $`\tau `$ Monte Carlo sweeps (MCS) 222Each Monte Carlo sweep consists in $`N`$ proposed spin flips. Each proposed spin flip is accepted or not accordingly to a standard Metropolis test. at each of the 1000 equidistant temperatures between $`T=1`$ and $`T=0`$. The highest temperature is such that the system very rapidly equilibrates. Typical values for $`\tau `$ are from 1 to $`10^3`$.

Notice that, for any fixed cooling schedule, the computational complexity of the simulated annealing method is linear in $`N`$. Then we expect it to be affected by metastable states of energy $`ϵ_D`$, which are present for $`p>p_d`$: the energy relaxation should be strongly reduced around $`ϵ_D`$ and eventually be completely blocked. Some results are plotted in Fig. 17 together with the theoretical prediction for $`ϵ_D`$. The good agreement confirm our picture: the algorithm gets stucked in metastable states, which have, in the great majority of cases, energy density $`ϵ_D`$.

Both “belief propagation” and local search algorithms fail to decode correctly between $`p_d`$ and $`p_c`$. This leads naturally to the following conjecture: no linear time algorithm can decode in this regime of noise. The (typical case) computational complexity changes from being linear below $`p_d`$ to superlinear above $`p_d`$. In the case of the binary erasure channel it remains polynomial between $`p_d`$ and $`p_c`$ (since optimal decoding can be realized with linear algebra methods). However it is plausible that for a general channel it becomes non-polynomial.

### III.3 Gradient descent and XORSAT

In this section the local procedure we consider is gradient descent (GD). GD is defined as follows. (1) Start from an initial randomly chosen configuration of the variables. Call $`E`$ the number of unsatisfied clauses. (2) If $`E=0`$ then stop (a solution is found). Otherwise, pick randomly one variable, say $`x_i`$, and compute the number $`E^{}`$ of unsatisfied clauses when this variable is negated; if $`E^{}E`$ then accept this change i.e. replace $`x_i`$ with $`\overline{x}_i`$ and $`E`$ with $`E^{}`$; if $`E^{}>E`$, do not do anything. Then go to step 2. The study of the performances of GD to find the minima of cost functions related to statistical physics models has recently motivated various studiesgdphys ; gdproof . Numerics indicate that GD is typically able to solve random 3-SAT instances with ratios $`\alpha <3.9`$ suedois ; cavi  close to the onset of clustering varia ; cavi ; par . We shall rigorously show below that this is not so for 3-XORSAT.

Let us apply GD to an instance of XORSAT. The instance has a graph representation explained in Fig. 18. Vertices are in one–to–one correspondence with variables. A clause is fully represented by a plaquette joining three variables and a Boolean label equal to the number of negated variables it contains modulo 2 (not represented on Fig. 18). Once a configuration of the variables is chosen, each plaquette may be labelled by its status, S or U, whether the associated clause is respectively satisfied or unsatisfied. A fundamental property of XORSAT is that each time a variable is changed, i.e. its value is negated, the clauses it belongs to change status too.

This property makes the analysis of some properties of GD easy. Consider the hypergraph made of 15 vertices and 7 plaquettes in Fig. 19, and suppose the central plaquette is violated (U) while all other plaquettes are satisfied (S). The number of unsatisfied clauses is $`E=1`$. Now run GD on this special instance of XORSAT. Two cases arise, symbolized in Fig. 19, whether the vertex attached to the variable to be flipped belongs, or not, to the central plaquette. It is an easy check that, in both cases, $`E^{}=2`$ and the change is not permitted by GD. The hypergraph of Fig. 19 will be called hereafter island. When the status of the plaquettes is U for the central one and S for the other ones, the island is called blocked. Though the instance of the XORSAT problem encoded by a blocked island is obviously satisfiable (think of negating at the same time one variable attached to a vertex $`V`$ of the central plaquette and one variable in each of the two peripherical plaquettes joining the central plaquette at $`V`$), GD will never be able to find a solution and will be blocked forever in the local minimum with height $`E=1`$.

The purpose of this section is to show that this situation typically happens for random instances of XORSAT. More precisely, while almost all instances of XORSAT with a ratio of clauses per variables smaller than $`\alpha 0.918`$ have a lot of solutions, GD is almost never able to find one. Even worse, the number of violated clauses reached by GD is bounded from below by $`\mathrm{\Psi }(\alpha )N`$ where

$$\mathrm{\Psi }(\alpha )=\frac{729}{1024}\alpha ^7e^{45\alpha }.$$
(18)
In other words, the number of clauses remaining unsatisfied at the end of a typical GD run is of the order of $`N`$. Our demonstration, inspired from gdproof , is based on the fact that, with high probability, a random instance of XORSAT contains a large number of blocked islands of the type of Fig. 19.

To make the proof easier, we shall study the following fixed clause probability ensemble. Instead of imposing the number of clauses to be equal to $`M(=\alpha N)`$, any triplet $`\tau `$ of three vertices (among $`N`$) is allowed to carry a plaquette with probability $`\mu =\alpha N/\left(\genfrac{}{}{0pt}{}{N}{3}\right)=6\alpha /N^2+O(1/N^3)`$. Notice that this probability ensures that, on the average, the number of plaquettes equals $`\alpha N`$. Let us now draw a hypergraph with this distribution. For each triplet $`\tau `$ of vertices, we define $`I_\tau =1`$ if $`\tau `$ is the center of a island, 0 otherwise. We shall show that the total number of islands, $`I=_\tau I_\tau `$, is highly concentrated in the large $`N`$ limit, and calculate its average value.

The expectation value of $`I_\tau `$ is equal to

$$E[I_\tau ]=\frac{(N3)\times (N4)\times \mathrm{}\times (N13)\times (N14)}{8\times 8\times 8}\times \mu ^A(1\mu )^B,$$
(19)
where $`A=7`$ is the number of plaquettes in the island, and

$$B=\left(\genfrac{}{}{0pt}{}{N}{3}\right)\left(\genfrac{}{}{0pt}{}{N15}{3}\right)7,$$
(20)
is the number of triplets with at least one vertex among the set of 15 vertices of the island that do not carry plaquette. The significance of the terms in Eq. (19) is transparent. The central triplet $`\tau `$ occupying three vertices, we choose 2 vertices among $`N3`$ to draw the first peripherical plaquette of the island, then other 2 vertices among $`N5`$ for the other peripherical plaquette having a common vertex with the latter. The order in which these two plaquettes are built does not matter and a factor $`1/2`$ permits us to avoid double counting. The other four peripherical plaquettes have multiplicities calculable in the same way (with less and less available vertices). The terms in $`\mu `$ and $`1\mu `$ correspond to the probability that such a 7 plaquettes configuration is drawn on the 15 vertices of the island, and is disconnected from the remaining $`N15`$ vertices. The expectation value of the number $`i=I/N`$ of islands per vertex thus reads,

$$\underset{N\mathrm{}}{lim}E[i]=\underset{N\mathrm{}}{lim}\frac{1}{N}\left(\genfrac{}{}{0pt}{}{N}{3}\right)E[I_\tau ]=\frac{729}{8}\alpha ^7e^{45\alpha }.$$
(21)
Chebyshev’s inequality can be used to show that $`i`$ is concentrated around its above average value. Let us calculate the second moment of the number of islands, $`E[I^2]=_{\tau ,\sigma }E[I_\tau I_\sigma ]`$. Clearly, $`E[I_\tau I_\sigma ]`$ depends only on the number $`\mathrm{}=0,1,2,3`$ of vertices common to triplets $`\tau `$ and $`\sigma `$. It is obvious that no two triplets of vertices can be centers of islands when they have $`\mathrm{}=1`$ or $`\mathrm{}=2`$ common vertices. If $`\mathrm{}=3`$, $`\tau =\sigma `$ and $`E_{\mathrm{}=0}E[I_\tau ^2]=E[I_\tau ]`$ has been calculated above. For $`\mathrm{}=0`$, a similar calculation gives

$`E_{\mathrm{}=0}`$ $`=`$ $`{\displaystyle \frac{(N6)(N7)\mathrm{}(N29)}{2^{18}}}\mu ^{14}(1\mu )^B^{}`$ (22)
$`B^{}`$ $`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{N}{3}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{N30}{3}}\right)14.`$
Finally, we obtain

$$E[i^2]=\frac{1}{N^2}\left[\left(\genfrac{}{}{0pt}{}{N}{3}\right)E_{\mathrm{}=3}+\left(\genfrac{}{}{0pt}{}{N}{3}\right)\left(\genfrac{}{}{0pt}{}{N3}{3}\right)E_{\mathrm{}=0}\right]=E[i]^2+O\left(\frac{1}{N}\right).$$
(23)
Therefore the variance of $`i`$ vanishes and $`i`$ is, with high probability, equal to its average value given by (21). To conclude, an island has a probability $`1/2^7=1/128`$ to be blocked by definition. Therefore the number (per vertex) of blocked islands in a random XORSAT instance with ratio $`\alpha `$ is almost surely equal to $`\mathrm{\Psi }(\alpha )`$ given by Eq. (18). Since each blocked island has one unsatisfied clause, this is also a lower bound to the number of violated clauses per variable. Notice however that $`\mathrm{\Psi }(\alpha )`$ is very small and bounded from above by $`\mathrm{1.5\; 10}^9`$ over the range of interest, $`0<\alpha <0.918`$. Therefore, one would in principle need to deal with billions of variables not to reach solutions and be in the true asymptotic regime of GD.

The proof is easily generalizable to gradient descent with more than one look ahead. To extend the notion of blocked island to the case where GD is allowed to invert $`R`$, and not only 1, variables at a time, it is sufficient to have $`R+1`$, and not 2, peripherical plaquettes attached to each vertex of the central plaquette. The calculation of the lower bound $`\mathrm{\Psi }(\alpha ,R)`$ to the number of violated clauses (divided by $`N`$) reached by GD is straightforward and not reproduced here. As a consequence, GD, even with $`R`$ simultaneous flips permitting to overcome local barriers, remains almost surely trapped at an extensive (in $`N`$) level of violated clauses for any finite $`R`$. Actually the lower bound $`\mathrm{\Psi }(\alpha ,R)N`$ tends to zero only if $`R`$ is of the order of $`\mathrm{log}N`$.

We stress that the statistical physics calculation of physical ‘complexity’ $`\mathrm{\Sigma }`$ (see Sec. III.2) predicts there is no metastable states for $`\alpha <0.818`$xorsat1 , while GD is almost surely trapped by the presence of blocked islands for any $`\alpha >0`$. This apparent discrepancy comes from the fact that GD is sensible to the presence of configurations blocked for finite $`R`$, while the physical ‘complexity’ considers states metastable in the limit $`R\mathrm{}`$ onlygiu .

### III.4 The WalkSAT procedure

The Pure Random WalkSAT (PRWSAT) algorithm for solving $`K`$-SAT is defined by the following rulespapa2 .

1. Choose randomly a configuration of the Boolean variables.
2. If all clauses are satisfied, output “Satisfiable”.
3. If not, choose randomly one of the unsatisfied clauses, and one among the $`K`$ variables of this clause. Flip (invert) the chosen variable. Notice that the selected clause is now satisfied, but the flip operation may have violated other clauses which were previously satisfied.
4. Go to step 2, until a limit on the number of flips fixed beforehand has been reached. Then Output “Don’t know”.

What is the output of the algorithm? Either “Satisfiable” and a solution is exhibited, or “Don’t know” and no certainty on the status of the formula is achieved. Papadimitriou introduced this procedure for $`K=2`$, and showed that it solves with high probability any satisfiable 2-SAT instance in a number of steps (flips) of the order of $`N^2`$papa2 . Recently Schöning was able to prove the following very interesting result for 3-SATscho . Call ‘trial’ a run of PRWSAT consisting of the random choice of an initial configuration followed by $`3\times N`$ steps of the procedure. If none of $`T`$ successive trials on a given instance has been successful (has provided a solution), then the probability that this instance is satisfiable is lower than $`\mathrm{exp}(T\times (3/4)^N)`$. In other words, after $`T(4/3)^N`$ trials of PRWSAT, most of the configuration space has been ‘probed’: if there were a solution, it would have been found. Though this local search algorithm is not complete, the uncertainty on its output can be made as small as desired and it can be used to prove unsatisfiability (in a probabilistic sense).

Schöning’s bound is true for any instance. Restriction to special input distributions allows to strengthen this result. Alekhnovich and Ben-Sasson showed that instances drawn from the random 3-Satisfiability ensemble described above are solved in polynomial time with high probability when $`\alpha `$ is smaller than $`1.63`$ben .

#### III.4.1 Behaviour of the algorithm

In this section, we briefly sketch the behaviour of PRWSAT, as seen from numerical experiments Parkes  and the analysis of notrewsat ; leurwsat . A dynamical threshold $`\alpha _d`$ ($`2.7`$ for 3-SAT) is found, which separates two regimes:

* for $`\alpha <\alpha _d`$, the algorithm finds a solution very quickly, namely with a number of flips growing linearly with the number of variables $`N`$. Figure 20A shows the plot of the fraction $`\phi _0`$ of unsatisfied clauses as a function of time $`t`$ (number of flips divided by $`M`$) for one instance with ratio $`\alpha =2`$ and $`N=500`$ variables. The curve shows a fast decrease from the initial value ($`\phi _0(t=0)=1/8`$ in the large $`N`$ limit independently of $`\alpha `$) down to zero on a time scale $`t_{res}=O(1)`$. Fluctuations are smaller and smaller as $`N`$ grows. $`t_{res}`$ is an increasing function of $`\alpha `$. This relaxation regime corresponds to the one studied by Alekhnovich and Ben-Sasson, and $`\alpha _d>1.63`$ as expectedben .
* for instances in the $`\alpha _d<\alpha <\alpha _c`$ range, the initial relaxation phase taking place on $`t=O(1)`$ time scale is not sufficient to reach a solution (Fig. 20B). The fraction $`\phi _0`$ of unsat clauses then fluctuates around some plateau value for a very long time. On the plateau, the system is trapped in a metastable state. The life time of this metastable state (trapping time) is so huge that it is possible to define a (quasi) equilibrium probability distribution $`p_N(\phi _0)`$ for the fraction $`\phi _0`$ of unsat clauses. (Inset of Fig. 20B). The distribution of fractions is well peaked around some average value (height of the plateau), with left and right tails decreasing exponentially fast with $`N`$, $`p_N(\phi _0)\mathrm{exp}(N\overline{\zeta }(\phi _0))`$ with $`\overline{\zeta }0`$. Eventually a large negative fluctuation will bring the system to a solution ($`\phi _0=0`$). Assuming that these fluctuations are independant random events occuring with probability $`p_N(0)`$ on an interval of time of order $`1`$, the resolution time is a stochastic variable with exponential distribution. Its average is, to leading exponential order, the inverse of the probability of resolution on the $`O(1)`$ time scale: $`[t_{res}]\mathrm{exp}(N\zeta )`$ with $`\zeta =\overline{\zeta }(0)`$. Escape from the metastable state therefore takes place on exponentially large–in–$`N`$ time scales, as confirmed by numerical simulations for different sizes. Schöning’s resultscho  can be interpreted as a lower bound to the probability $`\overline{\zeta }(0)>\mathrm{ln}(3/4)`$, true for any instance.

The plateau energy, that is, the fraction of unsatisfied clauses reached by PRWSAT on the linear time scale is plotted on Fig. 21. Notice that the “dynamic” critical value $`\alpha _d`$ above which the plateau energy is positive (PRWSAT stops finding a solution in linear time) is strictly smaller than the “static” ratio $`\alpha _c`$, where formulas go from satisfiable with high probability to unsatisfiable with high probability. In the intermediate range $`\alpha _d<\alpha <\alpha _c`$, instances are almost surely satisfiable but PRWSAT needs an exponentially large time to prove so. Interestingly, $`\alpha _d`$ and $`\alpha _c`$ coincides for 2-SAT in agreement with Papadimitriou’s resultpapa2 . Furthermore, the dynamical transition is apparently not related to the onset of clustering taking place at $`\alpha _s3.9`$.

#### III.4.2 Results for the linear phase $`\alpha <\alpha _d`$

When PRWSAT finds easily a solution, the number of steps it requires is of the order of $`N`$, or equivalently, $`M`$. Let us call $`t_{res}(\alpha ,K)`$ the average of this number divided by the number of clauses $`M`$. By definition of the dynamic threshold, $`t_{res}`$ diverges when $`\alpha \alpha _d^{}`$. Assuming that $`t_{res}(\alpha ,K)`$ can be expressed as a series of powers of $`\alpha `$, we find the following expansionnotrewsat

$$t_{res}(\alpha ,K)=\frac{1}{2^K}+\frac{K(K+1)}{K1}\frac{1}{2^{2K+1}}\alpha +\frac{4K^6+K^5+6K^310K^2+2K}{3(K1)(2K1)(K^22)}\frac{1}{2^{3K+1}}\alpha ^2+O(\alpha ^3).$$
(24)
around $`\alpha =0`$. As only a finite number of terms in this expansion have been computed, we do not control its radius of convergence, yet as shown in Fig. 22 the numerical experiments provide convincing evidence in favour of its validity.

The above calculation is based on two facts. First, for $`\alpha <1/(K(K1))`$ the instance under consideration splits into independent subinstances (involving no common variable) that contains a number of variables of the order of $`\mathrm{log}N`$ at most. Moreover, the number of the connected components containing $`m`$ clauses, computed with probabilistic arguments very similar to those of Section III.3, contribute to a power expansion in $`\alpha `$ only at order $`\alpha ^m`$. Secondly, the number of steps the algorithm needs to solve the instance is simply equal to the sum of the numbers of steps needed for each of its independent subinstances. This additivity remains true when one averages over the initial configuration and the choices done by the algorithm.

One is then left with the enumeration of the different subinstances with a given size and the calculation of the average number of steps for their resolution. A detailed presentation of this method has been given in a general case in clusters , and applied more specifically to this problem in notrewsat ; the reader is referred to these previous works for more details. Equation (24) is the output of the enumeration of subinstances with up to three clauses.

#### III.4.3 Results for the exponential phase $`\alpha >\alpha _d`$

The above small $`\alpha `$ expansion does not allow us to investigate the $`\alpha >\alpha _d`$ regime. We turn now to an approximate method more adapted to this situation.

Let us denote by $`C`$ an assignment of the boolean variables. PRWSAT defines a Markov process on the space of the configurations $`C`$, a discrete set of cardinality $`2^N`$. It is a formidable task to follow the probabilities of all these configurations as a function of the number of steps $`T`$ of the algorithm so one can look for a simpler description of the state of the system during the evolution of the algorithm. The simplest, and crucial, quantity to follow is the number of clauses unsatisfied by the current assignment of the boolean variables, $`M_0(C)`$. Indeed, as soon as this value vanishes, the algorithm has found a solution and stops.

A crude approximation consists in assuming that, at each time step $`T`$, all configurations with a given number of unsatisfied clauses are equiprobable, whereas the Hamming distance between two configurations visited at step $`T`$ and $`T+k`$ of the algorithm is at most $`k`$. However, the results obtained are much more sensible that one could fear. Within this simplification, a Markovian evolution equation for the probability that $`M_0`$ clauses are unsatisfied after $`T`$ steps can be written. Using methods similar to the ones in Section II.2, we obtain (see notrewsat  for more details and leurwsat  for an alternative way of presenting the approximation):

* the average fraction of unsatisfied clauses, $`\phi _0(t)`$, after $`T=tM`$ steps of the algorithm. For ratios $`\alpha >\alpha _d(K)=(2^K1)/K`$, $`\phi _0`$ remains positive at large times, which means that typically a large formula will not be solved by PRWSAT, and that the fraction of unsat clauses on the plateau is $`\phi _0(t\mathrm{})`$. The predicted value for $`K=3`$, $`\alpha _d=7/3`$, is in good but not perfect agreement with the estimates from numerical simulations, around $`2.7`$. The plateau height, $`2^K(1\alpha _d(K)/\alpha )`$, is compared to numerics in Fig. 21.
* the probability $`p_N(\phi _0)\mathrm{exp}(N\overline{\zeta }(\phi _0))`$ that the fraction of unsatisfied clauses is $`\phi _0`$. It has been argued above that the distribution of resolution times in the $`\alpha >\alpha _d`$ phase is expected to be, at leading order, an exponential distribution of average $`e^{N\zeta }`$, with $`\zeta =\overline{\zeta }(0)`$. Predictions for $`\overline{\zeta }(0)`$ are plotted and compared to experimental measures of $`\zeta `$ in Fig. 23. Despite the roughness of our Markovian approximation, theoretical predictions are in qualitative agreement with numerical experiments.

A similar study of the behaviour of PRWSAT on XORSAT problems has been also performed in notrewsat ; leurwsat , with qualitatively similar conclusions: there exists a dynamic threshold $`\alpha _d`$ for the algorithm, smaller both than the satisfiability and clustering thresholds (known exactly in this case xorsat2 ). For low values of $`\alpha `$, the resolution time is linear in the size of the formula; between $`\alpha _d`$ and $`\alpha _c`$ resolution occurs on exponentially large time scales, through fluctuations around a plateau value for the number of unsatisfied clauses. In the XORSAT case, the agreement between numerical experiments and this approximate study (which predicts $`\alpha _d=1/K`$) is quantitatively better and seems to improve with growing $`K`$.

## IV Conclusion and perspectives

In this article, we have tried to give an overview of the studies that physicists have devoted to the analysis of algorithms. This presentation is certainly not exhaustive. Let us mention that use of statistical physics ideas have permitted to obtain very interesting results on related issues as number partitioningmertens , binary search trees majum , learning in neural networks nn , extremal optimization boe

It may be objected that algorithms are mathematical and well defined objects and, as so, should be analysed with rigorous techniques only. Though this point of view should ultimately prevail, the current state of available probabilistic or combinatorics techniques compared to the sophisticated nature of algorithms used in computer science make this goal unrealistic nowadays. We hope the reader is now convinced that statistical physics ideas, techniques, … may be of help to acquire a quantitative intuition or even formulate conjectures on the average performances of search algorithms. A wealth of concepts and methods familiar to physicists e.g. phase transitions and diagrams, dynamical renormalization flow, out-of-equilibrium growth phenomena, metastability, perturbative approaches… are found to be useful to understand the behaviour of algorithms. It is a simple bet that this list will get longer in next future and that more and more powerful techniques and ideas issued from modern theoretical physics will find their place in the field.

Open questions are numerous. Variants of DPLL with complex splitting heuristics, random backtrackingsmarq  or applied to combinatorial problems with internal symmetriesliat  would be worth being studied. As for local search algorithms, it would be very interesting to study refined versions of the Pure WalkSAT procedure that alternate random and greedy steps defwsat ; selmantuning ; HoosStutz  to understand the observed existence and properties of optimal strategies. One of the main open questions in this context is to what extent performances are related to intrinsic features of the combinatorial problems and not to the details of the search algorithmlancas . This raises the question of how the structure of the cost function landscape may induce some trapping or slowing down of search algorithmspar . Last of all, the input distributions of instances we have focused on here are far from being realistic. Real instances have a lot of structure which will strongly reflect on the performances of algorithms. Going towards more realistic distributions or, even better, obtaining results true for any instance would be of great interest.

Acknowledgments. This work was partly funded by the ACI Jeunes Chercheurs “Algorithmes d’optimisation et systèmes désordonnés quantiques” from the French Ministry of Research.