Numerical Solution-Space Analysis of Satisfiability Problems

Numerical Solution-Space Analysis of Satisfiability Problems

Alexander Mann II. Institute of Physics, University of Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany    A.K. Hartmann Institute of Physics, University of Oldenburg, Carl von Ossietzky Straße 9-11, 26129 Oldenburg, Germany
July 15, 2019

The solution-space structure of the 3-Satisfiability Problem (3-SAT) is studied as a function of the control parameter (ratio of number of clauses to the number of variables) using numerical simulations. For this purpose, one has to sample the solution space with uniform weight. It is shown here that standard stochastic local-search (SLS) algorithms like “ASAT” and “MCMCMC” (also known as “parallel tempering”) exhibit a sampling bias. Nevertheless, unbiased samples of solutions can be obtained using the “ballistic-networking approach”, which is introduced here. It is a generalization of “ballistic search” methods and yields also a cluster structure of the solution space.

As application, solutions of 3-SAT instances are generated using ASAT plus ballistic networking. The numerical results are compatible with a previous analytic prediction of a simple solution-space structure for small values of and a transition to a clustered phase at , where the solution space breaks up into several non-negligible clusters. Furthermore, in the thermodynamic limit there are, for values of close to the SAT-UNSAT transition , always clusters without any frozen variables. This may explain why some SLS algorithms are able to solve very large 3-SAT instances close to the SAT-UNSAT transition.

satisfiability, combinatorial optimization problems, complexity, hierarchy, clusters, algorithms, stochastic local search, phase transition, Monte Carlo simulations, replica-symmetry breaking
89.75.Fb, 89.20.Ff, 75.40.Mg, 75.10.Nr, 64.60.De,05.90.+m

I Introduction

The application of notions, analytical approaches and numerical algorithms from statistical mechanics has lead to a better understanding Mertens (2002); Hartmann and Weigt (2005); Mézard and Montanari (2007) of NP-hard optimization problems Garey and Johnson (1979); Papadimitriou (1994). One main underlying question is, why these optimization problems are computationally hard. This means no fast algorithms are available, where the running times increase only polynomially with the problem size. The progress of gaining insight into this phenomenon has been considerable in particular for the typical-case complexity, where ensembles of random instances are studied as a function of control parameters. These ensembles often exhibit phase transitions where changes of the effective “hardness” of the problem can be observed. Often, these transitions are connected to changes of the structure of the solution space, comparable to energy landscapes in physics. In particular, one is interested in the question, how the change of the solution-space structures has an influence on the performance of exact and stochastic algorithms. For example, for the vertex-cover problem, which is defined on graphs, a clustering transition has been found analytically Weigt and Hartmann (2000) and numerically Barthel and Hartmann (2004); Hartmann et al. (2008) when increasing the edge density of Erdős-Rényi random graphs. This transition coincides with a change of the typical-case complexity from polynomial to exponential Bauer and Golinelli (2001). For other optimization problems, the situation is less clear, as for the satisfiability problem (SAT), which we study in this work.

As we will explain, exact enumeration of solutions works well in one region of the phase diagram, close to the SAT-UNSAT phase transition (see below), whereas Monte Carlo approaches perform well in the opposite part of the phase diagram, away from the SAT-UNSAT transitions. Unfortunately, the clustering transition is located right between these extreme parts, hence numerically difficult to study. We use a stochastic algorithm in combination with a correction of the sampling bias introduced by the stochastic algorithm to study the clustering phenomena.

The outline of the paper is as follows. In the second section, we give the necessary background on SAT and on clustering of solution landscapes. In the third section, we briefly explain the algorithms we use to sample solutions and show that they exhibit a bias. Next, we introduce ballistic networking and related methods, which we use to correct for the bias. In section five, we show the results we have obtained for random 3-SAT. Finally, we provide a conclusion and an outlook.

Ii Background

ii.1 Satisfiability

Satisfiability is one of the fundamental problems of computer science, and has attracted a lot of attention over the past years, also by physicists, due to its similarity to spin-glass problems. It is the first problem proven to belong to the class of NP-complete problems Cook (1971), a class of problems for which no algorithm has been found yet that exhibits a polynomial worst-case running time as a function of the problem size. Therefore it is still a challenge to find algorithms which perform well on typical instances and to understand the underlying structure of the solution space which may hinder the performance of algorithms.

Satisfiability belongs to the class of constraint satisfaction problems Du et al. (1997): Given Boolean variables and a Boolean formula describing a set of constraints, each of which forbids a certain assignment of values to some of the variables, you are to decide whether can be satisfied, i. e., whether there is an assignment such that all constraints are fulfilled simultaneously. In the -SAT formulation, is given in conjunctive normal form,

which describes a logical conjunction of constraints (clauses) each containing a disjunction of literals () which are either a variable or a logically negated variable .

A certain assignment of values to the variables is called a configuration in the following. If a configuration satisfies all clauses in it is called a solution. In the random -SAT ensemble each clause is chosen randomly and uniformly amongst the possible combinations in which no variable appears twice.

One defines a control parameter which is the number of clauses divided by the number of variables . For low the problem is typically satisfiable whereas for high values of there typically is no solution Mitchell et al. (1992); Selman and Kirkpatrick (1994). It has been proven rigorously Friedgut (1999) that the transition between the satisfiable and the unsatisfiable phase becomes sharp for . Whilst the position of the threshold for is known exactly Goerdt (1992), for larger there are only numerical estimates. In this paper we will stick to the case, where every clause contains exactly three literals. The satisfiability transition is located in this case at Mertens et al. (2006).

ii.2 Cluster phenomena

In addition to the SAT-UNSAT transition, analytical calculations Mézard et al. (2005); Krząkała et al. (2007) give rise to evidence that there are further (“structural”) phase transitions which refer to the formation of disconnected clusters of solutions for high values of the control parameter in the satisfiable phase. Formally, clusters in constraint satisfaction problems can be defined as extremal Gibbs measures which gives the following picture for Satisfiability: For small values of all solutions are contained in one connected component (cluster). When grows, more and more solutions disappear so that at some point the cluster decomposes into smaller clusters which initially, up to a threshold , make up only an exponentially small fraction of all solutions, whereas above many clusters contribute to the statisical behavior. Above a higher critical value we enter another type of clustered phase which is dominated by a small number of large clusters. The case of 3-SAT is special, as here , i. e., we directly enter the phase dominated by few clusters. The position of the dynamical threshold to the clustered phase is predicted to be at Krząkała et al. (2007).

