Boolean satisfiability [1] (k-SAT) is one of the most studied optimization problems, as an efficient (that is, polynomial-time) solution to k-SAT (for k ≥ 3) implies efficient solutions to a large number of hard optimization problems [2, 3] . Here we propose a mapping of k-SAT into a deterministic continuous-time dynamical system with a unique correspondence between its attractors and the k-SAT solution clusters. We show that beyond a constraint density threshold, the analog trajectories become transiently chaotic [4, 5, 6, 7] , and the boundaries between the basins of attraction [8] of the solution clusters become fractal [7, 8, 9] , signaling the appearance of optimization hardness [10] . Analytical arguments and simulations indicate that the system always finds solutions for satisfiable formulae even in the frozen regimes of random 3- SAT [11] and of locked occupation problems [12] (considered among the hardest algorithmic benchmarks); a property partly due to the systems hyperbolic [4, 13] character. The system finds solutions in polynomial continuous-time, however, at the expense of exponential fluctuations in its energy function. * The article appeared in Nature Physics 7, 966 (2011) † E-mail: ercsey.ravasz@phys.ubbcluj.ro ‡ E-mail: toro@nd.edu Boolean satisfiability [1] (k-SAT, k ≥ 3) is the quintessential constraint satisfaction problem, lying at the basis of many decision, scheduling, error-correction and bio-computational applications. k-SAT is in NP, that is its solutions are efficiently (polynomial time) checkable, but no efficient (polynomial time) algorithms are known to compute those solutions [2] . If such algorithms would be found for k-SAT, all NP problems would be efficiently computable, since k-SAT is NP-complete [2, 3] . In k-SAT there are given N Boolean variables {x 1 , . . . , x N }, x i ∈ {0, 1} and M clauses (constraints), each clause being the disjunction (OR, denoted as ∨) of k variables or their negation (x). One has to find an assignment of the variables such that all clauses (called collectively as a formula) are satisfied (TRUE = 1). When the number of constraints is small, it is easy to find solutions, while for too many constraints it is easy to decide that the formula is unsatisfiable (UNSAT). Deciding satisfiability, in the 'intermediate range', however, can be very hard: the worst-case complexity of all known algorithms for k-SAT is exponential in N . Inspired by the mechanisms of information processing in biological systems, analog computing received increasing interest from both theoretical [14, 15, 16] and engineering communities [17, 18, 19, 20, 21] . Although the theoretical possibility of efficient computation via chaotic dynamical systems has been shown previously [15] , nonlinear dynamical systems theory has not been exploited for NP-complete problems in spite of the fact that, as shown by Gu et al.[19], Nagamatu et al. [20] and Wah et al. [21] , k-SAT can be formulated as a continuous global optimization problem [19] , and even cast as an analog dynamical system [20, 21] . Here we present a novel continuous-time dynamical system for k-SAT, with a dynamics that is rather different from previous approaches. Let us introduce the continuous variables [19] s i ∈ [−1, 1] , such that s i = −1 if the i-th variable (x i ) is FALSE and s i = 1 if it is TRUE. We define c mi = 1 for the direct form (x i ), c mi = −1 for the negated form (x i ), and c mi = 0 for the absence of the i-th variable from clause m. Defining the constraint function (1 − c mi s i ) corresponding to clause m, we have K m ∈ [0, 1] and K m = 0 if and only if clause m is satisfied. The goal would be to find a solution s * with s * i ∈ {−1, 1} to E(s * ) = 0, where E is the energy function E(s) = M m=1 K m (s) 2 . If such s * exists, it will be a global minimum for E and a solution to the k-SAT problem. However, finding s * by a direct minimization of E(s) will typically fail due to non-solution attractors trapping the search dynamics. In order to avoid such traps, here we define a modified energy function V (s, a) = M m=1 a m K m (s) 2 , using auxiliary variables a m ∈ (0, ∞) similar to Lagrange multipliers [20, 21]  for any a ∈ (0, ∞) M . We now introduce a continuous-time dynamical system on Ω through: where ∇ s is the gradient operator with respect to s, and The initial conditions for s are arbitrary s(0) ∈ H N , however, for a they have to be strictly positive, a m (0) > 0 (e.g., a m (0) = 1). The k-SAT solutions s * ∈ V N are fixed points of (1-2), for any a ∈ (0, ∞) M . The k-SAT solution clusters are spanning piecewise compact, connected sets in Q N , and every point in them is a fixed point of (1-2) (Supplementary sect. A). System (1-2) has a number of key properties (see Supplementary Information). (i) The dynamics in s stays confined to H N . (ii) The k-SAT solutions s * ∈ V N are attractive fixed points of (1-2). In particular, every point s from the orthant of a k-SAT solution s * with the property |s| 2 ≥ N − 1 + (k − 1) 2 /(k + 1) 2 is guaranteed to flow into the attractor corresponding to s * . (iii) There are no limit cycles. (iv) For satisfiable formulae the only fixed point attractors of the dynamics are the global minima of V with V = 0. Note that in principle, the projection of the dynamics onto H N could be stuck in some point s, while da/dt = 0 indefinitely. This does not happen here, as shown in Supplementary sect. E. Moreover, analytical arguments supported by simulations indicate that the trajectory will leave any domain that does not contain solutions, see the discussion in Supplementary sect. E1. Note, that the constraint functions (hence their satisfiability) depend directly only on the location of the trajectory in , and not on the auxiliary variables. The dynamics in the a-space is simple expansion, and for this reason the features of the full phase space Ω lie within its projection onto H N . One can actually eliminate entirely the auxiliary variables from the equations by first solving (2) to give a m (t) = a m (0) exp . Another fundamental feature of (1-2) is that it is deterministic: for a given formula f , any initial condition generates a unique trajectory, and any set from H N has a unique preimage arbitrarily back in time. Hence, the characteristics of the solution space are reflected in the properties of the invariant sets [7] of the dynamics (1-2) within the hypercube H N . The deterministic nature of (1-2) allows us to define basins of attractions of solution clusters by colouring every point in H N according to which cluster the trajectory flows to, if started from there. These basins fill H N up to a set of zero (Lebesgue) measure, which forms the basin boundary [7] , from where the dynamics (by definition) cannot flow to any of the attractors. A k-SAT formula f can be represented as a hypergraph G(f ) (or equivalently, a factor graph) in which nodes are variables and hyperedges are clauses connecting the nodes/variables in the clause. Pure literals are those that participate in one or more clauses but always in the same form (direct or negated); hence they can always be chosen such as to satisfy those clauses. The core of G(f ) is the subgraph left after sequentially removing all the hyperedges having pure literals [23] . For simple formulae (such as those without a core), the dynamics of (1-2) is laminar flow and the basin boundaries form smooth, non-fractal sets (Fig.1a,c and Fig.2 top two rows) . Adding more constraints G(f ) develops a core, the spin equations (1) become mutually coupled, and the trajectories may become chaotic (Fig.1b,  Supplementary sect. F, Fig.12 ) and the basin boundaries fractal [7, 8, 9] (Fig. 1d, Fig. 2, Fig. 8 ). Therefore, as the constraint density α = M/N is increased within predefined ensembles of formulae (random k-SAT, occupation problems, k-XORSAT, etc.) a sharp change to chaotic behaviour is expected at a chaotic transition point α χ , where a chaotic core appears with non-zero statistical weight in the ensemble as N → ∞. As an example, let us consider 3-XORSAT. In this case, due to its inherently linear nature, it is actually better to work directly with the parity check equations as constraints, instead of their CNF form. The chaotic core here is a small finite hypergraph, and The colour indicates the energy E (colour bar) in a given point of the trajectory. While for easy formulae the trajectories exhibit laminar flow, for hard formulae they quickly become separated, showing a chaotic evolution. Taking a small 3-XORSAT instance with N = 15 (see Supplementary sect. G) we fix a random initial condition for all s i , except s 1 and s 2 which are varied on a 400 × 400 grid and we colour each point according to the solution they flow to for c) γ = 0.6 (instance shown on Fig. 8e ) and d) for γ = 0.8 (instance shown in Fig. 8f) ). Fig. 8 ). Note, a core can be non-chaotic, and thus the existence of a core is only a necessary condition for the appearance of chaos and in general the two transitions might not coincide. Further increasing the number of constraints (within any formula ensemble) unsatisfiability appears at the threshold value α s > α χ beyond which almost all formulae are unsatisfiable (UNSAT regime) [11, 12, 22, 24, 25, 26, 27, 28] . The closer α is to α s , the harder it is to find solutions, and beyond the so-called freezing transition point α f < α s (called the frozen regime) all known algorithms take exponentially long times or simply fail to find solutions [11, 12] . A variable is frozen if it takes on the same value for all solutions within a cluster, and a cluster is frozen if an extensive number of its variables are frozen. In the frozen regime all clusters are frozen and they are also far apart ( O(N ) Hamming distance) [11, 12] In LOPs all clusters are formed by exactly one solution, hence they are completely frozen and the frozen regime extends from the clustering (dynamical) transition point l d to the satisfiability threshold l s , and thus it is very wide [12] . An example LOP is random +1-in-3-SAT [12], made of constraints that have no negated variables and a constraint is satisfied only if exactly one of its variables is 1 (TRUE). In +1-in-3-SAT l d ∼ = 2.256, l s ∼ = 2.368, and beyond l d all known algorithms have exponential search times or fail to find solutions (here l = 3M/N ). As chaos is present for satisfiable formulae, that is, when system (1-2) has attracting fixed points, it is necessarily of transient type. Transient chaos [4, 5, 6, 7] is ubiquitous in systems with many degrees of freedom such as fluid turbulence [30] . It appears as the result of homoclinic/heteroclinic intersections of the invariant manifolds of hyperbolic (unstable) fixed points of (1-2) lying within the basin boundary [7, 8, 9] , leading to complex (fractal) foliations of the phase space (see Supplementary sect. F). We observed the prevalence of transient chaos in the whole region α χ < α < α s for all the problem classes we studied. Interestingly, the velocity fluctuations of trajectories in the chaotic regime are qualitatively similar to those of fluid parcels in turbulent flows as shown in Supplementary sect. K. Our findings suggest that chaotic behaviour may be a generic feature of algorithms searching for solutions in hard optimization problems, corroborating the observations by Elser et al. [10] using a heuristic algorithm based on iterated maps. In the following we show results on random 3-SAT and +1-in-3-SAT formulae in the frozen regime, however, the same conclusions hold for other ensembles that we tested. To investigate the complexity of computation by the flow (1-2), we monitored the fraction of problems p(t) not solved by continuous time t, as function of N and α. Figs. 3a,c show that even in the frozen phase, the fraction of unsolved problems by time t decays exponentially with t, that is, by a law p(t) = re −λ(N )t . The decay rate λ(N ) obeys λ(N ) = bN −β , with β 1.6 in both cases, see Fig. 3b,d. From these two equations, the continuous-time t(p, N ) needed to solve a fixed (1 − p)th fraction of random formulae (or to miss solving p-th fraction of them) is: Attractor basins. For a random 3-SAT instance with N = 50 we vary α by successively adding new constraints. Fixing a random initial condition for s i , i ≥ 3, we vary only s 1 and s 2 on a 400 × 400 grid, and we colour each point according to the solution (first column) or solution cluster (second column) they flow to. Each colour in a given column represents a solution or solution cluster respectively, however colours between columns are independent. The third column represents the analog search time t needed to find a solution (see colour bar) starting from the corresponding grid point. Maps are presented for values of α = 3.5, 3.7, 3.9, 4.1, 4.16, 4.2, 4.24. Easy formulae are characterized by smooth basin boundaries and small search times. Note that we only see the solutions (and clusters) that reveal themselves in the (s 1 , s 2 ) plane, others might not be seen. For hard formulae the boundaries and the search time maps become fractal. = 20, 30, 40, 50, 60, 80, 100, 125, 150 (colours). Averages were done over 10 5 instances for each N . For each instance the dynamics was started from one random initial condition. Black continuous lines show the decay p(t) = r exp(−λ(N )t). b) The decay rate follows λ(N ) = bN −β , with β 1.66. c) The fraction of problems p(t) unsolved by time t for +1-in-3-SAT at l = 2.34, for N = 20, 25, 30, 35, 40, 50, 60, 70, 80. For each instance the dynamics was started in parallel from 10 random initial conditions, averages were taken over 10 4 instances for each N . Black continuous lines show the same exponential decay as in a). d) The decay rate shows the same behaviour as in b) with exponent: β 1.68. indicating that the continuous time (CT) needed to find solutions scales as a power-law with N . Eq. (3) also implies power-law scaling for almost all hard instances in the N → ∞ limit (Supplementary sect. H). The length in H N of the corresponding continuous trajectories also scales as a power-law with N ( Supplementary Fig. 11b , sect. J). However, note that this does not mean that the algorithm itself is a polynomial-cost algorithm, as the energy function V can have exponentially large fluctuations. As the numerical integration happens on a digital machine, it approximates the continuous trajectory with discrete points. Monitoring the fraction of formulae left unsolved as function of the number of discretization steps n step in the frozen phase, we find exponential behaviour for n step (p, N ) (Supplementary sects. I,J, Fig. 10) . The difference between the continuous-and discrete-time complexities is due to the wildly fluctuating nature of the chaotic trajectories (see Fig. 1b and Methods) in the frozen phase. Compounding this, we also observe the appearance of the Wada property [7, 8] in the basin boundaries, Fig. 4 . A fractal basin boundary has Wada property if its points are simultaneously on the boundary of at least three colours/basins. (An amusing method that creates such sets uses four Christmas ball ornaments [7] .) Although the Wada property does not affect the true/mathematical analog trajectories, owing to numerical errors, it may switch the numerical trajectories between the basins. Since the clusters are far (O(N )) apart, the switched trajectory will flow towards another cluster into a practically opposing region of H N until it may come close again to the basin boundary etc., partially randomizing the trajectory in H N . Figure 4 : Wada property. Basin boundaries are shown for +1-in-3-SAT for an instance at N = 30 and l = 2.28. Fixing a random initial condition for all s i , i ≥ 3 we vary only s 1 and s 2 on a 200 × 400 grid and colour the points according to three different solutions they flow to. Successive magnifications illustrate the Wada property: the points on the basin boundaries are simultaneously on the boundary of all three basins implying that large enough magnifications will contain all three colours (although the blue-green boundary seems void of red in the third panel, panels four and five show that red is actually present). We conjecture that the power-law scaling of the continuous search times (3) is due in part to a generic property of the dynamical system (1-2), namely that it is hyperbolic [4, 6, 13], or near- hyperbolic. It has been shown that for hyperbolic systems the trajectories escape from regions far away from the attractors to the attractors at an exponential rate, for almost all initial conditions [4, 6, 13] . That is, the fraction of trajectories still searching for a solution after time t decays as e −κt ( Supplementary Fig. 13 ), where κ is the escape rate. Thus, κ −1 can be considered as a measure of hardness for a given formula. When taken over an ensemble at a given α, this property generates the exponential decay for p(t) with an average escape rate λ. The form of the energy function V incorporates the influence of all the clauses at all times, and in this sense the system (1-2) is a non-local search algorithm. As shown before, the auxiliary variables can be eliminated, however, they give a convenient interpretation of the dynamics. Namely, one can think of them as providing extra dimensions along which the trajectories escape from local wells, and their form (2) provides positive feedback that guarantees their escape. Clearly, these equations are not unique, and other forms based on the same principles may work just as well. To simulate (1-2), we use a 5-th order adaptive Cash-Karp Runge-Kutta method with monitoring of local truncation error to ensure accuracy. In order to keep the numerical trajectory within a tube of small, preset thickness around the true analog trajectory in Ω (Supplementary Fig. 9 ), the RK algorithm occasionally performs an exponentially large number of discretization steps n step . However, this only happens for hard formulae, when the analog trajectory has wild, chaotic fluctuations. For easy formulae both p(t) and p(n step ) decay exponentially as shown in Supplementary  Fig. 10a , inset.
For example, Ercsey-Ravasz and Toroczkai added exponentially growing auxiliary terms to each constraint, so that the network is able to escape from traps and converge to the k-SAT solution, possibly after a period of chaotic behavior REF .
Optimization hardness as transient chaos in an analog approach to constraint satisfaction
Nature Physics, vol. 7, p. 966-970, 2011
Physics Mathematics Computer Science