The backtracking survey propagation algorithm for solving random K-SAT problems

Discrete combinatorial optimization has a central role in many scientific disciplines, however, for hard problems we lack linear time algorithms that would allow us to solve very large instances. Moreover, it is still unclear what are the key features that make a discrete combinatorial optimization problem hard to solve. Here we study random K-satisfiability problems with K=3,4, which are known to be very hard close to the SAT-UNSAT threshold, where problems stop having solutions. We show that the backtracking survey propagation algorithm, in a time practically linear in the problem size, is able to find solutions very close to the threshold, in a region unreachable by any other algorithm. All solutions found have no frozen variables, thus supporting the conjecture that only unfrozen solutions can be found in linear time, and that a problem becomes impossible to solve in linear time when all solutions contain frozen variables.

Optimization problems with discrete variables are widespread among scientific disciplines and often among the hardest to be solved.The K-satisfiability (K-SAT) problem is a combinatorial discrete optimization problem of N Boolean variables, x = {x i } i=1,N , submitted to M constraints.Each constraint, called clause, is in the form of an OR logical operator of K literals (variables or their negations), and the problem is solvable when there exists at least one configuration of the variables, among the 2 N possible ones, that satisfies all constraints.The K-SAT problem for K ≥ 3 is a central problem in combinatorial optimization: it was among the first problems shown to be N P -complete [1,2,3] and is still very much studied.Efforts from the theoretical computer science community have been especially devoted to the study of the random K-SAT ensemble [4,5], where each formula is generated by randomly choosing M = αN clauses of K literals [6]; indeed formulas from this ensemble become extremely hard to solve when the clause to variable ratio α grows [6] and still the locally tree-like structure of the factor graph [7], representing the interaction network among variables, makes the random K-SAT ensemble a perfect candidate for an analytic solution.The study of random K-SAT problems and of the related solving algorithms is likely to shed light on the origin of the computational complexity and to allow for the development of improved solving algorithms.
Both numerical [8] and analytical [9,10] evidences suggest that a threshold phenomenon takes place in random K-SAT ensembles: in the limit of very large formulas, N → ∞, a typical formula has a solution for α < α s (K), while it is unsatisfiable for α > α s (K).It has been very recently proved in Ref. [11] that for K large enough the SAT-UNSAT threshold α s (K) exists in the N → ∞ limit and coincides with the prediction from the cavity method in statistical physics [12].A widely accepted conjecture is that the SAT-UNSAT threshold α s (K) exists for any value of K.The main problem is that finding solutions close to α s is very hard, and all known algorithms running in polynomial time fail to find solutions for α > α a .The algorithmic threshold α a depends on the specific algorithm, but it is well below α s for most algorithms.There are two main open questions: (i) finding improved algorithms having a larger α a , and (ii) understanding what is the theoretical upper bound to α a .In the present manuscript we present progresses on both issues.
The best prediction about the SAT-UNSAT threshold comes from the cavity method [12,13,14,15]: for example, α s (K = 3) = 4.2667 [14] and α s (K = 4) = 9.931 [15].Actually the statistical physics study of random K-SAT ensembles also provides us with a very detailed description of how the space of solutions changes when α spans the whole SAT phase (0 ≤ α ≤ α s ).Considering typical formulas in the large N limit and the vast majority of solutions in these formulas (i.e.typical solutions), we known that, at low enough α values, the solutions are many and well connected, such as to form a single cluster 1 .Increasing α not only the number of solutions decreases, but at α d the random K-SAT ensemble undergoes a phase transition and the space of solutions shatters in an exponentially large (in the problem size N ) number of clusters: two solutions belonging to different clusters have a Hamming distance O(N ).Defining the energy function E(x) as the number of unsatisfied clauses in configuration x, it has been found [12] that for α > α d there are exponentially many metastable states of positive energy, which may trap algorithms that look for solutions by energy relaxation (e.g.Monte Carlo simulated annealing).
Further increasing α, each cluster loses solutions and shrinks, but the most relevant change is in the number of clusters.The cavity method allows us to count clusters of solutions as a function of the number of solutions they contain [17] and from this very detailed description several other phase transitions have been identified [18,15].For example, at α c a condensation phase transition takes place, such that for α > α c the vast majority of solutions belong to a sub-exponential number of clusters, leading to effective long-range correlations among variables in typical solutions, which are hard to approximate by any algorithm with a finite horizon.In general α d ≤ α c ≤ α s hold.Most of the above picture of the solutions space has been proven rigorously in the large K limit [19,20].
Moving to the algorithmic side, a very interesting question is whether such a rich structure of the solutions space affects the performances of searching algorithms.While clustering at α d may have some impact on Monte Carlo based algorithms and on algorithms 1 In SAT problems we say 2 solutions are neighbors if they differ in the assignment of just one variable; in other problems (e.g.XOR-SAT [16]) this notation of neighboring needs to be relaxed, because a pair of solutions differing in just one variable are not allowed by construction.As long as the notion of neighborhood is relaxed to Hamming distances o(N ) all the picture of the solutions space based on statistical physics remains unaltered.
that try to sample solutions uniformly [21], many algorithms exist that can find at least one solution with α > α d [12,22,23].
A very likely conjecture is that the hardness of a formula is related to the existence of a subset of highly correlated variables, which is very hard to assign correctly altogether; the worst case being a subset of variables that can have a unique assignment.This concept was introduced with the name of backbone in Ref. [24].The same concept applied to solutions within a single cluster lead to the definition of frozen variables (within a cluster) as those variables taking the same value in all solutions of the cluster [25].It has been proven in Ref. [26] that the fraction of frozen variables in a cluster is either zero or lower bounded by (αe 2 ) −1/(K−2) ; in the latter case the cluster is called frozen.
According to the above conjecture, finding a solution in a frozen cluster is hard (in practice it should require a time growing exponentially with N ).So the smartest algorithm running in polynomial time should search for unfrozen clusters as long as they exist.Unfortunately counting unfrozen cluster is not an easy job, and indeed a large deviation analysis of their number has been achieved only very recently [27] for a different and simpler problem (bicoloring random regular hypergraphs).For random K-SAT only partial results are known, that can be stated in terms of two thresholds: for α > α r (rigidity) typical solutions are in frozen cluster (but a minority of solutions may still be unfrozen), while for α > α f (freezing) all solutions are frozen.It has been rigorously proven [28,29] that α f < α s holds strictly for K > 8.For small K, which is the interesting case for benchmarking solving algorithms, we only know that α f (K = 3) = 4.254 (9) from exhaustive enumerations in small formulas (N ≤ 100) [30] and α r (K = 4) = 9.883 (15) from the cavity method [15].In general The conjecture above implies that no polynomial time algorithm can solve problems with α ≥ α f , but also finding solutions close to the rigidity threshold α r is expected to be very hard, given that unfrozen solutions becomes a tiny minority.And this is indeed what happens for all known algorithms.Since we are interested in solving very large problems we only consider algorithms whose running time scales almost linearly with N and we measure performances of each algorithm in terms of its algorithmic threshold α a .
Solving algorithms for random K-SAT problems can be roughly classified in two main categories: (a) algorithms that search for a solution by performing a biased random walk in the space of configurations and (b) algorithms that try to build the solutions by assigning variables, according to some estimated marginals.WalkSat [31], focused Metropolis search (FMS) [22] and ASAT [23] belong to the former category; while in the latter category we find Belief Propagation guided Decimation (BPD) [21] and Survey Inspired Decimation (SID) [32].All these algorithms are rather effective in finding solutions to random K-SAT problems: e.g. for K = 4 we have α BPD a = 9.05, α FMS a 9.55 and α SID a 9.73 to be compared with a much lower algorithmic threshold α GUC a = 5.54 achieved by Generalized Unit Clause, the best algorithm whose convergence to a solution can be proven rigorously [33].Among the efficient algorithms above, only BPD can be solved analytically [21] to find the algorithmic threshold α BPD a ; for all the others we are forced to run extensive numerical simulations in order to measure α a .
At present the algorithm achieving the best performances on several constraint satisfaction problems is SID, which has been successfully applied to the random K-SAT problem [12] and to the coloring problem [34].The statistical properties of the SID algorithm for K = 3 have been studied in details in Refs.[35,32].Numerical experiments on random 3-SAT problems with a large number of variables, up to N = 3 × 10 5 , show that in a time that is approximately linear in N the SID algorithm finds solutions up to α SID a 4.2525 [35], that is definitely smaller, although very close to, α s (K = 3) = 4.2667.In the region α SID a < α < α s the problem is satisfiable for large N , but at present no algorithm can find solutions there.
To fill this gap we study a new algorithm for finding solutions to random K-SAT problems, the Backtracking Survey Propagation (BSP) algorithm.This algorithm (fully explained in the Methods section) is based, as SID, on the survey propagation (SP) equations derived within the cavity method [12,35,32] that provide an estimate on the total number of clusters N clus = exp(Σ).The BSP algorithm, like SID, aims at assigning gradually the variables such as to keep the complexity Σ as large as possibile, i.e. trying not to kill too many clusters [35].While in SID each variable is assigned only once, in BSP we allow unsetting variables already assigned such as to backtrack on previous non-optimal choices.In BSP the r parameter is the ratio between the number of backtracking moves (unsetting one variable) and the number of decimation moves (assigning one variable).r < 1 must hold and for r = 0 we recover the SID algorithm.Running times scale as 1/(1 − r), but the algorithm complexity remains practically linear in N for any r < 1.
The idea supporting backtracking [36] is that a choice made at the beginning of the decimation process, when most of the variables are unassigned, may turn to be suboptimal later on, and re-assigning a variable that is no longer consistent with the current best estimate of its marginal probability may lead to  and N = 50 000 (empty symbols).On instances not solved, a second run rarely (< 1%) finds a solution.a better satisfying configuration.We do not expect the backtracking to be essential when correlations between variables are short ranged, but approaching α s we know that correlations become long ranged and thus the assignment of a single variable may affect a huge number of other variables: this is the situation when we expect the backtracking to be essential.
This idea may look similar in spirit to the survey propagation reinforcement (SPR) algorithm [37], where variables are allowed to change their most likely value during the run, but in practice BSP works much better.In SPR, once reinforcement fields are large, the re-assignment of any variable becomes unfeasible, while in BSP variables can be re-assigned to better values until the very end, and this is a major advantage.