This value is compatible with recent numerical results Zhou and Ma (2009), where the cluster structure was investigated using the detection of community structures. Unfortunately, the sampling was performed using an algorithm, which does not exhibit uniform sampling of the solutions, see below. Anyway, there is no general rule how to translate the formal definition of clusters, which holds in the thermodynamic limit, to finite system sizes, hence other approaches than community structures are possible. For numerical studies often a very appealing approach is used, where a cluster is defined as the connected components in a graph where each solution is represented by a vertex and edges connect solutions differing in only one variable. This definition of a cluster will be used in this work as well. For every two solutions belonging to the same cluster there is therefore a “path” of configurations which all solve the SAT instance at hand. Unfortunately, this path can be long and peppered with many dead ends or loops which makes it very difficult to decide whether two configurations belong to the same cluster. The main problem when discussing clusters in high-dimensional discrete solution spaces like that of Satisfiability is that one is tempted to think of clusters as blob-like, well-seperated and homogeneous structures in configuration phase like, e. g., nano-clusters formed by agglomeration of atoms. The clusters which occur in high-dimensional discrete solution spaces are yet of a completely different nature in that they are more like fragmented and interweaved structures with lots of dead-ends, loops and holes which makes it difficult to speak of spatially seperated clusters.

The existence of a clustered phase has been proven rigorously for Mézard et al. (2005). In the language of statistical physics this clustering corresponds to one-step replica-symmetry breaking (1-RSB) Mézard and Zecchina (2002); Mézard et al. (2002). A further substructure in terms of another clustering of solutions taken from one cluster, giving a hierarchical structure of clusters is suspected where 1-RSB becomes instable and higher steps of replica symmetry breaking occur Montanari et al. (2004).

What makes cluster phenomena interesting from the algorithmic point of view is the question if (and if so in what way) clustering has an influence on the performance of local search heuristics. Usually it is assumed, that the existence of many clusters is an indication for a complicated “rugged” energy landscape, which then also gives rise to many local minima, hindering the performance of local search heuristics Montanari et al. (2004). In the same way, but with a slightly different focus, Krząkała et al. Krząkała et al. (2007) propose that the appearance of locally frozen variables in clusters is responsible for the slowdown of heuristic algorithms close to the SAT-UNSAT threshold. A locally frozen variable is a variable which takes the same value over all solutions belonging to one cluster. A cluster containing at least one frozen variable is called frozen. One defines the freezing transition as the smallest value of above which all solutions belong to frozen clusters.

To clarify the influence of phase transitions on the average computational hardness one can study the performance of stochastic algorithms as a function of the control parameter . Of particular interest is the algorithm-dependent value of up to which an algorithm shows linear-time performance, and compare this to threshold values of Ardelius and Aurell (2006). Studies of stochastic algorithms such as ASAT Ardelius and Aurell (2006), WalkSat Seitz et al. (2005) and ChainSat Alava et al. (2008) have shown however that those algorithms have linear behaviour up to values considerably beyond the clustering transition. This suggests that the cluster transition has no impact on the performance of local search algorithms as long as there are precautions against entropic traps. It is remarkable that ChainSat has this behaviour although it is greedy “in a weak sense” as it never allows steps which increase the number of unsatisfied clauses. Naively one would therefore expect it to get trapped very easily in local minima. The authors of Alava et al. (2008) interpret this as evidence for the belief that true local minima are very rare in high-dimensional search spaces. These results could also indicate that indeed it is more the (non-) existence of frozen clusters, which is responsible for the performance of local-search algorithms.

Nevertheless, for small instances there are always some frozen variables. Therefore in Ardelius and Zdeborová (2008) a different notion of frozen clusters via the whitening core is used. There one looks, for each solution, interatively for variables which can be flipped since they appear only in clauses satisfied by other variables or which contain variables already detected in the whitening core. The position of this freezing transition was then calculated by exact enumeration and clustering of all solutions for sufficiently small system sizes and is expected to lie at close to but below the satisfiability transition. In a similar way, Ardelius et al. (2007) finds a cluster condensation transition in the solutions generated by ASAT very close to , again these results rely on a non-uniform sampling of the solutions. Anyway, these results are compatible with the observation of a good performance of local-search algorithms close to the threshold .

ii.3 Algorithmic treatment of SAT

Algorithms for SAT include a broad spectrum, both stochastic and exact, from simple and straight-forward algorithms like RandomWalksat Papadimitriou (1991) and WalkSat Selman et al. (1996); Kautz et al. (1997) to complex algorithms like DPLL Davis et al. (1962) and message-passing algorithms such as Belief Propagation (BP) and Survey Propagation (SP) Mézard et al. (2002). For small systems exact enumeration of all solutions is possible using one of the numerous standard algorithms Gomes et al. (2008); Bayardo and Pehoushek (2000) such as the aforementioned DPLL. It can be shown Cocco and Monasson. (2001) that deterministic algorithms have longest average run times close to reflecting the difficulty of deciding whether a given SAT formula is satisfiable or not. The problem with exact enumeration is that it is limited to small systems due to hardware restrictions, especially because of the memory needed to store the huge number of configurations, as the number of solutions grows exponentially with system size. Furthermore the number of solutions is not a continuous function when crossing the satisfiability threshold, but it drops from a finite number to zero. This corresponds to a non-zero entropy at the phase transition, the entropy per variable grows approximately linearly with decreasing Monasson and Zecchina (1996); Biroli et al. (2000). In turn this means that even very close to the satisfiability threshold the number of solution grows exponentially, and quickly becomes so large that it is not feasible to enumerate all solutions even in this regime. From counting all solutions using DPLL for systems up to we can estimate the solution entropy near the phase boundary to be roughly between and which gives more than solutions already for , thus taking up at least 2 GB.

