Asymmetric Continuous-Time Neural Networks without Local Traps for Solving Constraint Satisfaction Problems

There has been a long history of using neural networks for combinatorial optimization and constraint satisfaction problems. Symmetric Hopfield networks and similar approaches use steepest descent dynamics, and they always converge to the closest local minimum of the energy landscape. For finding global minima additional parameter-sensitive techniques are used, such as classical simulated annealing or the so-called chaotic simulated annealing, which induces chaotic dynamics by addition of extra terms to the energy landscape. Here we show that asymmetric continuous-time neural networks can solve constraint satisfaction problems without getting trapped in non-solution attractors. We concentrate on a model solving Boolean satisfiability (k-SAT), which is a quintessential NP-complete problem. There is a one-to-one correspondence between the stable fixed points of the neural network and the k-SAT solutions and we present numerical evidence that limit cycles may also be avoided by appropriately choosing the parameters of the model. This optimal parameter region is fairly independent of the size and hardness of instances, this way parameters can be chosen independently of the properties of problems and no tuning is required during the dynamical process. The model is similar to cellular neural networks already used in CNN computers. On an analog device solving a SAT problem would take a single operation: the connection weights are determined by the k-SAT instance and starting from any initial condition the system searches until finding a solution. In this new approach transient chaotic behavior appears as a natural consequence of optimization hardness and not as an externally induced effect.