Results
In this section we show the results of extensive numerical simulations using BSP.In each plot having α on the abscissa, the right end of the plot coincides with the best estimate of α s , in order to provide an immediate indication of how close to the SAT-UNSAT threshold the algorithm can work.
Probability of finding a SAT assignment The standard way to study the performances of a solving algorithm is to measure the fraction of instances it can solve as a function of α.We show in Fig. 1 such a fraction for BSP run with three values of the r parameter (r = 0, 0.5 and 0.9) on random 4-SAT problems of two different sizes (N = 5 000 and N = 50 000).The probability of finding a solution increases both with r and N , but an extrapolation to the large N limit of these  The upper panel shows Σ res /N res as a function of α, for random 4-SAT problems solved by BSP with r = 0.9.This order parameter is mostly size-independent (few deviations appear only at the right end) and vanishes at α BSP a ≈ 9.9, slightly beyond the rigidity threshold α r = 9.883 (15), marked by a vertical line.The lower panel shows that the same linear extrapolation holds for other values of r, the black line being the fit to r = 0.9 data shown in the upper panel, but SID without backtracking (r = 0) has a much lower algorithmic threshold, α SID a ≈ 9.83.
data is unlikely to provide a reliable estimation of the algorithmic threshold α a .
Order parameter and algorithmic threshold In order to obtain a reliable estimate of α a we look for an order parameter vanishing at α a and having very little finite size effects.We identify this order parameter with the quantity Σ res /N res , where Σ res and N res are respectively the complexity (i.e.log of number of clusters) and the number of unassigned variables in the residual formula.As explained in Methods, BSP assigns and re-assigns variables, thus modifying the formula, until the formula simplifies enough that the SP   9) measured on small problems [30], and marked by the vertical line.A linear fit to the residual complexity (red line) extrapolates to zero slightly beyond the SAT-UNSAT threshold, at α BSP a ≈ 4.268, strongly suggesting BSP can find solutions in the entire SAT phase for K = 3 in the large N limit.The black line is a linear fit vanishing at α s .fixed point has only null messages: the residual formula is defined as the last formula with non-null SP fixed point messages.We have experimentally observed that the BSP algorithm (as the SID one [35]) can simplify the formula enough to reach the trivial SP fixed point only if the complexity Σ remains strictly positive during the whole decimation process.In other words, on every run where Σ becomes very close to zero or negative, SP stops converging or a contradiction is found.This may happen either because the original problem was unsatisfiable or because the algorithm made some wrong assignments incompatible with the few available solutions.Thanks to the above observation we have that Σ res ≥ 0 and thus a null value for the mean residual complexity signals that the BSP algorithm is not able to find any solution, and thus provides a valid estimate for the algorithmic threshold α BSP a .As we see in the data shown in the upper panel of Fig. 2 the mean value of the intensive mean residual complexity Σ res /N res is practically size-independent and a linear fit provides a very good data interpolation: tiny finite size effects are visible in the largest N datasets only close to the dataset right end.The linear extrapolation predicts α BSP a ≈ 9.9 (for K = 4 and r = 0.9), which is slightly above the rigidity threshold α r = 9.883 (15) computed in Ref. [15].This means that BSP is able to find solutions in a region of α where the majority of solutions is in frozen clusters and thus hard  to find.We show below that BSP actually finds solutions in atypical unfrozen clusters, as it has been observed already in very smart algorithms solving other kind of constraint satisfaction problems [38,39].
The effectiveness of the backtracking can be appreciated in the lower panel of Fig. 2, where the order parameter Σ res /N res is shown for r = 0 and r = 0.5, together with linear fits to these datasets and to the r = 0.9 dataset (black line).We observe that the algorithmic threshold for BSP is much larger (on the scale measuring the relative distance from the SAT-UNSAT threshold) that the one for SID (i.e.r = 0 dataset).
For random 3-SAT the algorithmic threshold of BSP, run with r = 0.9, practically coincide with the SAT-UNSAT threshold α s (see Fig. 3), thus providing a strong evidence that BSP can find solutions in the entire SAT phase.The estimate for the freezing threshold α f = 4.254(9) obtained in Ref. [30] with N ≤ 100 is likely to be too small, given that all solutions found by BSP for N = 10 6 are unfrozen, even for α > α f .
Computational complexity As explained in Methods, the BSP algorithm performs O(f −1 (1−r) −1 ) elementary steps, where at each step f N variables are either assigned [with prob.1/(1 + r)] or released [with prob.r/(1 + r)].At the beginning of each step, the algorithm solves the SP equations with a mean number η of iterations2 , and then sort local marginals.The total number of elementary operations to solve an instance grows as f −1 (1−r) −1 (ηN +N log N ), i.e. the algorithm running time is at most O(N log N ).Fig. 4 shows that η is actually a small number changing mildly with α and N both for K = 3 and K = 4.The main change that we observe is in the fluctuations of η that for K = 3 become larger beyond α f (marked by a vertical line).We expect η to eventually grow as O(log(N )), but for the sizes studied we do not observe such a growth.Moreover, given that the sorting of local marginals does not need to be strict (i.e. a partial sorting [40] running in O(N ) time can be enough), we have that in practice the algorithm runs in a time almost linear in the problem size N .
Whitening procedure Given that the BSP algorithm is able to find solutions even slightly beyond the rigidity threshold α r , it is natural to check whether these solutions have frozen variables or not.We concentrate on solutions found for random 3-SAT problems with N = 10 6 , since the large size of these problems makes the analysis very clean, and because we have a good number of solutions beyond the best estimate for the freezing threshold α f = 4.254(9) [30].
On each solution found we run the whitening procedure (first introduced in [41] and deeply discussed in [42,26]), that identifies frozen variables by assigning the joker state to unfrozen variables, i.e. variables that can take more than one value without violating any clause and thus keeping the formula satisfied.At each step of the whitening procedure, a variables is considered unfrozen (and thus assigned to ) if it belongs only to clauses which either involve a variable or are satisfied by another variable.The procedure is continued until all variables are or a fixed point is reached: non-variables at the fixed point correspond to frozen variables in the starting solution.
We uncover that all solutions found by BSP are converted to all-by running the whitening procedure, thus showing that solutions found by BSP have no frozen variables.This is somehow expected, according to the conjecture discussed in the Introduction: finding solutions in a frozen cluster would take an exponential time, and so the BSP algorithm actually finds solutions at very large α values by smartly focusing on the sub-dominant unfrozen clusters.
The whitening procedure leads to a relaxation of the number of non-variables as a function of the number of iterations t that follows a two steps relaxation process [25] with an evident plateau (see upper panel in Fig. 5), that becomes longer increasing α towards the algorithmic threshold.The time for leaving the plateau, scales as the time τ (c) for reaching a fraction c on non-variables (with c smaller than the plateau value).The latter has large fluctuations from solution to solution, as shown in the central panel of Fig. 5 for c = 0.4 (very similar, but shifted, histograms are obtained for other c values).However, after leaving the plateau, the dynamics of the whitening procedure is the same for each solution.Indeed plotting the mean  The upper panel shows the mean fraction of non-variables during the whitening procedure applied to all solutions found by the BSP algorithm.The line is the SP prediction for the fraction of frozen variables in typical solutions at α = 4.262, and shows that solutions found by BSP are atypical.The central panel shows the histograms of the whitening times, defined as the number of iterations required to reach a fraction 0.4 of nonvariables.Increasing α both the mean value and the variance of the whitening times grow.The lower panel shows the same data used for the upper one, but averaged at fixed τ (0)−t, the time to reach the all-fixed point.Errors are tiny since the whitening dynamics below the plateau is the same for each solution.fraction of non-variables as a function of the time to reach the all-configuration, τ (0) − t, we see that fluctuations are strongly suppressed and the relaxation is the same for each solution (see lower panel in Fig. 5).