To overcome these limitations one can turn to stochastic algorithms which, starting at an arbitrary configuration, do successive changes either completely randomly or based on a heuristic evaluating information about the local configurational neighbourhood. Stochastic algorithms are not guaranteed to find a solution, even if solutions do exist, but they can be significantly faster than deterministic algorithms. It is thus possible to obtain solutions for much larger systems, but on the other hand stochastic algorithms can never prove that there is no solution, tests for unsolvability can only be done by using deterministic algorithms.

In this paper we study the cluster structure numerically for , which requires unbiased sampling of the solution space. Different types of sampling algorithms are studied and shown to be biased. We therefore present an algorithm which uses a different approach to create a survey of the cluster structure of Satisfiability instances from which it is then possible to derive unbiased samples. It is an improvement on the “Ballistic Search” algorithm which has originally been applied to spin glasses Hartmann (2000a, b); Hed et al. (2001). The main advantage of the algorithm is that it is able to provide an overview of the cluster structure of the solution space without having to enumerate all solutions which is no longer possible already for moderate numbers of variables.

Iii Sampling algorithms

iii.1 Bias in SLS algorithms

If stochastic local search (SLS) algorithms found all solutions with the same probability one could use them directly to probe the solution space. Unfortunately, this is not the case as we will see in this section. Later on in this article we will use ASAT as solution generator so we use it here for an exemplary presentation of the bias in SLS algorithms.

ASAT is a simplified variant of Focused Metropolis Search Seitz et al. (2005) and was first described in 2006 Ardelius and Aurell (2006). It starts at a random configuration and in each step picks a variable from an unsatisfied clause. This variable is flipped if either this decreases the number of unsatisfied clauses, or otherwise with a constant probability which is a tuning parameter of the algorithm. ASAT has run times linear in the system size at least up to on 3-SAT. For instances of moderate size like those we study here, it can well be used beyond this point Ardelius et al. (2007).

The test procedure is very simple: For a randomly chosen small instance we run ASAT again and again starting each time from a different randomly chosen configuration and count how often each solution returned by ASAT is found. If there were no bias we would expect the histogram of solution multiplicities to be flat except for statistical fluctuations around a plateau value, i. e. the histogram should resemble a Gaussian error function.

Fig. 1 shows the resulting histogram, in comparison to a histogram filled with the same number of random integers drawn from a truely flat distribution over the range corresponding to the number of solutions of the SAT instance showing what the distribution should look like if there were no bias. Clearly there is a strong bias favouring some solutions over others. To quantify the deviations we use a test, and calculate the -value giving the probability that an unbiased sampling process yielded a sample deviating at least as much than the one at hand. The -values numerically are smaller than (i. e., the resolution of our double numbers).

To test whether this bias can be corrected in a simple way, we did a further check, where instead of using the solutions returned by ASAT directly, for each solution found by ASAT a solution from the same cluster was generated using a Monte Carlo (MC) search starting at the ASAT result. The outcome of this modification is shown in the inset of Fig. 1 for the same SAT instance as before. The distribution now clearly has 5 plateaus corresponding to the 5 clusters of the solution space and looks much flatter but exhibits still some bias. One sees that most of the ASAT solutions stem from the smallest cluster, hence the sampling does not respect the cluster size. Hence, ths bias can be strongly decreased by additional MC simulations, but not completely. Further checks showed that the bias persists independently of the system size.

Figure 1: Multiplicities of solutions found by ASAT in runs for a randomly chosen instance with , , compared to an unbiased distribution. Inset: Multiplicities of the ASAT solutions after an additional -MC step with 10 sweeps.

Since we want to study clustering properties of the solution ensemble we need to remove the bias completely and sample solutions in proportion to the cluster sizes. To ensure this, we will perform reweighting using the Ballistic-Search algorithm as described in section IV. Before we come to the Ballistic Search, we will show in the next section that MCMCMC, another important sampling method, fails on sampling SAT solutions uniformly as well.

iii.2 Bias in MCMCMC

The Metropolis-Coupled Markov Chain Monte Carlo (MCMCMC) method, first proposed in 1991 by Geyer Geyer (1991), also known as Parallel Tempering Hukushima and Nemoto (1996); Newman and Barkema (1999), is a powerful and versatile tool, commonly used in biophysics and statistical physics to perform equilibrium simulations and to generate unbiased samples in large configuration spaces. MCMCMC uses a set of replicas of single instances, simulated in parallel at different temperatures and linked by global updates in which replicas are swapped pair-wise with an acceptance probability depending on their energy difference and temperature spacing (Metropolis-Hastings criterion), thus facilitating the tunneling through barriers seperating local minima of the phase space Moreno et al. (2003).

To study the performance of MCMCMC on SAT we employ a histogram test similar to the one described in section III.1 for the performance of ASAT. For several values of scattered over the satisfiable phase, the number of variables is chosen such that expected number of solutions is 1000. This, e. g., results for the smallest value of considered here in a system size , while for the highest value of , is feasible.

We apply a straight-forward implementation of MCMCMC to a set of 50 instances for each value of the control parameter , where we use 15 temperatures, the lowest, at which the samples are taken, being initially , the highest such that the corresponding energy is found to be approximately which is the expected energy of a completely random configuration. Every 1000 steps the temperatures are adjusted to drive the replica exchange rate between neighboring temperatures towards 50 lowest temperature fixed at . The procedure chosen to adjust the temperatures leads to a distribution of temperatures where for the lowest temperatures the exchange rates indeed reach 50 whereas the highest temperatures all gather in the random phase. This can be seen as an indication that the number of temperatures used is sufficient to allow the replicas to travel between the constrains, i. e. the highest temperatures are indeed located in the “paramagnetic phase”. We take one sample every second sweep to generate a total of samples. Only successful sampling steps are counted, i. e. those where the energy of the configuration at is zero.

