On the Empirical Time Complexity of Random 3-SAT at the Phase Transition

1 Introduction

The propositional satisfiability problem (SAT) was the first problem proven to be NP-complete [Cook, 1971], and no known algorithm solves the problem efficiently in the sense of better-than-exponential scaling of running time with instance size in the worst case. SAT also has a broad range of practical applications, and despite its discouraging worst-case time complexity, many SAT solvers have been constructed that perform well on theoretically and practically interesting sets of benchmark instances. The performance of SAT solvers are regularly assessed in SAT competitions [SatCompetition.org, 2014] , which feature benchmark instances from various applications, instances that are hand-crafted to challenge solvers, and randomly generated instances.

Following seminal work by [Cheeseman et al., 1991] and [Mitchell et al., 1992], randomly generated 3-SAT instances at the so-called solubility phase transition, where 50% of the instances generated at a given size are satisfiable, have been widely studied; to this day, they are regarded as a model for the difficulty of the NP-complete 3-SAT problem and represent one of the most prominent benchmarks for SAT solvers. It is conjectured, yet unproven, that the location of the solubility phase transition for uniform random 3-SAT, in terms of the ratio of the number of clauses and variables, m/n, converges towards a limiting value as n approaches infinity. [Crawford and Auton, 1996] studied the location of the phase transition for random 3-SAT as a function of n and provided a widely used formula, based on empirical data for n = 20 . . . 300.

Here, we consider the question how the performance of prominent, high-performance SAT solvers on uniform random 3-SAT at the solubility phase transition scales with n. Prompted in part by earlier, inconclusive results suggesting polynomial scaling of incomplete SAT solvers based on stochastic local search (SLS) [Gent and Walsh, 1993; Parkes and Walser, 1996; Gent et al., 1997], we characterise the empirical scaling of the performance of prominent incomplete, SLS-based and complete, DPLL-based SAT solvers. We are interested in determining whether there are qualitative differences in the scaling behaviour of SLS-based and DPLL-based solvers, in the scaling of DPLL-based solvers on satisfiable and unsatisfiable instances, and in the scaling of different solvers within the two major groups (SLS-based and DPLL-based). Our investigation leverages and extends a recent approach for empirical scaling analysis by [Hoos and Stützle, 2014] that uses resampling techniques to statistically assess the degree to which observed algorithm behaviour is consistent with given scaling models (see Sec. 3).

To avoid 'falling off' the phase transition, which could bias our scaling results towards easier instances, we re-examine the model in [Crawford and Auton, 1996] for the location of the phase transition point, and derive a new model that agrees better with observed data for n > 300 and with recent results from statistical physics (see Sec. 4). Using sets of random 3-SAT instances generated according to this model, our main findings are as follows (see Sec. 5):

• The median running times of the three prominent SLS-

based solvers we studied, WalkSAT/SKC [Selman et al., 1994], BalancedZ [Li et al., 2014] and probSAT [Balint and Schöning, 2014; 2012], scale polynomially (with an exponent of about 3), and exponential scaling models are rejected with 95% confidence. Furthermore, we found no evidence that higher percentiles of the distributions of running times over sets of instances for fixed n may scale exponentially.

• The median running times of the three DPLL-based solvers we considered, kcnfs [Dequen and Dubois, 2004], march_hi [Heule and van Maaren, 2009] and march_br [Heule, 2013], exhibit exponential scaling (with a base of about 1.03), and polynomial models are rejected with 95% confidence.
• For all three DPLL-based solvers, the median running times when solving only satisfiable and only unsatisfiable instances, respectively, clearly exhibit exponential running time scaling, and the respective scaling models differ mainly by a constant factor.
• While the scaling models for the SLS-based solvers are very similar to each other, the two march-variants scale significantly better than kcnfs.

As we will discuss in the final part of this paper (Sec. 6), these results shed new light on random 3-SAT at the phase transition and on the fundamental differences in the behaviour of SLS- and DPLL-based solvers; they also invite thorough empirical scaling studies for other problems and solvers.