Critical exponent for the whitening time divergence
In order to quantify the increase of the whitening time approaching the algorithmic threshold, and inspired by critical phenomena, we check for a power law divergence as a function of (α a − α) or Σ res , which are linearly related.In Fig. 6 we plot in a double logarithmic scale the mean whitening time τ (c) as a function of the residual complexity Σ res , for different choices of the fraction c of non-variables defining the whitening time.Data points are interpolated via the power law τ (c) = A(c) + B(c)Σ −ν res , where the critical exponent ν is the same for all the c values.Joint interpolations return the following best estimates for the critical exponent: ν = 0.269 (5) for K = 4 and ν = 0.281 (6) for K = 3.The two estimate turn out to be well compatible within errors, thus suggesting a sort of universality for the critical behavior approaching the algorithmic threshold α BSP a .Nonetheless a word of caution is needed since the solutions we are using as starting points for the whitening procedure are atypical solutions (otherwise they would likely contain frozen variables and would not flow to the all-configuration under the whitening procedure).So, while finding universal critical properties in a dynamical process is definitely a good news, how to relate it to the behavior of the same process on typical solutions it is not obvious (and indeed for the whitening process starting from typical solutions one would expect the naive mean field exponent ν = 1/2, which is much larger than the one we are finding).

Discussion
We have studied the Backtracking Survey Propagation (BSP) algorithm for finding solutions in very large random K-SAT problems and provided numerical evidence that it works much better than any previously available algorithm.That is, BSP has the largest algorithmic threshold known at present.The main reason for its superiority is the fact that variables can be re-assigned at any time during the run, even at the very end.In other solving algorithms that may look similar, as e.g.survey propagation reinforcement [37], re-assignment of variables actually takes place mostly at the beginning of the run, and this is far less efficient in hard problems.Even doing a lot of helpful backtracking, the BSP running time is still O(N log N ) in the worst case, and thanks to this it can be used on very large problems with millions of clauses.
For K = 3 the BSP algorithm finds solutions practically up to the SAT-UNSAT threshold α s , while for K = 4 a tiny gap to the SAT-UNSAT threshold still remains, but the algorithmic threshold α BSP a seems to be located beyond the rigidity threshold α r in the large N limit.Beating the rigidity threshold, i.e. finding solutions in a region where the majority of solutions belongs to clusters with frozen variables, is hard, but not impossible.Even under the assumption that finding frozen solutions takes an exponential time in N , very smart polynomial time algorithms can look for a solution in the sub-dominant unfrozen clusters [38,39].BSP belongs to this category, as we have shown that all solutions found by BSP have no frozen variables.
One of the main questions we tried to answer with our extensive numerical simulations is whether BSP is reaching (or approaching closely) the ultimate threshold α a for polynomial time algorithms solving large random K-SAT problems.Under the assumption that frozen solutions can not be found in polynomial time, such an algorithmic threshold α a would coincide with the freezing transition at α f (i.e. when the last unfrozen solution disappears).Unfortunately for random K-SAT the location of α f is not known with enough precision to allow us to reach a definite answer to this question.It would be very interesting to run BSP on random hypergraph bicoloring problems, where the threshold values are known [43] and a very recent work has shown that the large deviation function for the number of unfrozen clusters can be computed [27].
It is worth noticing that the BSP algorithm is easy to parallelize, since most of the operations are local and do not require any strong centralized control.Obviously the effectiveness of a parallel version of the algorithm would largely depend on the topology of the factor graph representing the specific problem: if the factor graph is an expander, then splitting the problem on several cores may require too much inter-core bandwidth, but in problems having a natural hierarchical structure the parallelization may lead to further performance improvements.
The backtracking introduced in the BSP algorithm helps a lot in correcting errors made during the partial assignment of variables and this allow the BSP algorithm to reach solutions at large α values.Clearly the price we pay is that a too frequent backtracking makes the algorithm slower.A natural direction to improve this class of algorithms would be to used biased marginals focusing on solutions which are easier to be reached by the algorithm.For example in the region α > α r the measure is concentrated on solutions with frozen variables, but these can not be really reached by the algorithm.The backtracking thus intervenes and corrects the partial assignment until a solution with unfrozen variables is found by chance.If the marginals could be computed from a new biased measure which is concentrated on the unfrozen clusters, this could make the algorithm go immediately in the right direction and much less backtracking would be hopefully needed.