The histograms with the resulting distribution look pretty close to those drawn from a flat distribution (not shown here). We again use the -values obtained from a test to quantify the deviations and find that in most cases MCMCMC gives reasonably flat distributions, hence this method appears to exhibit on the first sight a much lower sampling bias. Nevertheless, there also are a number of histograms having a significant bias exhibiting very small -values. The higher the value of the control parameter is chosen, the larger becomes the spread of the distribution of -values towards extremely small values. This can be seen from Fig. 2 where the distribution of -values, integrated over all system sizes and values of the control parameter, is shown for the two cases where the instances exhibit one or more than one cluster of solutions, respectively. (The number of clusters can easily be calculated exactly for these small instances.) In the case that the solution space consists of only one connected component, the -value distribution is flat showing that MCMCMC works unbiased as expected. The presence of clustering on the other hand leads to a bias or imbalance in the sampling process resulting in a strong peak at small -values.

Figure 2: Bias of MCMCMC: Dependence of -value on the number of clusters, integrated over all and . Dashed line: -values for instances with only one cluster, solid line: -values for instances with more than one cluster.

Since we are interested in particular in those instances which exhibit many clusters, MCMCMC turns out to be not suitable as well, since all instances have to be sampled correctly. Note that for larger system sizes, the number of instances having just one cluster, where MCMCMC seems to work well, will strongly decrease. Hence, for large system sizes, MCMCMC will exhibit a bias for basically all instances of interest. To create an unbiased sample we need a different method which will be presented in section IV.

Iv Ballistic search

Here, as mentioned in Sec. II.2, we are using the neighbor-based definition of clusters: two solutions are considered to be in the same clusters, if there exists a path in solution space consisting of single-variable flips.

We use Ballistic Search, which has been introduced in the year 2000 as a method for studying ground-state properties of spin glasses Hartmann (2000a). The approach is able to provide a survey of the cluster landscape using stochastic algorithms, in particular without the need to enumerate all ground states as it is usually necessary when one aims at clustering. The sheer number of ground states forbids exact enumeration when studying spin glasses, and, as mentioned above, the same holds for SAT. We therefore use this method which relies on generating a survey of the most important clusters.

The survey consists of a set where each element represents one cluster and consists of a (small) set of solutions from the cluster and an estimate of the size of cluster . The survey should cover all clusters, or at least all but those which are negligibly small. One can then sample the whole solution space with correct weights by generating the desired number of solution samples from the representative sets of solutions for each cluster according to the respective cluster sizes.

Two main ingredients form the basis of the Ballistic Search algorithm. Firstly, the above described data structure storing small sets of representative solutions for each cluster instead all solutions. Secondly, a “ballistic path search” is used to analyse the cluster space and generate the survey from a given set of solutions. The basic operation of this procedure is that we have to determine for any given pair , of solutions, whether they belong to the same cluster. This has to work under the assumption that for the case of , belonging to the same cluster, a complete nearest-neighbor path of solutions between and is not contained in the set of already found solutions. Instead, one searches stochastically for paths between and by starting at one solution and subsequently changing a randomly chosen free variable. A variable is called free if its value can be changed without violating any constraint, so that one never leaves the solution cluster. This is repeated until either the target solution is reached or no free variable is left, because as an additional stop condition every variable shall be touched at most once. (Because of this additional constraint the path search is called “ballistic”.) This implies that in a successful ballistic path search the number of steps taken is always equal to the Hamming distance of the solutions, i. e., the number of variables in which the two solutions differ.

Figure 3: Between those 6 solutions (black circles) from the same cluster (light grey) Ballistic Search has found only 2 paths. The apparent number of clusters is 4. We need to increase the density of solutions to make Ballistic Search more efficient.

Figures 3 and 4 give a graphical description of how the Ballistic Search algorithm works. We start with some randomly generated solutions depicted as black circles in Fig. 3, all of which belong to the same cluster which is drawn in grey in a 2-dimensional cartoon of the N-dimensional configuration space. For the sake of simplicity we assume here that all solutions belong to the same cluster, the generalization to more than one cluster is obvious. Running the ballistic path search we find that some of the solutions can be connected by paths drawn as lines in the picture, i. e., for these solutions the algorithms correctly finds that they belong to the same cluster. The problem is that for low solution densities the average distance between solutions is large, and the efficiency of the ballistic-path search strongly decreases with larger distances Hartmann (2000a). Therefore we only find a few paths and the apparent number of clusters in this example is larger than its true value, cf. Fig. 3.

Figure 4: Adding solutions (grey) has yielded a correct identification of the cluster, because more paths (grey) have been found, now connecting all solutions.

What we need to do is to increase the number of solutions by rerunning ASAT. For few added solutions, the measured number of clusters will increase, since only few additional paths within clusters are detected, less than the number of added solutions. When generating even more solutions, we will find that the apparent number of clusters at some point no longer increases, but instead it decreases as more and more paths between solutions are found, until finally all solutions are correctly assigned to the same cluster as shown in Fig. 4.

iv.1 Ballistic networking

Figure 5: This cluster has a more complex structure than the one in Fig. 4, illustrated by the additional holes. Here adding even more new solutions (grey) than before does not work. Still not all solutions are recognized as belonging to the same cluster.

Our studies have shown that the simple ballistic-path search algorithm has very low efficiency when applied to Satisfiability which can be attributed to a high complexity of the solutions space or large sizes of the clusters. We therefore developed a refinement of the algorithm named “ballistic networking”, which is a very general extension of the ballistic path search so that it can readily be applied to other problems.

Figure 6: Ballistic Networking improves on the result of Ballistic Search by not adding solutions randomly, but adding solutions from the same cluster using a MC search (arrows). Now all solutions are found to belong to one cluster.

The idea of the algorithmic refinement is to increase the probability of identifying two solutions, origin and target , as belonging to the same cluster using ballistic path search by again increasing the number of solutions. Instead of using ASAT to generate more solutions we generate additional solutions by performing independent MC simulations starting at and , respectively. Hence, we are sure that the additional solutions belong to the same cluster as their respective “parent” solution. We then try to find connections using the ordinary ballistic path search between all pairs of solutions, where one solution belongs to the origin and one to the target. If at least one path is found, it is clear that and belong to the same cluster. We apply this test to all pairs of solutions which have not yet been found to belong to the same cluster. An artist’s view of this improvement is given in figures 5 and 6. Fig. 5 shows how the standard ballistic search fails due to a more complex structure of the cluster, although even more solutions than before have been used. In Fig. 6 the solutions found by the MC search are drawn as points connected to their parent solution by arrows. We can see that the number of successful ballistic path searches (gray lines) does not have to be very high, but still is enough to correctly identify the cluster.