2 Background and RelatedWork

SAT is one of the most intensely studied problems in the computer science literature and beyond, and 3-SAT is arguably the most prominent NP-complete decision problem [Cook, 1971]. Interest in phase transition phenomena in combinatorial problems and in uniform random 3-SAT specifically rose sharply when [Cheeseman et al., 1991] demonstrated that the hardest instances are found around a critical value of an order parameter, where a transition from predominantly soluble to mostly insoluble instances occurs. Uniform random k-SAT instances are generated by constructing uniformly and independently at random m clauses, each of which is obtained by sampling, again uniformly and independently at random, 3 of n variables, and negating each of them with probability 1/2 [Mitchell et al., 1992] (duplicate clauses are eliminated); the order parameter is the clauses/variables ratio, m/n. It is believed, but not yet proven, that for k ≥ 3, the location of the phase transition point of uniform random k-SAT converges to a fixed threshold value as n approaches infinity (it is known, however, that the phase transition is provably sharp for n → ∞ [Friedgut and Bourgain, 1999]). Assuming threshold values exist for small k, an accurate theoretical upper bound was proven by [Franco and Paull, 1983] and was later improved to 2k log 2 - (1 + log 2) /2 + ok (1) by [Kirousis et al., 1998]. [Achlioptas and Peres, 2004] achieved a major breakthrough in proving a lower bound, which was recently improved to 2k log 2-(1 + log 2) /2-ok (1) [CojaOghlan, 2014].

In a prominent study on the empirical difficulty of SAT instances, [Mitchell et al., 1992] demonstrated that instances drawn from the phase transition region of uniform random 3SAT instances tend to be the most difficult for a simple DPLL solver. Similar results were shown by [Yokoo, 1997] for an SLS-based solver on the satisfiable phase of uniform random 3-SAT. [Crawford and Auton, 1996] studied the phase transition region of uniform random 3-SAT empirically, developed a model for the location of the phase transition point and presented additional evidence for the difficulty of the SAT instances found there.

[Gent and Walsh, 1993] studied the empirical behaviour of GSAT, one of the earliest and most prominent SLS-based SAT solvers [Selman et al., 1992], and its two variants, DSAT and HSAT. They noted that the scaling of the average number of variable flips required by these solvers for solving phase transition random 3-SAT instances was consistent with lessthan-linear exponential scaling with n and did not rule out a polynomial scaling model with degree ≈ 3. Later, [Parkes and Walser, 1996] presented empirical evidence that the scaling of the average number of flips required by a more recent, prominent SLS-based SAT solver, WalkSAT/SKC [Selman et al., 1994], on the same class of instances might scale either as a power function with a slowly growing exponent, or as a sub-exponential function. Furthermore, [Gent et al., 1997] found that the 90th percentile of the number of flips of GSAT for random 3-SAT appears to grow no faster than n4 . However, in all cases, performance was measured in variable flips (which become more expensive as n increases) and the results are based on limited curve fitting only, with a vaguely defined notion of 'good fit'. [Coarfa et al., 2003] used simple curve fitting to study the empirical scaling of median running time for three complete solvers on random 3-SAT and observed exponential scaling above certain solver-dependent density thresholds. In contrast, in this work, we use a more advanced and statistically sound approach for assessing scaling models. Unlike these earlier studies, we also challenge our scaling models by assessing their predictions on larger instances sizes than those used for fitting the models.

A significant advance in the methodology for studying the empirical scaling of algorithm performance with input size was achieved by [Hoos, 2009]. Their method uses standard numerical optimisation approaches to automatically fit scaling models, which are then challenged by extrapolation to larger input sizes. Most importantly, it uses a resampling approach [Efron and Tibshirani, 1993] to assess the models and their predictions in a statistically meaningful way. [Hoos and Stützle, 2014] used this approach to characterise the scaling behaviour of Concorde, a state-of-the-art complete algorithm for the travelling salesperson problem (TSP), demonstrating that the scaling of Concorde's performance on 2D Euclidean TSP instances with n cities agrees well with a model of the form a · b√n . In this work, we apply the same empirical methodology to SAT, and extends it in two useful ways.