Introduction
The most common approach in non-conventional computation is to treat dynamical systems as algorithms. Physical and biological systems are capable of achieving their functions and reaching their optimal state with incredible speed. Computer science and information technology tries to learn from nature and especially now, when CMOS technology reaches its limits at the small scale [1] (e.g. [2]), there is a hastened search for novel computational paradigms.
The use of analog dynamical systems for computation received increasing interest in the last three decades both in the theoretical and in engineering communities. Differential equations, continuous maps, and several neural network models have been employed to perform various computational tasks. In this approach, the dynamical systems are designed in a way to converge to attractors that are interpreted as the output of the computation [3][4][5][6]. Siegelmann, Orponen, Moore and others have recently provided a computational complexity theory for analog systems [5][6][7][8][9][10][11]. A fundamental discovery also by Siegelmann was to show, that in principle, computation beyond the Turing limit is possible. She used the strongly chaotic analog shift map as an example, and proved that it has computational power beyond the Turing machine (super-Turing) [7]. Technology has also developed devices imitating nervous system-like processing, such as the Cellular Neural/Nonlinear Network (CNN) [12,13], or analog VLSI devices [14,15]. These can solve a large variety of problems in robotics, sensory computing (vision, hearing, bionic eyeglasses) etc. The CNN is an array of analog dynamical cells performing parallel continuous-time processing, effectively solving a system of coupled ordinary differential equations (ODEs) with programmable coupling parameters. Roska and Chua proved the CNN to be at least as universal as a Turing machine [16].
Neural network models have originally been developed and investigated with the purpose of modeling brain function, however their capability of solving optimization problems has also been explored. One of the earliest works was presented by Hopfield and Tank who used neural network models to solve the traveling salesman problem [17][18][19]. Despite the evidently chaotic nature of brain activity and the theoretical results showing the power of chaotic analog dynamical systems [7], neural networks designed to solve combinatorial optimization problems mainly avoided chaotic dynamics, focusing on simple converging systems with Lyapunov dynamics (symmetric Hopfield networks, symmetric CNN etc.). In this approach the neural network minimizes a Lyapunov function (energy function) by converging directly to a local minimum [18,19], and this analog process is used as a basic step of the algorithm. However, for finding the global minimum classical techniques typical in digital computing are required (such as simulated annealing etc.) [8,9,[20][21][22], and the algorithm becomes quite estranged from the original purpose of analog computing.
The need for more complex dynamics has been realized and some chaotic neural network models solving optimization problems were presented by Chen and Aihara [23,24]. These are also called as chaotic simulated annealing methods. In this approach usually the discrete-time symmetric Hopfield network is used and local traps are avoided by introducing a deterministic chaotic dynamics with a bifurcation parameter that is gradually decreased during the annealing process [23]. This method has been further improved in different ways [25][26][27], however the need for careful tuning of parameters has not been eliminated and direct correspondence between global optimum and the final output has not been achieved.
In our recent papers [28,29] we have shown that optimization hardness is strongly interrelated with chaotic/turbulent dynamics, implying that designing analog dynamical systems with an output that corresponds directly to the solutions of a hard problem (global minima of the energy) will necessarily show transiently chaotic behavior. Here we show that asymmetric continuous-time neural networks can be designed to solve hard problems simply due to their structure, without requiring a step-by-step algorithm similar to those used in digital computers. In this case transient chaotic behavior appears as a natural consequence of optimization hardness and not as an artificially added tool.
At a recent conference [30] we presented a continuous-time asymmetric neural network (CTANN) model designed to solve Boolean satisfiability (k-SAT) (description of the model provided below). k-SAT is one of the most studied constraint satisfaction problems lying at the basis of many decision, scheduling, errorcorrection and bio-computational applications. Our model can be transformed to solve a large variety of constraint satisfaction problems, because k-SAT is NP-complete, meaning that every problem in NP can be transformed into this form in polynomial time (as function of the system size) [31,32]. The NP class contains the set of optimization problems whose solutions (once given) are easily checked to satisfy the constraints, however, finding those solutions in case of hardest problems takes exponentially long search-times. For details on the computational complexity of NPcomplete problems see [32].
Here we explore in details the properties of this CTANN model. We show how we can achieve a one-to-one correspondence between the k-SAT solutions and the stable fixed points of our CTANN. Simulations on 3-SAT, 4-SAT and 5-SAT problems show that the two important parameters do not need careful tuning during the computational process. For a given k their optimal values are fairly independent of the size and other properties of the system. We even find a common area when comparing the optimal regions for different values of k. This way non-solution traps (such as limit cycles) can be avoided and after a transiently chaotic phase the system converges to a solution.

Results
This section is organized as follows. First we briefly introduce k-SAT and summarize previous analog approaches including the continuous-time dynamical system introduced in [28]. Next we present our CTANN model discussing its key mathematical properties and finally we present numerical evidence on the effectiveness of the model.

Boolean satisfiability problems
In k -SAT there are given N Boolean variables, x i [f0,1g and a propositional formula F , which is the conjunction (AND) of M clauses (constraints) C m . Each clause is the disjunction (OR, denoted by _) of k variables (x i ) or their negation (x i ). In 3-SAT a constraint could be for example C 1~x1 _ x 5 _ x 7 . The formula may be encoded as a matrix c mi : where m~1, . . . ,M, i~1, . . . ,N. The goal is to find an assignment of the variables such that all clauses are satisfied (TRUE).
Performance of algorithms is usually tested on random k-SAT instances, where each clause includes a randomly selected set of k variables. Considering that each variable could be included in its normal or negated form, these k variables can form 2 k possible clauses. We always randomly choose one of these. The simplest measure to characterize hardness of random k-SAT formulae is the constraint density: a~M=N. In the easy-SAT region there are few constraints/clauses (small a) and it is easy to find solutions. For too many constraints (large a) it is easy to decide that the formula is unsatisfiable (UNSAT). There is an intermediate range (hard-SAT), however, where deciding satisfiability can be very hard: the worst-case complexity of any known algorithm for k-SAT (k §3) is exponential in N [31]. Using statistical physics methods it has been shown that changing the constraint density a the solution space goes through several phase transitions. The hardness of problems is related to these different phases and the hardest instances appear just before the satisfiability threshold [33][34][35].

Previous analog approaches
There have been several attempts to solve k-SAT by mapping the Boolean variables onto a continuous space. In [36,37] k-SAT is formulated as a global optimization problem in the x[½0,1 N continuous space, however it is solved with various local search and backtracking methods characteristic to digital computing (not continuous-time dynamical systems). Refs. [38][39][40] use Lagrange programming neural networks for k-SAT. They define a Lagrange function using a linear combination of the individual constraints, with coefficients serving as Lagrange multipliers. The employed ODEs are similar to a Hopfield neural network, whose trajectories however, cannot guarantee that the corresponding dynamics does not get trapped by non-solution fixed points.
In [28] we provided a new and exact mapping of Boolean satisfiability into a set of ODEs with a unique correspondence between its set of attractors and the k-SAT solutions. This eliminates the key weaknesses of previous attempts. The Boolean variables are mapped to the s[½{1,1 N continuous space. s i [½{1,1, i~1, . . . ,N, such that s i~{ 1 if the i th Boolean variable (x i ) in the SAT problem is 0 (FALSE) and s i~1 when it is 1 (TRUE). Each constraint can be formulated as a function K m (s)[½0,1 which is 0 if and only if the constraint is satisfied. Accordingly an energy function can be defined as Finding a solution to the SAT problem (if it exists) is equivalent to finding the global minima, s Ã , of this function (E(s Ã )~0). A dynamical system defined via e.g., a simple gradient descent to find the global minimum of E(s), however, will typically get trapped in local minima where E(s)w0. The continuous trajectories approach these attracting non-solution fixed-points at an exponential rate (the vector field is analytic) and hence, a corresponding exponential extraction is needed from these regions by an algorithm that does not get stuck. To achieve that, we modified the energy function by introducing auxiliary variables for each constraint (a m [½1,?),m~1, . . . ,M), acting similar to Lagrange multipliers: V (s,a)~P M m~1 a m K 2 m (s). The dynamics of s is defined as a gradient descent on the energy surface and the role of the auxiliary variables is to provide extra dimensions along which the trajectory escapes from local wells. The dynamics ensures that whenever a constraint is not satisfied, the respective auxiliary variable grows exponentially, modifying the energy function and ultimately extracting the trajectory from the local minima/wells [28]. Due to the unbounded auxiliary variables exhibiting exponential growth when needed, this system achieves polynomial continuous-time efficiency, however at the cost of exponentially large fluctuations in the energy function V (s,a). This study has also shown that the hardness of (solvable) problems appears as chaotic dynamics, however, it is of transient type [41,42] as the system still finds the solution.
While using unbounded auxiliary variables one can avoid local traps and achieves polynomial efficiency in the analog search times, the question is whether one can design a continuous-time dynamical system for k-SAT using only bounded variables (implementation friendly), but preserving as many of the desirable features of the system as possible. At a recent conference [30] we presented an implementation friendly model for solving k-SAT (see below), which is a cellular neural network model similar to those used in CNN computers and also similar to Hopfield models. Here we rigorously define the parameter region where a one-to-one correspondence between solutions and fixed-points can be achieved. We also show that the optimal region of parameters (where the system is most efficient) is independent of the properties of the problem. Most importantly, our system does not get trapped in local minima, finding the solution is one single continuous-time process which does not need ''intervention'' and tuning of parameters.

Continuous-time asymmetric neural network for k-SAT
Continuous-time recurrent neural networks in general are defined as: where x i is the state value, or activation potential of the cell, f (x) is the output function of the neuron (usually a sigmoid), u i is the input, or bias of the neuron and w ij are connection weights. Cellular neural networks have the same form, however in real implementations the cells are placed on a square lattice, and so far only neighbors can influence each other. We defined our continuous-time asymmetric neural network model on a bipartite graph with two types of nodes (cells) (Fig. 1A) [30]. One type (''s-type'') represents the variables of k-SAT, whose state value will be denoted by s i , i~1, . . . ,N and their output function defined via (see Fig. 1B): When the Boolean variable is true (x i~1 ) we will have f (s i )~1, when it is false (x i~0 ) then f (s i )~{1. However, during the dynamics we allow any continuous value f (s i )[½{1,1. For simplicity we say that f (s) is a solution of k-SAT, whenever x~½f (s)z1=2 is a solution. The input of these cells will not be needed, we fix them as u i~0 Vi. The self-coupling parameter will be a constant value w ii~A , this being one of the important parameters of the model. The clauses are represented by the second type of cells with state value a m , m~1, . . . ,M and output function (Fig. 1C): These variables, or ''a-type'' cells will play a similar role as the auxiliary variables in [28,29]. They determine the impact a clause has at a given moment on the dynamics of the s variables. For this reason g(a m )~0 will correspond to the clause being true, and g(a m )~1 to the clause being false. The second important parameter of the model is the self-coupling of these cells w mm~B . Their input will be u m~u~1 {k where k represents the number of variables in the clause (k~3 for 3-SAT, k~4 for 4-SAT, etc.). As we will see later, this is needed in order to achieve the correspondence between k-SAT solutions and stable fixed points. The connection weights between the cells are determined by the c mi matrix elements of the given k-SAT problem. The dynamical system is defined via: This neural network is asymmetrical: the influence of a clause on a variable (with connection weight c mi ) is exactly the opposite of the influence of the variable on the same clause (with weight {c mi ). We cannot assign a Lyapunov function (or energy function) to this system, the dynamics is not a simple gradient descent, and in case of hard problems it can show complex chaotic dynamics.

Important theorems
Here we list some important theorems showing the properties of the model. The proofs are presented at the end of the paper (in section Proof of Theorems).
Theorem 1 Variables remain bounded: If initially Ds i (0)Dƒ1 and 0ƒa m (0)ƒ1, then the state values of cells s i (t) and a m (t) remain bounded for all tw0, Vi,m: Theorem 2 Every k-SAT solution has a corresponding stable fixed point: Given a k- When proving this third theorem it was relatively easy to see that it holds for 1vBv2 [30]. Simulation results, however, indicate that the optimal value of parameter B is larger than 2. Indeed, here we present a more general proof showing that the correspondence between stable fixed points and solutions is preserved in the larger interval 2ƒBv2½ k 2 z2 (see Proof of Theorems).

Numerical results
The theorems presented above guarantee that all stable fixed points of the system correspond to k-SAT solutions. However, there is no guarantee that there are no other attractors -such as limit cycles or chaotic attractors -in the system. The existence or non-existence of such attractors is very difficult to show analytically, but simulation results indicate that parameters (A,B) have an optimal region fairly independent of the properties of the problems, where the dynamics avoids getting trapped in nonsolution attractors and finds a k-SAT solution.
We performed the simulations using the fifth-order adaptive Runge-Kutta method. In Figure 2 we plot the time evolution of a few s-type and a-type variables and the energy function E(f (s))~P M m~1 K 2 m (f (s)) mentioned above (for details see [28]) for a large 3-SAT problem with N~1000 variables and a~4 constraint density. While our neural network does not explicitly use an energy function (like Hopfield networks do), we use this function to monitor the evolution of the trajectory in its search for a solution. This strongly depends on the parameters (A,B). In Figs. 2A, B, C we show a case when there is transient chaotic dynamics, but finally a solution is found (A~1:54,B~2:18). In spite of the fluctuations the energy consistently decreases until finding the solution where E~0. There are parameter values A and B, however, where the dynamics gets trapped in limit cycles, with an example shown in Figs. 2D, E, F (A~1:1,B~1:1). In such cases some of the s i and a m variables remain constant but others follow complicated periodic orbits. The energy is very noisy in the simulation, but we can see it gets trapped and fluctuates in a narrow interval. This usually happens when the A and B values are small (see also below). In these situations the dynamics gets out very slowly from the subspace (s,a)[½{1,1 N S ½ 0,1 M , inside which the dynamics is linear and limit cycles can easily occur. These are not necessarily stable limit cycles, it can happen that the dynamics escapes after a very long time. Similar phenomenon of extremely long transient oscillations have been observed in CNN systems with a particular ring shape [43].
We needed to investigate how the efficiency of the system in finding solutions depends on the (A,B) parameters. As mentioned above the standard way of testing algorithmic performance is to use random SAT instances. In Figs. 3, 4 we show maps covering the A[(1,2), B[(1,3) parameter region, depicting the performance of the system and how it changes as the system size varies. We ran 100 random 3-SAT (Fig. 3) and 4-SAT problems (Fig. 4) for each point of the maps (out of the randomly generated instances we use only the satisfiable ones). The resolution of the maps is 0.02. This means that preparing these maps is equivalent to solving 50|100|100~5|10 5 instances for each N. Because these are computationally costly, we had to use relatively small instances. Constraint densities from the hardest regions of 3-SAT (a~4:25) and 4-SAT (a~9:55) were used [33][34][35]. In the first column of Figs. 3, 4 the color indicates the fraction of solved problems in the given time t max~5 000, . . . ,15000. In the second column we show the average continuous-time (not the simulation running time) the system takes to solve them (see color bars). When the solution is not found we include t max in the average. The maps show a peculiar shape consistent while changing the system size. The bottom left corner and the top middle is a parameter region where solutions are hard to find. Our observations also indicate that this is caused by limit cycles in the bottom left corner, and extremely long chaotic transients (super-transients) at the top of the map. In the middle, however there is a large region where the solution is found efficiently. As the system size increases this middle region gets lighter in the first (red) column and darker in the second one (blue). This is because the time needed to solve problems increases with the system size. In order to solve the same fraction of problems (to achieve the same red shade on the map) we would need to greatly increase the simulation times for the larger systems (see also Fig. 7), which is too costly in our case. The statistics being based on only 100 random instances, instead of searching for the optimal (A,B) parameter setting we indicate the 4% of the whole map (orange squares) where the fraction of solved problems is the largest (when the fractions are identical comparisons are made based on time values). This optimal region is fairly consistent while changing the size of problems (also see Fig. 6A, B).
We also checked how the map and the optimal parameter region changes as varying the constraint density a, and by this the hardness of problems. In the first row of Fig. 5 we show the maps for 3-SAT problems with N~40 and a~3:5,4:0,4:25 (from left to right). In the last column the frames of the optimal regions of the three maps are placed on each other showing an excellent match. Maps in the second row show results obtained on 4-SAT problems with N~20 and a~8:5,9:0,9:55 and in the third row on 5-SAT problems with N~20 and a~15,18,20:8 (the hardest region in 5-SAT being around a~20:8 [34]). Again the optimal areas show a good match. Because for larger k the B parameter can have larger values (see Theorem 3), for 5-SAT we show the maps on the B [(1,4) interval. However, we see that the optimal region still remains at the lower values close to the optimal parameter range of 3-SAT and 4-SAT.
These maps indicate that the optimal parameter region is surprisingly consistent while changing the size (N) and hardness (a) of problems. On Fig. 6A, B, C we draw on top of each other the frames of these optimal areas found on maps of all simulations performed on 3-SAT, 4-SAT and 5-SAT problems, respectively. On Fig. 6D we compare 3-SAT (black frame), 4-SAT (red) and 5-SAT (green) by showing the optimal regions in case of N~20 and constraint densities a~4:25,9:55,20:8 these being the hardest SAT phases for each. This indicates that there is a smaller region (in the middle) which seems to be part of the optimal regions of all k.
Because the variables remain bounded and there is no extra energy introduced into the system (contrary to the system in [28]) the dynamics naturally has an exponential continuous-time complexity in the hard-SAT region. In Fig. 7A, B, C we plot the fraction p(t) of problems which remain unsolved after a time t for various sizes (N~20,30, . . . ,100) of randomly chosen 3-SAT, 4-SAT and 5-SAT instances with constraint densities in the hardest a regions. We chose a parameter setting (A~1:4,B~2:24) from the common part of the optimal regions shown in Fig. 6D. The distributions are decreasing as a power law (p(t)*t {b(k,N) ), where the power b(k,N) depends on k and the size N of the k-SAT instances. b(N) is again a power law, b(N)*N {c (c%0:95,1:76,2:23 for k~3,4,5) indicating that the time complexity of the model is exponential for solving a fixed fraction of problems (Fig. 7D). This power-law decrease of p(t) shows that the probability of not finding the solution goes to zero (not a positive constant), supporting the claim that in this optimal parameter region (common for all k) the dynamics does not get trapped in limit cycles and a solution is always found after a transiently chaotic period.

Discussion
Solving NP-complete problems is a key test for any nonconventional computation. Here we presented an asymmetric continuous-time neural network that can efficiently solve Boolean satisfiability without getting trapped in non-solution attractors, and without requiring careful parameter tuning during the dynamical process. In particular, it has the following key properties: 1) It has a deterministic continuous-time dynamics. 2) All variables remain bounded.
3) The dynamics can be implemented with analog circuits (has almost the same form as used in CNN computers). 4) There is a one-to-one correspondence between the stable fixedpoints of the system and the solutions of the k-SAT problem.
Numerical simulations show that our method works consistently and efficiently on 3-SAT, 4-SAT and 5-SAT problems. A careful study of the dynamics as function of the two important parameters of the system shows that their optimal interval has a peculiar shape surprisingly consistent when changing the size and hardness of SAT instances. Comparing the optimal parameter regions for the different k~3,4,5 SAT classes we find a common parameter range which seems to work efficiently for each k-SAT instance. This assures that the system does not need careful choosing of parameters depending on the properties of SAT formulae.
While there are parameter intervals where limit cycles frequently occur (mainly the smaller values of A and B), statistics done with a parameter setting (the same for all k) chosen from the optimal region shows that here the dynamics does not get trapped in long cycles. The distribution of transient times shows a clear power-law decay in contrast with the distributions obtained with non-optimal parameters (not shown on the figures), where it goes to a positive constant value indicating that a part of problems are not solved because the dynamics gets trapped in long oscillations.
Previous approaches mainly concentrated on symmetric -and dominantly on discrete-time -neural networks. Most of the time gradient descent dynamics was used possibly combined with annealing processes. Here we have shown that asymmetric continuoustime neural networks can be designed to solve constraint satisfaction problems on their own, without additional annealing processes. This dynamics does not get trapped in non-solution attractors and transient chaotic behavior appears as an unavoidable byproduct of optimization hardness [28,29].
On an analog device this algorithm would take a single operation: the connection weights are based on the c mi matrix corresponding to the given k-SAT instance (the input of the operation) and starting from any initial condition the system searches until finding a solution, without the need of any further intervention by the user.
Our model is implementation friendly, being similar to neural networks used in analog CNN computers. However, when considering analog computation, an important question -which needs to be investigated -is the effect of noise on the dynamics. Preliminary studies show that similarly to other transiently chaotic systems [41,42], the p(t) distribution of transient times (and thus the efficiency of finding the solution) is not sensitive to noise. Actually noise effects may even help avoiding long transient oscillations, thus extending the optimal parameter region.

Proof of Theorem 1
Let us recall the dynamics of s i (t): This is a first-order ODE and its solution can be written as: It follows that where we used the facts that f (s i ) j jƒ1, 0ƒg(a m )ƒ1 and the initial condition s i (0) j jƒ1. It is also easy to see that 0ƒ Ð t 0 e {(t{t) dtƒ1. For proving Eq. (8) we recall the dynamical equation: which has the solution: Because there are exactly k variables in a clause and f (s i ) j jƒ1, then P i c mi f (s i (t)) j j ƒk. Recall also that g(a m )ƒ1, and the initial condition a m (0)ƒ1. Using this we can write: For proving the lhs of Eq. (8) we will use that Bw0, g(a m ) §0, a m (0) §0, P i c mi f (s i (t)) ƒk and 0ƒ Ð t 0 e {(t{t) dtƒ1: Proof of Theorem 2 We presented this proof in the conference paper [30], however we briefly recall it here to make it easier to readers to follow the next proof.
Given the definitions, it follows that if f (s Ã i ) satisfies the clause C m then c mi f (s Ã i )~1 and if it does not satisfy it, when c mi =0, then c mi f (s Ã i )~{1. Accordingly, the sum P i c mi f (s Ã i ) in Eq. (9) can take kz1 possible values: {k, {kz2, . . . ,k{2, k. Only the value of 2k corresponds to the clause C m not being satisfied, in all other cases there is at least one variable satisfying the constraint.
Because by assumption f (s Ã i ) is a solution, we must have: Using Eqs. (9,24) and including the values of s Ã i , a Ã m , g(a Ã m ) into the dynamical equations (5) and (6) we get ds Ã i =dt~0 and da Ã m =dt~0, confirming that we have a fixed point.
To prove stability, we will show the following: Starting in any point (s,a)~(s Ã zE,a Ã zd) in a compact vicinity of the fixed point, such that Ds i {s Ã i D~DE i DvEvA{1 Vi and Da m {a Ã m D~Dd m Dvdv1 Vm, the square distance from the fixed point: is decreasing, that is dR(s(t),a(t))=dtv0, until the dynamics reaches the fixed point where dR(s(t),a(t))=dt~0.
From Eq. (9) we know that Ds Ã i D~ADf (s Ã i )D~Aw1 Vi. We have two cases: Because the distance from the fixed point cannot increase along any of the axes: dE 2 i =dt~{E 2 i , dd 2 i =dt~{d 2 i , the conditions set for E i , d i remain valid for all i,m and the distance continues to

Proof of Theorem 3
It is easy to see that similarly to CNN models our system has the following property (the proof for CNN can be found in [12]): if Aw1 then in a stable fixed point Ds Ã i Dw1 (Df (s Ã i )D~1) Vi; if Bw1 then a Ã m v0 or a Ã m w1 (g(a Ã m )[f0,1g) Vm. (If these conditions do not hold, there is always an unstable direction along which the dynamics can escape from the fixed point, see [12].) Being in a stable fixed point s Ã : Multiplying Eq. (27) with f (s Ã i ) we get: g(a Ã n )~0u where we used Aƒ2 and introduced the notation J (z) i and J ({) i for the two parts of the sum, the first (second) part includes the clauses which are satisfied (not satisfied) by variable s Ã i . As discussed in the previous theorem, for an unsatisfied clause C m the sum P i c mi f (s Ã i )~{k. Inserting this into the dynamical equation (6) if in the fixed point we have an unsatisfied clause C m , then: Because g(a Ã m ) must be 0 or 1 and Bw1, this can hold if and only if g(a Ã m )~1 and a Ã m~B z1 for unsatisfied constraints. When parameter Bv2 it can be shown that contrary to unsatisfied clauses, satisfied clauses must have g(a Ã m )~0 (see [30]). However, this second statement is no longer true when B can have values larger than 2: there can be satisfied constraints for which g(a Ã m )~1. So let us denote as M (q) (q[Z,0ƒqƒk), the number of clauses for which g(a Ã m )~1 and there are exactly q variables satisfying the clause. (Here M (0) is exactly the number of unsatisfied constraints.) If there are M (q) clauses with g(a Ã m )~1 and satisfied by exactly q variables, from the definitions introduced in Eq. (28) we get: g(a Ã m )~X