Figure 7: Comparison of ballistic path search without and with additional solutions (“ballistic networking”) for and . We show the probability for finding a path between two solutions generated from the same cluster, as function of the total number of ballistic path searches between all pairs of parent and children configurations averaged over 1000 runs. The case “+0” corresponds to the original ballistic path search.

Indeed this procedure improves the performance of the search so much that it outweighs the additional effort of having to carry out ballistic path searches instead of one. Fig. 7 shows a comparison of the performance of ballistic path search without and with additional solutions. The case “+0” corresponds to the bare ballistic path search. The horizontal axis shows the number of ballistic path searches which have to be carried out in the worst case, i. e., when no connecting path is found. We found to be a suitable range for the system sizes under study.

When creating the additional solutions to test whether a pair of solutions , belongs to the same cluster, to improve the success probability, one can think of introducing a bias into the MC search which pushes the additional solutions derived from the first solution closer to the second solution and vice versa. Indeed we found that such a bias has a positive influence on the success probability of the ballistic path search. Yet we did not use this bias in the implementation, because the positive influence comes at a high cost. For each pair of solutions to be tested a dedicated biased set of additional solutions has to be generated which cannot be reused when comparing either or to a third solution. The necessary computational effort for generating each time new biased configurations by far outweighs the positive effect of the bias.

Figure 8: Comparison of the actual number of configurations per cluster to the number estimated by Monte Carlo integration. Each point corresponds to one cluster.

The second part of the cluster survey consists of the sizes of the clusters. An exact calculation of the cluster size is possible, but takes too long, since it typically grows exponentially with . We have therefore examined several different estimation methods with respect to their reliability in giving a correct estimate for the cluster size by comparing the estimated cluster size to the exact cluster size on a random ensemble of clusters for different values of and small system sizes . The best method known to us has been found to be the estimation of the cluster size using a Monte Carlo integration as it has been used in Hartmann and Ricci-Tersenghi (2002) in an application to spin glasses. Fig. 8 shows a comparison of the actual number of configurations in one cluster to the number as estimated by Monte Carlo integration. We used several combinations of and where the total number of solutions (over all clusters) did not exceed , such that all clusters could be calculated exactly, and afterwards for each cluster the MC estimation was run.

iv.2 Implementation

Combining ballistic networking and the cluster size estimation the full algorithm is comprised of two alternating steps. The first step is to generate a given number (of the order of 1000) of solutions of the Satisfiability instance at hand using the ASAT algorithm. The tuning parameter of ASAT is chosen to be which is the optimal value as given in Ardelius and Aurell (2006). In the second step Ballistic Networking of the solutions found by ASAT is done as described above to create the cluster survey, and then the sizes of the clusters are estimated. Afterwards ASAT is run again and another set of new solutions is created. The cluster survey is then updated using Ballistic Networking on the new solutions and the solutions representing the so far found clusters in the existing survey. Here new clusters may be found, and if so, their size is estimated. This is repeated until the cluster survey is considered complete, i. e., no more relevant clusters are found.

From the cluster survey for each instance a set of unbiased solutions can be generated using the cluster-size estimates. For each solution to be generated for a given instance, first a cluster from the survey is selected with a probability proportional to the cluster size. One solution is selected from the set of representative solutions and starting from this solution a MC search is performed finally giving the solution to be used in the analysis.

Defining a good stopping criterion is a crucial point of the algorithm. As the cluster number in Satisfiability can be rather large, we decided not to generate all clusters, but all except for those which contain only a neglegibly number of solutions. For this purpose we monitor the total cluster weight . We run the algorithm until the total cluster weight has not increased by more than over the last half of solutions included in the clustering process. We store the order in which the solutions have been generated by ASAT and label each cluster with a number telling the position of the earliest solution which has been found to belong to this cluster.

When trying to optimize the number of new solutions added in each round one has to consider two competing effects: On the one hand adding solutions — as in ordinary Ballistic Search — may reveal that two clusters actually are parts of the same cluster, connected maybe by only a small path in configuration space which has been too hard to find with fewer solutions. On the other hand increasing the number of solutions makes Ballistic Networking slower and, even worse, increases the probability of false new clusters which in turn can lead to an on-going increase of the cluster number and total cluster weight and thus to a failure of the stopping criterion.

The system sizes which can be reached using the method described above depend, of course, on the control parameter . For small all solutions are contained in only one large cluster where there are many possible paths between configurations so that Ballistic Search is very efficient and system sizes of a few hundred variables are possible. For high values in the solvable phase, the number of solutions is small so that in this regime Ballistic Search still is rather efficient due to the small extent of the clusters and relatively large system can be done.

Figure 9: Dependence of of ballistic path search on and for . As the typical distance of two solutions depends on , so do the ranges shown in the plot. For small distances , a path is usually found, while for large distances, the probability depends on the value of the control parameter , which can be understood by a more and more complex structure of the solution clusters.

Fig. 9 shows the dependence of the success probability of ballistic path search on the Hamming distance between the configurations for , for different values of the control parameter . Up to the probability decreases strongly with increasing , as the clusters develop more holes. Above this point the curve is approximately independent of . We also find that the probability decreases weakly with increasing system size (not shown). The fact that the average distance between solutions decreases with makes Ballistic Networking most difficult in the intermediate regime around . Here the number of ballistic path searchs needed to find a connection between two solutions from the same cluster is highest. The cluster structure seems to be such that there are many “dead ends” in which the search may get stuck. Together with the high number of clusters, which enters quadratically in the running time, this limits the reachable system size. All in all, Satisfiability instances of up to variables were doable in reasonable time over the whole range of interest while for smaller intervals of the control parameter, we also studied .

V Results

Here, we study the behavior of random 3-SAT instances as a function of the parameter . This is meant in the sense that we generate an instance using a given number of variables and a set of (arbitrarily ordered) clauses (). We chose , where is the largest value of the control parameter we want to consider. We can study the behavior of each instance as function of by considerung each time exactly the clauses for . Also, we can average over these distances for each value of the control parameter.