Methods
The Backtracking Survey Propagation (BSP) algorithm is an extension of the Survey Inspired Decimation (SID) algorithm for finding solutions to random K-satisfiability problems.A detailed description of the SID algorithm can be found in Refs.[12,13,32].Here, we briefly recall SID and describe the novelties leading to the BSP algorithm [36].
Survey Inspired Decimation (SID) The SID algorithm is based on the survey propagation (SP) equations derived by the cavity method [12,13], that can be written in a compact way as ma→i = j∈∂a\i m j→a , where ∂ a is the set of variables in clause a, and ∂ + ia (resp.∂ − ia ) is the set of clauses containing x i , excluding a itself, satisfied (resp.not satisfied) when the variable x i is assigned to satisfy clause a.
Physical interpretation of the SP equations: ma→i represents the fraction of clusters where clause a is satisfied solely by variable x i (that is, x i is frozen by clause a), while m i→a is the fraction of clusters where x i is frozen to an assignment not satisfying clause a.
The SP equations impose 2KM self-consistency conditions on the 2KM variables {m i→a , ma→i } living on the edges of the factor graph [7], that are solved in an iterative way, leading to a message passing algorithm (MPA) [4], where outgoing messages from a factor graph node (variable or clause) are functions of the incoming messages.Once the MPA reaches a fixed point {m i→a , m a→i } that solves the SP equations, the number of clusters can be estimated via the complexity where K a is the length of clause a (initially K a = K) and ∂ + i (resp.∂ − i ) is the set of clauses satisfied by setting x i = 1 (resp.x i = −1).The SP fixed point messages also provide information about the fraction of clusters where variable x i is forced to be positive (w + i ), negative (w − i ) or not forced at all (1 − w + i − w − i ) The SID algorithm then proceed by assigning variables (decimation step).According to SP equations, assigning a variable x i to its most probable value (i.e., setting x i = 1 if w + i > w − i and viceversa), the number of clusters gets multiplied by a factor, called bias With the aim of decreasing the lesser the number of cluster and thus keeping the largest the number of solutions in each decimation step, SID assigns/decimate variables with the largest b i values.In order to keep the algorithm efficient, at each step of decimation a small fraction f of variables is assigned, such that in O(log N ) steps of decimation a solution can be found.After each step of decimation, the SP equations are solved again on the subproblem, which is obtained by removing satisfied clauses and by reducing clauses containing a false literal (unless a zero-length clause is generated, and in that case the algorithm returns a failure).The complexity and the biases are updated according to the new fixed point messages, and a new decimation step is performed.
The main idea of the SID algorithm is that fixing variables which are almost certain to their most probable value, one can reduce the size of the problem without reducing too much the number of solutions.The evolution of the complexity Σ during the SID algorithm can be very informative [35].Indeed it is found that, if Σ becomes too small or negative, the SID algorithm is likely to fail, either because the iterative method for solving the SP equations no longer converges to a fixed point or because a contradiction is generated by assigning variables.In these cases the SID algorithm returns a failure.On the contrary, if Σ always remains well positive, the SID algorithm reduces so much the problem, that eventually a trivial SP fixed point, m i→a = ma→i = 0, is reached.This is a strong hint that the remaining subproblem is easy and the SID algorithm tries to solve it by WalkSat [31].
A careful analysis of the SID algorithm for random 3-SAT problems of size N = O(10 5 ) shows that the algorithmic threshold achievable by SID is α SID a = 4.2525 [35], which is close, but definitely smaller than the SAT-UNSAT threshold α s = 4.2667.
The running time of the SID algorithm experimentally measured is O(N log(N )).However it is not excluded than an extra factor log(N ) could be present in the convergence time of the iterative solution to SP equations, but for sizes up to N = O(10 7 ) has not been observed [32].
Backtracking Survey Propagation (BSP) Willing to improve the SID algorithm to find solutions also in the region α SID a < α < α s , one has to reconsider the way variables are assigned.It is clear that the fact the SID algorithm assigns each variable only once is a strong limitation, especially in a situation where correlations between variables becomes extremely strong and long-ranged.In difficult problems it can easily happen that one realizes that a variable is taking the wrong value only after having assigned some of its neighbours variables.However the SID algorithm is not able to solve this kind of frustrating situations.
The BSP algorithms tries to solve this kind of problematic situations by introducing a new backtracking step, where a variable already assigned can be released and eventually re-assigned in a future decimation step.It is not difficult to understand when it is worth releasing a variable.The bias b i in terms of the SP fixed point messages { m a→i } a∈∂ i arriving in i can be computed also for a variable x i already assigned: if the bias b i , that was large at the time the variable x i was assigned, gets strongly reduces by the effect of assign-ing other variables, then it is likely that releasing the variable x i may be beneficial in the search for a solution.So both variables to be fixed in the decimation step and variables to be released in the backtracking step are chosen according to their biases b i : variables to be fixed have the largest biases and variables to be released have the smallest biases.
The BSP algorithm then proceeds similarly to SID, by alternating the iterative solution to the SP equations and a step of decimation or backtracking on a fraction f of variables in order to keep the algorithm efficient (in all our numerical experiments we have used f = 10 −3 ).The choice between a decimation or a backtracking step is taken according to a stochastic rule (unless there are no variables to unset), where the parameter r ∈ [0, 1) represents the ratio between backtracking steps to decimation steps.Obviously for r = 0 we recover the SID algorithm, since no backtracking step is ever done.Increasing r the algorithm becomes slower by a factor 1/(1−r), because most variables are reassigned O(1/(1−r)) times before the BSP algorithm finishes, but its time complexity remains O(N log N ) in the problem size.
The BSP algorithm can stop for the same reasons the SID algorithm does: either the SP equations can not be solved iteratively or the generated subproblem has a contradiction.Both cases happen when the complexity Σ becomes too small or negative.On the contrary if the complexity remain always positive the BSP eventually generate a subproblem where all SP messages are null and on this subproblem WalkSat is called.
The code used for collecting data shown in this study is available from the authors upon request.