3 Experimental Setup andMethodology

For our study, we selected three SLS-based solvers, WalkSAT/SKC [Selman et al., 1994] (UBCSAT version 1.2

[Tompkins and Hoos, 2005]), BalancedZ [Li et al., 2014] (V2014.05.07) and probSAT [Balint and Schöning, 2014] (version sc14), and three DPLL-based solvers, kcnfs [Dequen and Dubois, 2004] (V2004), march_hi [Heule and van Maaren, 2009] and march_br [Heule, 2013] (version sat+unsat). We chose these, because WalkSAT/SKC and kcnfs are two classical solvers that are widely known in the community, and the others showed excellent performance in recent SAT competitions.

All sets of uniform random 3-SAT instances used in our study were obtained using the generator which produced SAT competition instances [Balint et al., 2012]. To separate satisfiable from unsatisfiable instances with n ≤ 500, we used CSHCrandMC [Malitsky et al., 2013], the winner of Random SAT+UNSAT Track in the 2013 SAT Competition. For larger instances, for which the use of complete solvers is impractical, we performed long runs (43200 CPU sec) of BalancedZ and treated those instances not solved in those runs as unsatisfiable. Since BalancedZ never failed to solve any of our instances known to be satisfiable, and because the cut-off time we used for these long runs was at least 20 times of the maximum running time observed on any satisfiable instance of the same size, we have high confidence that we did not miss any satisfiable instances. We note that, even if we had misclassified some satisfiable instances, the quantile performance measures would not be much affected (if at all).

For collecting running time data for our solvers, we used the Compute Canada / Westgrid cluster orcinus (DDR), each equipped with two Intel Xeon E5450 quad-core CPUs at 3.0 GHz with 16GB of RAM running 64-bit Red Hat Enterprise Linux Server 5.3. To avoid measurement noise caused by memory bottlenecks, we used only one core per CPU to run our solvers and imposed a memory limit of 1GB per run.

After observing floor effects for WalkSAT/SKC due to small CPU times required for solving small instances, we calculated running times based on the number of variable flips performed, and on estimates of CPU time per variable flip for each problem size obtained from long runs on unsatisfiable instances. The CPU time estimates thus obtained closely agreed with measured CPU times for short and long runs, but are considerably more accurate. For all other solvers, which were less affected by such floor effects, we used runsolver [Roussel, 2011] to record CPU times (and to enforce CPU time limits).

For our scaling analysis, we considered two parametric models:

• Exp [a, b] (n) := a · bn (2-parameter exponential);
• Poly [a, b] (n) := a · nb (2-parameter polynomial).

Our approach could be easily extended to other scaling models, but, as we will show in the following, these two models jointly characterise the scaling observed in all our experiments, and thus we saw no need to consider different or more complex models. For fitting parametric scaling models to observed data, we used the non-linear least-squares LevenbergMarquardt algorithm.

Models were fitted to performance observations in the form of quantiles, chiefly the median of the distributions of running times over sets of instances for given n. Compared to the mean, the median has two advantages: it is statistically more stable and immune to the presence of a certain amount of timed-out runs. Considering the stochastic nature of the SLS-based solvers, we performed 5 independent runs per instance and used the median over those 5 running times as the running time for the respective instance. To assess the fit of a given scaling model to observed data, we used root-meansquare error (RMSE).