v.1 Hierarchical clustering

For the analysis Hartmann (2009) of the behavior of 3-SAT as a function of the control parameter , we start by looking at the hierarchical structure of a set of solutions sampled for a typical 3-SAT instance. We have used “Ward’s algorithm” Ward Jr. (1963); Jain and Dubes (1988), an agglomerative hierarchical matrix updating algorithm, on the set to extract a hierarchical clustering from which we can then draw a visual representation of the solution space.

Ward’s algorithm has been applied in many different fields ranging from RNA secondary structures over optimization problems to spin glasses Higgs (1996); Hed et al. (2001); Hartmann et al. (2008); Zhou and Ma (2009). It is an iterative procedure where initially each configuration comprises a single item cluster. In each step those two clusters are merged which have minimal distance with respect to an effective distance measure chosen such that the sum of the variances in each cluster is minimized. After each merger, the distances of the remaining clusters to the new cluster have to be calculated, for details see, e. g., Ref. Jain and Dubes (1988). Finally, one reorders the configurations according to the hierarchy obtained in the iterative merging process, and draws a color-coded visualization of the distance matrix.

Figure 10: The hierarchical structure resulting from Ward’s algorithm visualized both as tree structure (dendrogram) and distance matrix, for and (top), and (bottom). Darker grey scales correspond to smaller distances.

Next, we present some results for a typical instance. We chose one which exhibits its SAT-UNSAT transition close to the numerical estimate of the ensemble average given in Mertens et al. (2006). Fig. 10 shows the color-coded distance matrices and the dendrogram which were generated for three different values of . The difference in the solution landscape and cluster structure between the phases is clearly visible. For low the Ward matrix is featureless and homogeneously grey. All solutions belong to one single cluster and the phase space shows no specific features. In the intermediate range one sees box-like structures along the diagonal in a darker grey. These correspond to clusters, because darker means smaller Hamming distance and the solutions inside a cluster are closer to each other than to other solutions. Some of these boxes show a substructure which can be interpreted as the solutions from this cluster themselves forming sub-clusters. This is consistent with the theoretical prediction of replica-symmetry breaking beyond 1-RSB in the intermediate range Montanari et al. (2004). Nevertheless, as mentioned in the introduction, it is to be expected that most of the clusters are not relevant in the thermodynamic limit and a small number of clusters contains almost all solutions. For higher values of the substructures inside the clusters become washed out whereas the first-level cluster structure becomes more pronounced as the cluster become smaller. In the replica symmetry breaking framework this would be interpreted as a vanishing of higher level RSB above a certain threshold, but this cannot be deduced from looking at single instances, of course.

v.2 Averaged quantities

Figure 11: Complexity as a function of and for several system sizes between and . Lines have been drawn to guide the eye. The error bars give statistical errors.

The complexity is defined as the logarithm of the number of clusters normalized to the system size. Fig. 11 shows the complexity as a function of and for system sizes up to variables averaged over 200–500 instances for each value of and each value of . The number of clusters was taken directly from the cluster surveys created using the Ballistic Networking method described above.

For the “easy” part of the satisfiable phase, where the value of the control parameter is small, there is only one cluster, thus the complexity is zero. In an intermediate range the number of clusters grows peaking at a value which is strongly affected by finite-size effects and then becomes smaller again. This behaviour reflects the theoretical prediction of one single cluster in the low regime “crumbling” into smaller pieces when is increased and the clustered phase is reached. For even higher the vanishing of solutions leads to the disappearance of clusters and the cluster number decreases again. Clearly, the peak of the complexity curves seems to converge towards with increasing system size, in accordance with the analytic predictions Krząkała et al. (2007).

Looking at this plot one has to keep in mind that the stopping criterion used in the algorithm is based on the number of solutions covered by the clusters that were found so far. In a phase with a large number of small clusters we will miss small clusters if they only comprise a negligible part of the solutions (in the sense of the stopping criterion) and therefore underestimate the number of clusters. It is therefore natural that the complexity found here is lower than the one given in Ardelius and Zdeborová (2008). After all the complexity shown in the graph is only a lower bound for the true complexity respecting all clusters. In phases dominated by few and large clusters its value should nevertheless be close to the true value.

Fig. 12 displays the fraction of the solutions contained in the largest cluster. For this value seems to increase with growing system size. At it exhibits a minimum, while for it decreases slightly with growing system size, but it is larger than the value found right at . These results are also compatible with the analytical prediction Krząkała et al. (2007), which states that only for a range more than one cluster is relevant in the thermodynamic limit. Nevertheless, we cannot deduce from the data, since system sizes we can reach are rather limited, whether for all values of this growing fraction converges to one or to a smaller value.

Figure 12: Fraction of the weight of the largest cluster with respect to the total weight of all clusters.

Next, we have a closer look at the average structure of the solution space. As mentioned in the discussion of the Ward matrices, solutions belonging to the same cluster are more similar to each other, i. e., closer in terms of the Hamming distance, than pairs of solutions which belong to different clusters. The cluster structure is thus reflected in the set of all pairwise overlaps, where the overlap of two solutions and for which the Boolean variables take the values given by and is defined as . Fig. 13 shows the overlap distribution for and several values of . For each system size at least 1000 instances have been processed with the algorithm described in section IV.2 and 500 solutions have been generated from each cluster survey.

Two peaks are visible. One peak is lying close to and due to the overlap of solutions belonging to the same cluster. With larger system size it moves slightly to lower values and becomes sharper. The second peak at about is not discernible for the smallest system size, but only evolving with larger system sizes and only visible weakly against an also growing background. Note that a pure two-peak structure would correspond to the picture of one-step broken replica symmetry (1-RSB) Mézard et al. (1984). Nevertheless, the result is not fully clear here, since in addition to the peaks, there is also a contious part between the two peaks. On the other hand, it is clear that the overlap converges to zero for values of the overlap smaller than 0.5. This speaks against a full-level RSB structure of the solutions space. Note that we have found similar results for other values of the control parameter (not shown).