Figure 1 :
Figure1: Fraction of random 4-SAT instances solved by BSP.Such a fraction is computed over 100 instances at each α value for N = 5 000 (solid symbols) and N = 50 000 (empty symbols).On instances not solved, a second run rarely (< 1%) finds a solution.

Figure 2 :
Figure2: BSP algorithmic threshold on random 4-SAT problems.The upper panel shows Σ res /N res as a function of α, for random 4-SAT problems solved by BSP with r = 0.9.This order parameter is mostly size-independent (few deviations appear only at the right end) and vanishes at α BSP

Figure 3 :
Figure 3: BSP algorithmic threshold on random 3-SAT problems.Same as Fig. 2 for K = 3. BSP does not show any change at the freezing threshold α f = 4.254(9) measured on small problems[30], and marked by the vertical line.A linear fit to the residual complexity (red line) extrapolates to zero slightly beyond the SAT-UNSAT threshold, at α BSP

Figure 4 :
Figure 4: BSP convergence time.The mean number η of iterations to reach a fixed point of SP equations grows very mildly with α and N , both for K = 3 (upper panel) and K = 4 (lower panel).

Figure 5 :
Figure 5: Whitening random 3-SAT solutions.The upper panel shows the mean fraction of non-variables during the whitening procedure applied to all solutions found by the BSP algorithm.The line is the SP prediction for the fraction of frozen variables in typical solutions at α = 4.262, and shows that solutions found by BSP are atypical.The central panel shows the histograms of the whitening times, defined as the number of iterations required to reach a fraction 0.4 of nonvariables.Increasing α both the mean value and the variance of the whitening times grow.The lower panel shows the same data used for the upper one, but averaged at fixed τ (0)−t, the time to reach the all-fixed point.Errors are tiny since the whitening dynamics below the plateau is the same for each solution.

Figure 6 :
Figure6: Critical exponent for the whitening time divergence.The whitening time τ (c), defined as the mean time needed to reach a fraction c of nonvariables in the whitening procedure, is plotted in a double logarithmic scale as a function of Σ res for random 3-SAT problems with N = 10 6 (upper dataset) and random 4-SAT problems with N = 5 × 10 4 (lower dataset).The lines are power law fits with exponent ν = 0.281(6) for K = 3 and ν = 0.269(5) for K = 4.