Closely following [Hoos, 2009; Hoos and Stützle, 2014], we computed 95% bootstrap confidence intervals for the performance predictions obtained from our scaling models, based on 1000 bootstrap samples per instance set and 1000 automatically fitted variants of each scaling model. We extended their approach in two ways. Firstly, noting that observed running time statistics are also based on measurements on sets of instances, we calculate bootstrap percentile intervals for those, in addition to the point estimates used in the original approach. This way, we capture dependency of these statistics on the underlying sample of instances. Secondly, we determine to which extent a solver A1 shows scaling behaviour different from another solver A2, by comparing the observed running time statistics of A1 to the bootstrap confidence intervals obtained from the scaling model of A2. If the latter do not contain the former, we can reject the hypothesis that the performance of A1 is consistent with the scaling model of A2. Both extensions can obviously be used in the scaling analysis of solvers for problems other than SAT.

4 Location of Phase Transition

[Crawford and Auton, 1996] modelled the location of the phase transition point based on extensive experiments on uniform random 3-SAT instances with up to 300 variables as:

where n is the number of variables and mc the critical number of clauses, at which about 50% satisfiable instances are obtained. The parametric form of this model was derived by [Kirkpatrick and Selman, 1994], using finite-size scaling. While this model does provide a good fit for n ≤ 300, in preliminary experiments, we found that for n > 300, it increasingly underestimates mc. Furthermore, the model in Eq. 1 predicts that, as n grows, mc/n approaches 4.258, which is in contradiction with a more recent result by [Mertens et al., 2006], who used the heuristic one-step replica symmetry breaking cavity method from statistical physics to estimate the value of mc/n as 4.26675 ± 0.00015.

Because in the empirical scaling study that follows, we wanted to be sure to not drift off the phase transition point (which could bias the scaling models, especially, as the phase transition, in terms of mc/n, is known to become increasingly steep as n grows), we decided to revisit and improve the Crawford & Auton's model.

We first chose several m/n ratios for each n ∈ {300, 400, . . . , 1400}, around the predictions from Eq. 1, loosely adjusted based on results from preliminary experiments. Next, we generated sets of 600 uniform random 3SAT instances for each m/n ratio. We separated out the satisfiable instances using a hybrid complete solver, CSHCrandMC [Malitsky et al., 2013], with a cutoff of 3600 CPU

seconds per run. Up to n = 500, we solved all instances within that cutoff. Beyond, it would have been impractical to run any complete solver sufficiently long to prove unsatisfiability, and we therefore heuristically assumed that the instances not solved by CSHCrandMC are unsatisfiable. As mentioned earlier, we later used much longer runs of BalancedZ to challenge these putative unsatisfiable instances further, and among the thousands of instances for which this was done, only one was found to be satisfiable. Nevertheless, in what follows, we are well aware of the fact that we may underestimate the fraction of satisfiable instances in our sets.

We note that, even at large numbers of instances sampled at the correct (unknown) value of mc/n, we should expect to get sets with a fraction of satisfiable instances varying according to a binomial distribution. Based on this observation, we determined 95% confidence intervals for the fraction of satisfiable instances observed. We then rejected the m/n values for which we empirically observed fractions of satisfiable instances outside of the respective 95% confidence interval as valid estimates of mc/n. In this way, we obtained bounds on the location of the phase transition point, mc/n. In many cases, the confidence intervals for our initial sets of 600 instances were too wide to yield bounds, and we subsequently increased the number of instances in those sets until for every n, we obtained an upper and lower bound on mc/n. Those bounds and the respective set sizes are shown in Table 1.

These results provide further evidence for our earlier observation that for larger n, Crawford & Auton's model becomes inaccurate, as the respective predictions are below the lower bounds from our analysis. We note that our lower bounds are valid, as they could only increase if some of our putatively unsatisfiable instances were in fact satisfiable. For the same reason, the upper bounds may be inaccurate.

Interestingly, and notably different from Crawford & Auton's model, our results suggest that mc/n as a function of n is not monotonic, but first drops, before slowly increasing again, possibly towards a threshold value. This observation, in combination with the limiting value found by [Mertens et al., 2006] led us to choose a model of the form