Figure 13: Overlap of solutions in clusters found by Ballistic Networking, for . For the curves are essentially zero.

v.3 Freezing transition

To complete the picture we also studied the freezing transition, which as mentioned in section II.2 is defined as the smallest above which all solutions belong to frozen clusters and has been found to lie at .

To check directly whether a cluster contains frozen variables, we need to generate and compare all solutions from this cluster, therefore cluster surveys do not help here. Using exact algorithms we find that for the system sizes we can reach, for all near the SAT-UNSAT transition there are always frozen variables in all clusters. This is probably due to too small systems sizes.

Thus, we followed a different approach. For each instance, taken at and and for system sizes up to , we generated a solution using ASAT, which belongs with high probability to the largest cluster. Then we performed a very long simulation starting from this solution, and measured the fraction of variables which have never flipped while performing this random walk inside the solution cluster. We extrapolated this fraction to a large number of MC steps, yielding , see Fig. 14. With increasing system size, seems to converge to zero, see inset of Fig. 14. Hence for and the largest clusters seems to contain no frozen variables in the thermodynamic limit. This is compatible with , meaning that in the thermodynamic limit no frozen clusters occur below this value of .

Figure 14: Fraction of never-flipped variables for the “largest” (i. e. first detected) cluster during MC simulations for and different system sizes . The solid lines indicate fits to functions of the form Inset: Extrapolated values as function of system size in a double-logarithmic scale, for and . The “?” marks a point where the extrapolation failed and an upper limit was estimated by eye.

Vi Conclusion

In this work we have shown that stochastic local search algorithms cannot be expected to produce correctly weighted samples of the solution space of Satisfiability. The same holds true for MCMCMC which is widely used to sample configuration spaces in many fields of application, when the SAT solutions are spread over several clusters.

A new type of algorithm has been presented and used for studying clustering phenomena in the solution landscape of the Satisfiability problem. It is an improved version of the Ballistic Search algorithm which has been successfully used for studying spin-glasses. Its guiding principle is to generate a survey of clusters of solutions represented by small sets of solutions rather than enumerating and clustering all solutions which is unfeasible already for moderate system sizes. By using a different approach, Ballistic Networking, in the reconstructing process of the cluster structure the efficiency of the Ballistic Search could be improved so that its performance becomes reasonably high when used on Satisfiability. The method presented here is general enough to be suitable for many other problems. Of course, it would be natural to study Satisfiability for using Ballistic Networking, but the efficiency for ballistic path search seems to be still much lower than for the case of which sets very restricting limits on the system sizes which can be reached. Nevertheless, the approach presented here should be useful for many disordered systems like other types of combinatorial optimization problems.

In the case of Satisfiability, the range of low values of (where many solutions exist but belong to only one cluster) can be studied by MCMCMC. Furthermore, the case of high values of close to the SAT-UNSAT transition (fewer solutions contained in several clusters) can be studied using exact enumeration of all solutions. In contrast, the algorithm presented here allows to study the full satisfiable phase, but it is limited to moderate system sizes in the intermediate range. Nevertheless it is the only reliable method to generate unbiased samples in this regime.

Using the method described here the ensemble properties of Satisfiability with moderate system size could be studied and analytic predictions about the cluster structure could be tested. To this aim we first did a visual inspection of the cluster landscape using a graphical representation in terms of Ward distance matrices. These show the expected structural differences of the different phases of Satisfiability. Furthermore we had a look at the complexity measure over the whole spectrum in the easy phase, and the overlap distribution of the solutions for particular values of . Our findings are in good agreement with the theoretical predictions and previous numerical studies using other methods.