Finally, taking the data points from [Crawford and Auton, 1996] up to n = 300 (which we believe to be of high quality) and the mid-points between our bounds from Table 1, we fitted this 4-parameter model, resulting in:

Figure 1 shows the predictions obtained from this model, along with the bounds and mid-points between them from Table 1. This model was subsequently used to generate the instance sets used in the scaling analysis described in the following. While it is possible that our model is biased by the fact that we may have missed satisfiable instances for larger n, it provides a better basis than previously available for generating instances at or very near the solubility phase transition of uniform random 3-SAT.

5 Empirical Scaling of Solver Performance

We first fitted our two parametric scaling models to the median running times of the six solvers we considered, as described in Section 3. For each SLS-based solver, both models were fitted using median running times for n = 200, 250, . . . 500 (support) and later challenged with median running times for n = 600, 700, . . . 1000. For each DPLLbased solver, we used support data for n = 200, 250, . . . 400 and challenge data for n = 450, . . . 550, as for even larger n we could no longer complete sufficiently many runs to estimate even median running times. This resulted in the models, shown along with RMSEs on support and challenge data, shown in Table 2.

We note that the RMSEs on support data, i.e., the data used for fitting the models, often, but not always provide a good indication of their predictive power. Based on the latter, in the form of RMSEs on challenge data, we see clear indications that the median running times of all three SLS-based solvers are overall more consistent with our polynomial scaling model, while those of the DPLL-based solvers are more consistent with our exponential scaling model - even when only considering performance on satisfiable instances. We also studied Lingeling [Biere, 2014] (version ayv) and found that its median running time for n = 400 is more than 1000 times larger than that of kcnfs, and that this performance gap increases with n.

1 This type of model is consistent with current theoretical thinking that, for small n, second-order terms may cause non-monotonic behaviour of the function mc(n). (Dimitris Achlioptas, personal communication.)

But how much confidence should we have in these models? Are the RMSEs small enough that we should accept them? 2 To answer this question, we assessed the fitted models using the bootstrap approach outlined in Section 3. The results of this analysis, shown in Table 3, clearly show that observed median running times for the SLS-based solvers are consistent with our polynomial scaling model and inconsistent with the exponential model (as illustrated in Figure 2 for WalkSAT/SKC), while the opposite holds for the DPLLbased solvers. This is especially striking for the larger challenge instances sizes. Limited experiments for even larger instances sizes (up to n = 2000) produced further evidence consistent with the polynomial scaling models for all three SLS-based solvers (data not shown). Although the running time of the SLS-based solvers are still quite moderate even at these instances sizes, the much longer runs required to filter out satisfiable instances with high confidence make it difficult to further increase n.

Next, we used the same bootstrap confidence intervals to investigate the significance of differences observed between the scaling models for different solvers. Using this approach, we cannot reject the hypothesis that the differences reflected in the constants for the polynomial models of our three SLSbased solvers seen in Table 2 are insignificant. Furthermore,

2 We note that RMSE values are sensitive to the absolute magnitude of the differences in the given data. It is therefore difficult to judge different models solely based on absolute RMSE.

n

Figure 2: Fitted models of the median running times of WalkSAT/SKC. Both models are fitted with the median running times of walksat solving the SAT instances from the set of 1200 random instances of 200, 250, ..., 500 variables, and are challenged by median running times of 600, 700, ..., 1000 variables.

we observed that the median running times of CCASat [Cai and Su, 2012] (V2014), another prominent SLS-based solver, are consistent with the same polynomial models. Likewise, the bootstrap confidence intervals for march_hi and march_br largely overlap with each other, and we cannot detect sufficient evidence for statistically significant differences in scaling behaviour. On the other hand, the observed running times for kcnfs are inconsistent with the scaling models for march_hi and march_br, confirming that the performance of the march solvers indeed scales significantly better than that of kcnfs. Examining these results in detail, we believe that by substantially increasing the number of instances for each value of n, the bootstrap confidence intervals can likely be narrowed to also demonstrate significant scaling differences between all pairs of SLS-based solvers.

Using the same approach, we investigated the significance of the differences in scaling for each DPLL-based solver when applied to satisfiable and unsatisfiable instances only, respectively. Intuitively, we would expect unsatisfiable instances to be harder, and the differences between the respective scaling models are significant, in that the observed running times for solving unsatisfiable instances are generally inconsistent with the scaling model for the same solver on only satisfiable instances (as illustrated in Figure 3 for march_hi). To further investigate whether the performance difference can be attributed to the constant factor a in the exponential models, we fitted a model of the form a · bnsat of the median running times for solving unsatisfiable instances, where a is a free parameter and bsat is the value of parameter b from the scaling model for satisfiable instances. For all three DPLLbased solvers, the observed running times for solving unsatisfiable instances is consistent with these constrained models, suggesting that, indeed, the scaling on satisfiable and unsatisfiable instances differs only by a constant factor.

Finally, we investigated the question whether our observations regarding the qualitative differences in the scaling of median running times of SLS- and DPLL-based solvers also hold when considering higher quantiles of the distribution of running times. For the DPLL-based solvers, we found no evidence for a substantial change in the ratios between the median and higher quantiles as n increases; thus, there is no reason to assume any difference in the scaling models of a higher quantile and the median other than a constant factor. For the SLS-based solvers, we noted weak evidence that the ratio between the higher quantiles and the median increases with n, but the observed scaling of the 0.9-quantile is still found to be consistent with a polynomial model and inconistent with an exponential model (data not shown).

6 Conclusions

In this work, we presented an empirical analysis of the scaling behaviour of several prominent SAT solvers on phase transition uniform random 3-SAT instances. We were surprised to find solid support for low-degree polynomial scaling of the performance of all SLS-based SAT solvers we studied, which stands in stark contrast to the exponential performance scaling of the DPLL-based solvers we considered, even when these were run on satisfiable instances only. As expected, we did find evidence for significant differences in the performance scaling models of DPLL-based solvers on satisfiable and unsatisfiable instances, but these differences could be attributed to a constant factor, suggesting that they cannot be leveraged to obtain reasonably accurate guesses on the unsatisfiability of instances based on long runs of those solvers. However, the qualitative differences between the scaling of SLS- and DPLL-based solvers can be exploited, in that for growing n, it appears to be possible to separate satisfiable from unsatisfiable instances with high and further increasing accuracy based on long runs of SLS-based or hybrid solvers. Furthermore, the polynomial scaling of even high quantiles of the distribution of running times for SLS-based solvers across instance sets suggests that random 3-SAT instances from the solubility phase transition are likely not capturing the difficulty we expect to encounter in the worst case when solving an NP-hard problem. We note that, considering their sizes, these instances are still very hard, compared to almost all types of structured SAT instances, and may therefore still be useful as a benchmark.

To the best of our knowledge, ours is the first study that uses a statistically sound way to assess the scaling of SAT solver performance with instance size, and to discriminate between different scaling models. In the context of our study, we have extended the methodology introduced by [Hoos, 2009; Hoos and Stützle, 2014] in two ways. The empirical scaling analysis we have performed here can easily be applied to other SAT solvers, other distributions of SAT instances (as long as reasonably large sets of instances for each n can be obtained), and to other problems. We believe that doing so can and will produce interesting and useful results that can inspire the design of algorithms and benchmark instance generators as well as, hopefully, theoretical work.

Acknowledgements

We thank the annoymous reviewers for their thoughtful comments and suggestions. We also thank Dimitris Achlioptas, Paul Beame, Henry Kautz, Donald Knuth, Bart Selman and Moshe Vardi for valuable advice regarding previous work on the 3-SAT phase transition and empirical performance of SAT algorithms. This work was supported by Compute Canada, by the Institute for Computing, Information and Cognitive Systems (ICICS) at UBC, and by NSERC, through a Discovery Grant to H.H.