We have profited a lot from discussions with M. Alava, J. Ardelius, E. Aurell, P. Kaski, F. Krzakala, M. Mézard, A. Montanari, P. Orponen, S. Seitz and L. Zdeborová. This work was supported financially by the VolkswagenStiftung within the program “Nachwuchsgruppen an Universitäten”. We furthermore wish to acknowledge the allocation of computer time by the Gesellschaft für wissenschaftliche Datenverarbeitung mbH Göttingen, by the Institute of Theoretical Physics (University of Göttingen), and by the Cluster for Scientific Computing “GOLEM” (University of Oldenburg).


  • Mertens (2002) S. Mertens, Computing in Science and Engineering 4, 31 (2002),
  • Hartmann and Weigt (2005) A. K. Hartmann and M. Weigt, Phase Transitions in Combinatorial Optimization Problems (Wiley-VCH, 2005).
  • Mézard and Montanari (2007) M. Mézard and A. Montanari, Constraint Satisfaction Networks in Physics and Computation (Oxford University Press, 2007), preprint, URL
  • Garey and Johnson (1979) M. R. Garey and D. S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness (W. H. Freeman, 1979).
  • Papadimitriou (1994) C. H. Papadimitriou, Computational Complexity (Addison-Wesley, Reading, 1994).
  • Weigt and Hartmann (2000) M. Weigt and A. Hartmann, Phys. Rev. Lett. 84, 6118 (2000).
  • Hartmann et al. (2008) A. K. Hartmann, A. Mann, and W. Radenbach, J. Phys.: Conf. Series 95, 012011 (2008).
  • Barthel and Hartmann (2004) W. Barthel and A. K. Hartmann, Phys. Rev. E 70, 066120 (2004).
  • Bauer and Golinelli (2001) M. Bauer and O. Golinelli, Eur. Phys. J. B 23, 339 (2001).
  • Cook (1971) S. A. Cook, in Conference Record of the 3rd Annual ACM Symposium on Theory of Computing (ACM, 1971), pp. 151–158.
  • Du et al. (1997) D. Du, J. Gu, and P. M. Pardalos, eds., Satisfiability Problem: Theory and Applications, vol. 35 of Dimacs Series in Discrete Mathematics and Theoretical Computer Science (American Mathematical Society, 1997), DIMACS Workshop, March 11-13, 1996.
  • Mitchell et al. (1992) D. Mitchell, B. Selman, and H. Levesque, in Proc. 10th Natl. Conf. Artif. Intell. (AAAI-92) (AAAI Press/MIT Press, Cambridge, Massachusetts, 1992), p. 440.
  • Selman and Kirkpatrick (1994) B. Selman and S. Kirkpatrick, Science 264, 1297 (1994).
  • Friedgut (1999) E. Friedgut, J. Amer. Math. Soc. 12, 1017 (1999).
  • Goerdt (1992) A. Goerdt, in Proceedings of the 17th International Symposium on Mathematical Foundations of Computer Science, MFCS’92 (Prague, Czechoslovakia, August 24-28, 1992), edited by I. M. Havel and V. Koubek (Springer-Verlag, 1992), vol. 629 of LNCS, pp. 264–274.
  • Mertens et al. (2006) S. Mertens, M. Mézard, and R. Zecchina, Rand. Struc. Algor. 28, 340 (2006).
  • Krząkała et al. (2007) F. Krząkała, A. Montanari, F. Ricci-Tersenghi, G. Semerjian, and L. Zdeborová, PNAS 104, 10318 (2007).
  • Mézard et al. (2005) M. Mézard, T. Mora, and R. Zecchina, Phys. Rev. Lett. 94, 197205 (2005).
  • Zhou and Ma (2009) H. Zhou and H. Ma, Phys. Rev. E 80, 066108 (2009), URL
  • Mézard and Zecchina (2002) M. Mézard and R. Zecchina, Phys. Rev. E 66, 056126 (2002).
  • Mézard et al. (2002) M. Mézard, G. Parisi, and R. Zecchina, Science 297, 812 (2002).
  • Montanari et al. (2004) A. Montanari, G. Parisi, and F. Ricci-Tersenghi, J. Phys. A 37, 2073 (2004).
  • Ardelius and Aurell (2006) J. Ardelius and E. Aurell, Phys. Rev. E 74, 037702 (2006).
  • Seitz et al. (2005) S. Seitz, M. Alava, and P. Orponen, JSTAT p. P06006 (2005).
  • Alava et al. (2008) M. Alava, J. Ardelius, E. Aurell, P. Kaski, S. Krishnamurthy, P. Orponen, and S. Seitz, PNAS 105, 15253 (2008).
  • Ardelius and Zdeborová (2008) J. Ardelius and L. Zdeborová, Phys. Rev. E 78, 040101(R) (2008).
  • Ardelius et al. (2007) J. Ardelius, E. Aurell, and S. Krishnamurthy (2007), preprint, eprint arXiv:cond-mat/0702672.
  • Papadimitriou (1991) C. H. Papadimitriou, in Proceedings of the 32nd Annual IEEE Symposium on Foundations of Computer Science, FOCS’91 (San Juan, Puerto Rico, October 1-4, 1991) (IEEE Computer Society Press, Los Alamitos-Washington-Brussels-Tokyo, 1991), pp. 163–169.
  • Selman et al. (1996) B. Selman, H. Kautz, and B. Cohen, in Cliques, Coloring, and Satisfiability, edited by D. S. Johnson and M. A. Trick (American Mathematical Society, 1996), DIMAC Series in Discrete Mathematics and Theoretical Computer Science, pp. 521–532.
  • Kautz et al. (1997) H. Kautz, B. Selman, and Y. Jiang, in Satisfiability Problem: Theory and Applications, edited by D. Du, J. Gu, and P. M. Pardalos (American Mathematical Society, 1997), vol. 35 of Dimacs Series in Discrete Mathematics and Theoretical Computer Science, pp. 573–586.
  • Davis et al. (1962) M. Davis, G. Logemann, and D. Loveland, Communications of the ACM 5, 394 (1962), URL
  • Gomes et al. (2008) C. P. Gomes, H. Kautz, A. Sabharwal, and B. Selman, in Handbook of Knowledge Representation (Elsevier, 2008), vol. 3 of Foundations of Artificial Intelligence, pp. 89–134.
  • Bayardo and Pehoushek (2000) R. J. Bayardo and J. D. Pehoushek, in Proceedings of the 17th AAAI (AAAI Press, Menlo Park, 2000), pp. 157–162.
  • Cocco and Monasson. (2001) S. Cocco and R. Monasson., Phys. Rev. Lett. 86, 1654 (2001).
  • Monasson and Zecchina (1996) R. Monasson and R. Zecchina, Phys. Rev. Lett. 76, 3881 (1996).
  • Biroli et al. (2000) G. Biroli, R. Monasson, and M. Weigt, Eur. Phys. J. B 14, 551 (2000).
  • Hartmann (2000a) A. K. Hartmann, J. Phys. A 33, 657 (2000a).
  • Hartmann (2000b) A. K. Hartmann, Eur. Phys. J. B 13, 539 (2000b).
  • Hed et al. (2001) G. Hed, A. K. Hartmann, D. Stauffer, and E. Domany, Phys. Rev. Lett. 86, 3148 (2001).
  • Geyer (1991) C. Geyer, in 23rd Symposium on the Interface between Computing Science and Statistics (Interface Foundation North America, Fairfax, 1991), p. 156.
  • Newman and Barkema (1999) M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics (Oxford University Press, 1999).
  • Hukushima and Nemoto (1996) K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
  • Moreno et al. (2003) J. J. Moreno, H. G. Katzgraber, and A. K. Hartmann, Int. J. Mod. Phys. C 14, 285 (2003).
  • Hartmann and Ricci-Tersenghi (2002) A. K. Hartmann and F. Ricci-Tersenghi, Phys. Rev. B 66, 224419 (2002).
  • Hartmann (2009) A. K. Hartmann, Practical Guide to Computer Simulations (World Scientific, 2009).
  • Ward Jr. (1963) J. H. Ward Jr., J. Amer. Stat. Assoc. 58, 236 (1963).
  • Jain and Dubes (1988) A. Jain and R. Dubes, Algorithms for Clustering Data, Advanced Reference Series (Prentice Hall, 1988).
  • Higgs (1996) P. G. Higgs, Phys. Rev. Lett. 76, 704 (1996).
  • Mézard et al. (1984) M. Mézard, G. Parisi, N. Sourlas, G. Toulouse, and M. Virasoro, Phys. Rev. Lett. 52, 1156 (1984).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description