Estimating the Density of States of Boolean Satisfiability Problems on Classical and Quantum Computing Platforms
aa r X i v : . [ c s . D M ] O c t Estimating the Density of States of Boolean SatisfiabilityProblems on Classical and Quantum Computing Platforms
Tuhin Sahai Anurag Mishra
United Technologies Research Center,Berkeley, CA 94705.
José Miguel Pasini
United Technologies Research Center,East Hartford, CT 06118.
Susmit Jha
Computer Science LaboratorySRI International,Menlo Park, CA 94025.
Abstract
Given a Boolean formula φ ( x ) in conjunctive normal form (CNF), the density of statescounts the number of variable assignments that violate exactly e clauses, for all values of e . Thus, the density of states is a histogram of the number of unsatisfied clauses over allpossible assignments. This computation generalizes both maximum-satisfiability (MAX-SAT) and model counting problems and not only provides insight into the entire solutionspace, but also yields a measure for the hardness of the problem instance. Consequently,in real-world scenarios, this problem is typically infeasible even when using state-of-the-artalgorithms. While finding an exact answer to this problem is a computationally intensivetask, we propose a novel approach for estimating density of states based on the concentra-tion of measure inequalities. The methodology results in a quadratic unconstrained binaryoptimization (QUBO), which is particularly amenable to quantum annealing-based solu-tions. We present the overall approach and compare results from the D-Wave quantumannealer against the best-known classical algorithms such as the Hamze-de Freitas-Selby(HFS) algorithm and satisfiability modulo theory (SMT) solvers.
1. Introduction
The density of states (DOS) for a given Boolean formula not only provides insight into thecomplete solution space but also serves as an accurate measure of the difficulty or hard-ness of the problem instance. The ability to compute DOS for optimization and feasibilityproblems has critical applications in system-requirements engineering of complex aerospaceproducts. It provides a metric for requirements engineers to compare constraints, prescribedrequirements (Ferrante, Scholte, Pinello, Ferrari, Mangeruca, & Liu, 2016; Johnson et al.,1998; Costello & Liu, 1995), and requirements decompositions (Kirkman, 1998). This com-putation is particularly germane to the design and optimization of complex aerospace sys-tems (Sommerville, 2005). Current classical methods for computing density of states (Wang& Landau, 2001; Ermon, Gomes, & Selman, 2010, 2011) have limited scalability. While thefocus of this paper is on Boolean formulae, we note that constrained programming and fea-1ibility problems can be easily mapped to equivalent Boolean satisfiability instances (Walsh,2000; Tamura, Taga, Kitagawa, & Banbara, 2009).The DOS problem in the standard k -satisfiability ( k -SAT) setting can be elucidated asfollows: instead of deciding whether a given logical formula is satisfiable or not, one aimsto compute the entire histogram of the number of clauses satisfied over all possible variableassignments . Note the DOS, for a given instance of an optimization or decision problem,captures its hardness (distributions with a low footprint for all satisfied clauses are harderto compute or satisfy). The DOS histogram sheds light on the fundamental nature of thefeasible solution set and difficulty of solving the overall optimization problem. As statedearlier, these problems frequently arise when constructing complex systems for aerospaceand defense applications (Leveson & Weiss, 2009).The lack of methodical approaches that enable the comparison of competing safety-critical system requirements, while optimizing performance, stymie the development of next-generation complex systems. Note there are often multiple paths to decompose the overallsystem safety requirements down to subsystems requirements. Some of these decomposi-tions may lead to costly design and redesign cycles to achieve desired levels of performance.Decompositions that have a higher DOS in the satisfiable range result in greater freedom tooptimize performance and, consequently, result in quicker design cycles and fewer redesigns.The ability to quickly estimate the DOS of satisfiability problems will enable the specifi-cation engineer to ensure the prescribed requirements are satisfiable, internally consistent,and amenable to design space exploration very early in the design requirement step.In this work, we aim to construct novel approaches for rapidly computing the DOS fora SAT problem (Boolean formula) (Biere, Heule, & van Maaren, 2009). Our approximateapproach to estimate DOS of SAT instances exploits the concentration of measure inequal-ities (Boucheron, Lugosi, & Massart, 2013). These inequalities provide bounds on the tailsof the distributions of random functions and have been used to construct the theory ofgeneralization in machine learning (Abu-Mostafa, Magdon-Ismail, & Lin, 2012), computeoptimal bounds on uncertainty (Owhadi, Scovel, Sullivan, McKerns, & Ortiz, 2013), certifysystems (Leyendecker, Lucas, Owhadi, & Ortiz, 2010), compute bias of statistical estima-tors (Gourgoulias, Katsoulakis, Rey-Bellet, & Wang, 2017), and derive results in randommatrix theory (Tao, 2012).In this paper, we make the following contributions:1. We introduce a novel approach to estimate the density of states for SAT problems byusing concentration of measures (McDiarmid’s inequality), and bound the deviationof the number of unsatisfied clauses (energy) from the expected (mean) number ofunsatisfied clauses for uniformly distributed assignments.2. The deviation of the energy function from its expected value depends on its diameter(function variability), which can be computed by solving an optimization (maximiza-tion) problem (Owhadi et al., 2013). We show this maximum deviation computationcan be posed in the form of a quadratic unconstrained binary optimization (QUBO)that is particularly amenable to quantum annealers and results in tight bounds onthe DOS histogram. 2. We demonstrate our approach on classical platforms by computing the diameter andassociated concentration of measure bounds using Selby’s implementation (Selby, 2013,2014) of the Hamze-de Freitas- Selby (HFS) algorithm (Hamze & de Freitas, 2004).4. We use satisfiability modulo theory (SMT) solvers (Bjørner, Phan, & Fleckenstein,2015; Sebastiani & Trentin, 2015; Dutertre, 2014) to solve the QUBO formulationas an alternative approach on classical platforms. The solutions from SMT solversprovide tighter estimates but require significantly higher computational effort and donot scale.5. We then compare the classical results to the computations on the D-Wave quan-tum annealer, a commercially available noisy intermediate-scale quantum (NISQ) de-vice (Preskill, 2018). We find the D-Wave machine provides higher-quality solutionswhen compared to the HFS algorithm, and scales better than SMT solvers. We fur-ther note the search for useful problems that are appropriate for present day NISQdevices is a very active area of research within quantum computing (Preskill, 2018).We propose the DOS computation task as a potential test problem that can be usedto benchmark current- and next-generation quantum annealers against their classicalcounterparts .The rest of the paper is organized as follows. In section 2, we introduce the SAT problemand its associated density of states. We then discuss the existing state-of-the-art methodsfor computing DOS and highlight their limitations. Section 3 introduces concentration ofmeasure inequalities in the context of DOS computations for the SAT problem. Section 4presents our novel formulation of the DOS problem as a QUBO for estimating the associatedenergy histogram. In section 5, we present results comparing state-of-the-art algorithms forcomputing density of states with our proposed concentration of measures approach on classi-cal and quantum platforms. For the concentration of measures approach, we further presentcomparisons between the performance of the D-Wave machine with the HFS algorithm andthe Z3 SMT solver. We conclude in section 6 by summarizing our results and presentingkey challenges.
2. Background
Consider a k -SAT formula φ ( x ) of N binary variables and m clauses, φ : { , } N → { , } ,written in the conjunctive normal form (CNF) (Biere et al., 2009) as follows, φ ( x ) = m ^ i =1 C i = m ^ i =1 ( x i ∨ x i ∨ . . . ∨ x i k ) , (1)where x i l is the l th literal in clause C i . A SAT formula is said to be satisfiable if there existsan assignment for the binary variables x such that φ ( x ) = 1 (true). It is well known that thesatisfiability problem is NP-complete (Cook, 1971). A critical parameter associated withthe satisfiability problem is the clause density α = m/N (Biere et al., 2009). In particular,the probability that a random k -SAT instance is satisfiable undergoes a phase transition asa function of α ( N → ∞ ) (Xu & Li, 2000; Biere et al., 2009). The MAX-SAT problem (andthe corresponding weighted version) (Krentel, 1988; Chieu & Lee, 2009) requires one to find3hat assignment (or assignments) that maximize the number (or the cumulative weights) ofsatisfied clauses. Consider a SAT formula φ , then every assignment x can be mapped to an“energy” Φ( x ) such that, Φ( x ) = m X i =1 C i , (2)where C i = 1, if the i -th clause evaluates to true. In other words, the goal under the MAX-SAT problem is to find the assignment for x such that the number of satisfied clauses (orenergy) is maximized. Using De Morgan’s laws, one can easily show that,Φ( x ) = m − m X i =1 k Y l =1 f ( x i l ) , (3)where f ( x i l ) = ( x i l , if x i l is negated in the clause(1 − x i ) , otherwise . (4)Using the above formula, it is easy to see that if each literal x i l were random with equalprobability for values { , } , then the expected number of satisfied clauses is, E [Φ( x )] = m (2 k − k . (5)Thus, even though the satisfiability is NP-complete, a random assignment is expected tosatisfy a large fraction of the clauses. Note that for the 3-SAT, the above formula reducesto, E [Φ( x )] = 7 m/ . (6)One can use this expected (mean) value of the number of satisfied clauses to estimate theDOS using concentration of measure inequalities (see section 3). For SAT instances thatarise from specific application domains (thus, not random), one can estimate the expectednumber of satisfied clauses by sampling over the independent variables in the Booleanformula.The DOS d ( e ) of a SAT formula φ is equal to the number of assignments x for whichΦ( x ) = e . In other words, it is the histogram of the number assignments as a function of e satisfied clauses. Note the value of the number of satisfied clauses e lies between 0 and m where m is the total number of clauses in the SAT formula. Following the terminologyfrom the physics community, we will also call this the energy of SAT formula. Since thetotal number of possible assignments is 2 N , one can define the normalized density of statesas follows, p ( e ) = d ( e )2 N . (7)The normalized DOS acts as a discrete probability distribution. Note it is not necessarythat all energies e have a valid assignment. For example, if the SAT formula cannot besatisfied then p ( m ) = 0. As explained previously, for a random 3-SAT, the mean of p ( e ) vs e is at 7 m/ σ is inverselyproportional to its DOS. Thus, the sampling approach effectively converges to a flat-visithistogram (captured by a flatness parameter). In (Ermon et al., 2010, 2011), the authors testthe algorithm on multiple benchmark examples. We use this method to compute the DOSfor a series of random SAT instances and compare the resulting histograms to estimatedDOS using our concentration of measure-based approach (as outlined in section 3).Fig. 1 describes a simple example of the DOS problem for a Boolean satisfiability prob-lem with N = 100 and α = 4 .
0, and shows an example output of the MCMC-FlatSATalgorithm. In our experiments, we found that for k -SAT instances close to the phase tran-sition (Monasson, Zecchina, Kirkpatrick, Selman, & Troyansky, 1999), the mixing times ofthe Markov chain (Levin & Peres, 2017) increase significantly.0 20 40 60 80 100 1200204060 Number of unsat clauses D e n s i t y o f s t a t e s ( L og2 S c a l e ) Figure 1: This example Boolean satisfiability problem has N = 100 variables, and the clausedensity ratio α is 4 .
0. The x-axis is the number of UNSAT clauses and the y-axisprovides a numerical value for the number of assignments with the correspondingnumber of UNSAT clauses. The DOS over all possible variable assignments arecaptured in the histogram.
3. SAT and Concentration of measure inequalities
In this section, we describe our approach based on a concentration of measure inequalityfor estimating DOS of satisfiability problems, and its relationship to the quadratic uncon-strained binary optimization (QUBO) problem.The concentration of measure phenomena bounds the deviation of functions of randomvariables around their mean (Boucheron et al., 2013). There are a host of inequalitiesassociated with various situations, see (Boucheron et al., 2013; Tao, 2010) for more details.5or our setting, we use McDiarmid’s inequality (McDiarmid, 1989) summarized in thetheorem below.
Theorem. (McDiarmid’s inequality) Let x , x , . . . , x N be random variables taking valuesin the range R , R , . . . , R N and let F : R × R × . . . × R N → R be a function with theproperty that if one freezes all but the i -th variable, then F ( x , x , . . . , x N ) fluctuates by atmost D i > , D i = sup x ,x ,...,x i , ˆ x i ,...,x N | F ( x , x , . . . , x i , . . . , x N )) − F ( x , x , . . . , ˆ x i , . . . , x N )) | . (8) Then the probability that F deviates from its expected value is given by, P [ | F ( x , x , . . . , x N ) − E [ F ( x , x , . . . , x N )] | ≥ ǫ ] ≤ C exp( − c ǫ D ) , (9) where D = qP i D i is called the diameter and C, c are constants.Proof.
See (Boucheron et al., 2013; Tao, 2010).So instead of computing the DOS using the MCMC approach outlined in section 2,one can exploit McDiarmid’s inequality to compute bounds on the histogram of number ofsatisfied (or unsatisfied) clauses. In the setting of the k -SAT problem, the x , x , . . . , x N inthe McDiarmid’s inequality are replaced by the variables x present in the logical formula.That is, x , x , . . . , x N are the unique set of Boolean variables that occur in the formula.The MCMC computation is now replaced by the set of optimization problems for computingthe diameter as shown in Eqn. 8.For ease of presentation, we focus on the 3-SAT problem instead of the generic k -SATformulation. Polynomial time reductions from k -SAT to 3-SAT make this translation non-restrictive. Every k -SAT instance can be converted to a 3-SAT instance by introducingadditional (ancillary) variables. We now show that the diameter computations for the 3-SAT problem give rise to a QUBO (Boros, Hammer, & Tavares, 2007; Rieffel & Polak, 2011)problem that is particularly amenable to quantum annealers (Kochenberger, Hao, Glover,Lewis, Lü, Wang, & Wang, 2014).
4. QUBO formulation for diameter computations
To estimate D i in Eqn. 8 for the 3-SAT setting, consider the following form for the SATformula, φ ( x ) = m ^ i =1 C i = m ^ i =1 ( x i ∨ x i ∨ x i ) , where x i l is the l -th literal in clause C i . It is easy to check that the number of satisfiedclauses can be expressed as:Φ( x ) = m X i =1 C i = m X i =1 ( x i + x i + x i − x i x i − x i x i − x i x i + x i x i x i ) . (10)6hile this expression has a cubic term, the cubic term disappears when computing thediameter in Eqn. (8) (shown later). We now state our central result that formulates theestimation of diameters D i needed for computing the DOS as a quadratic unconstrainedBoolean optimization problem. Theorem.
The diameter D i for the variable x i in McDiarmid’s inequality can be computedby solving the following optimization problem, D i = max x \ x ( i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X p ∈ S + i [1 − x p − x p + x p x p ] − X p ∈ S − i [1 − x p − x p + x p x p ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , where S + i and S − i are the sets of clauses in which x ( i ) appears in direct and negated forms,respectively.Proof. To compute the diameter D i , pick the i -th variable of x denoted as x ( i ) and computethe worst-case variation in the number of satisfied clauses. Since Φ( x ) is a sum over differentclauses, only those clauses that include x ( i ) in their literal set will contribute to the diameter. S + i and S − i are the sets of clauses in which x ( i ) appears in direct and negated formsrespectively, that is, S + i = { p : C p = x ( i ) ∨ x p ∨ x p } ,S − i = { p : C p = ¬ x ( i ) ∨ x p ∨ x p } , are the set of clauses in which the variable x ( i ) appears in the corresponding literal sets aseither x ( i ) or ¬ x ( i ), respectively. Furthermore, because ∨ is commutative, we can assumewithout loss of generality that the variable x ( i ) appears as the first literal of the clause.Therefore, the expression for the number of satisfied clauses in the 3-SAT instance is: ∀ p ∈ S + i Φ p = x ( i ) + x p + x p − x ( i ) x p − x ( i ) x p − x p x p + x ( i ) x p x p , ∀ p ∈ S − i Φ p =(1 − x ( i )) + x p + x p − (1 − x ( i )) x p − (1 − x ( i )) x p − x p x p + (1 − x ( i )) x p x p . (11)Let S i = { , . . . , m } \ ( S + i ∪ S − i ) be the set of clauses within which x ( i ) does not occur,thus, Φ( x ) = X p ∈ S i Φ p + X p ∈ S + i Φ p + X p ∈ S − i Φ p . The first term in the above sum is not affected by changing x ( i ) and does not contribute tothe diameter and cancels in the subtraction in Eqn. 8. Now, since x ( i ) can only take one oftwo values, { , } , the number of clauses satisfied by setting x ( i ) to 1 in S + i is X p ∈ S + i [1], and7n S − i is X p ∈ S − i x p + x p − x p x p (computed using Eqn. 11). Symmetrically, the number ofclauses satisfied by setting x ( i ) to 0 in S + i is X p ∈ S − i [1], and in S + i is X p ∈ S + i x p + x p − x p x p . D i is the maximum deviation between the two, that is, D i = max x \ x ( i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X p ∈ S + i [1] + X p ∈ S − i [ x p + x p − x p x p ] − X p ∈ S + i [ x p + x p − x p x p ] − X p ∈ S − i [1] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . We get the following optimization problem to compute D i by collecting the terms for sum-ming over S + i and S − i , D i = max x \ x ( i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X p ∈ S + i [1 − x p − x p + x p x p ] − X p ∈ S − i [1 − x p − x p + x p x p ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (12)This result makes sense intuitively, because the expression inside each bracket is logicallyequivalent to ¬ ( x p ∧ x p ), and if either of the other literals is true, the disjunctive clause C p remains false regardless of x ( i ), and therefore does not contribute to the diameter.The expression inside the absolute value in (12) is a quadratic form in x \ x ( i ). Notethat Eqn. 12 can easily be cast into a purely quadratic form x T Q x as the linear terms canbe absorbed into the diagonal of the matrix because x p = x p for binary variables. Remark 1.
The computation for D i involves the maximization of an absolute value. To ad-dress the absolute value, we simply perform two separate maximizations as follows sup x | f ( x ) | =max { sup x f ( x ) , − sup x f ( x ) } . Thus, we compute the two maximizations and choose thelarger result to obtain the diameter. Note that for a -SAT instance with N unique vari-ables, one needs to perform N optimizations. Remark 2.
Besides providing a novel approach for estimating the DOS of k -SAT problems,the diameter computation can be used to benchmark optimization algorithms and computingplatforms. In particular, by comparing the value of the computed diameter by different ap-proaches, one can quantify their performance. Higher diameter values correspond to “better”solutions of the optimization problem. Remark 3.
Using the density of states, one can extract the probability of a random as-signment being at least ǫ (in terms of energy or number of clauses satisfied) away from the verage value ¯ E = P me =0 p ( e ) e . Thus, P h (cid:12)(cid:12)(cid:12) E − ¯ E (cid:12)(cid:12)(cid:12) ≥ ǫ i = m X ¯ E + ǫ p ( e ) . (13) This quantity can be computed numerically from d ( e ) . In our experiments, we solved this QUBO formulation using quantum and classicalcomputing methods. We describe these in detail in the next section.
5. Results
In our experiments, we attempt to answer the following questions.• How does the proposed DOS approach of using concentration of measures and QUBOcompare with the baseline MCMC-FlatSAT approach (Ermon et al., 2010, 2011)?• How do the quantum and classical implementations of the proposed DOS approachcompare with one another?To analyze the performance of the proposed approach, we generate 20 random 3-SAT in-stances for every possible combination of the following sizes ( N = [30 , , , , , , α = [4 . , . , . , α values span thephase transition at α ≈ .
24 that demarcates “easy” and “hard” instances of the satisfiabilityproblem (Monasson et al., 1999). Thus, in total, we generate 560 random 3-SAT instances,and for each instance we compute the baseline DOS using MCMC-FlatSAT. These resultsare then compared with the proposed concentration of measure inequality approach. Ad-ditionally, the 2 N diameter computations for each instance are performed classically usingthe HFS algorithm and the D-Wave quantum annealer. The performance of the D-Wavedevice is then compared with the classical results. We also implemented an SMT-basedoptimization approach on classical platforms, and compared the D-Wave results with thestandard classical solver. We implemented the MCMC-FlatSAT (Ermon et al., 2010, 2011) algorithm in C++. De-pending on the mixing time of the Markov chain (Levin & Peres, 2017), there was a largevariation in the performance of the code. The computation time ranged from hours to sev-eral days (in some instances the code took 3 −
10 days to converge). The computations wereperformed for all of the 560 instances as outlined above. A few instances of the resultingdensity of states are shown in Fig. 2. In general, the MCMC approach was found to havea high computational cost.
We used the Z3 SMT solver (Bjørner et al., 2015) to encode the QUBO problem as abitvector problem exploiting the fixed range of discrete values that can be taken by the9 50 100 150020406080 Number of unsat clauses D e n s i t y o f s t a t e s ( L og2 S c a l e ) (a) D e n s i t y o f s t a t e s ( L og2 S c a l e ) (b) D e n s i t y o f s t a t e s ( L og2 S c a l e ) (c) Figure 2: Density of states for random satisfiability instances with varying size and clausedensity. The X-axis is the number of unsat clauses, Y-axis is the DOS showingnumber of assignments in log scale. (a) N = 125 , α = 5 .
0, (b) N = 150 , α = 4 . N = 200 , α = 4.diameters. The resulting problem is a pseudo-Boolean optimization problem that we solveiteratively using satisfiability solving by binary search (between 0 and the number of clausesin which the variable occurs) over the optimization goal. We allow the SMT solver atimeout of 100 seconds for every trial to find a larger diameter. The 560 instances took8 days to compute. The SMT solver found better solutions for the QUBO compared toquantum annealing, and hence placed more accurate bounds on the histogram. However,the scalability declined sharply with the increase in the number of variables. In particular,we found that the Z3 solver took multiple hours to days to complete several instances. Thediameters computed using the SMT solver can be seen in Fig. 6 for a few values of N and α . 10 -80 -60 -40 -20 0 20 40 60 80 P | E − ¯ E | ≥ | ǫ | $ -20 -10 N = 125 , α = 5 .
00, Instance 6, D = 97 . Numerical DataFit with c=71.23 (a) ǫ -100 -50 0 50 100 P | E − ¯ E | ≥ | ǫ | $ -15 -10 -5 N = 150 , α = 4 .
25, Instance 3, D = 94 . Numerical DataFit with c=59.26 (b) ǫ -100 -50 0 50 100 P | E − ¯ E | ≥ | ǫ | $ -20 -10 N = 150 , α = 4 .
25, Instance 15, D = 94 . Numerical DataFit with c=64.98 (c) ǫ -150 -100 -50 0 50 100 150 P | E − ¯ E | ≥ | ǫ | $ -20 -10 N = 200 , α = 4 .
00, Instance 9, D = 104 . Numerical DataFit with c=57.59 (d)
Figure 3: Comparison of the density of states computed by MCMC and the concentrationof measure bounds. 11 lause Density ( α ) O p t i m a l c N=30N=50N=75N=100N=125N=150N=200
Figure 4: Concentration of measure constant as a function of α . The different lines corre-spond to different values of N. We use the D-Wave 2X (DW2X) annealer (Johnson, Amin, Gildert, Lanting, Hamze, Dick-son, Harris, Berkley, Johansson, Bunyk, et al., 2011) located at the USC InformationScience Institute in Marina del Rey as our quantum platform for computations. ThisDW2X processor is an 1152-qubit quantum annealing device made using superconductingflux qubits (Bunyk, Hoskinson, Johnson, Tolkacheva, Altomare, Berkley, Harris, Hilton,Lanting, Przybysz, & Whittaker, 2014). It has 1098 functional qubits that function at 12mK. The annealer implements the transverse Ising Hamiltonian, H ( s ) = A ( s ) X i σ xi + B ( s ) X i h i σ zi + X ij J ij σ zi σ zj , (14)where s = t/t f is the normalized time, t f is the total evolution time, and A ( s ) and B ( s )are the annealing schedules that modulate the transverse field and Ising field strength,respectively. The total annealing time t f can be set in the range [5 , µ s. The couplingstrengths J ij between qubits i and j can be set in the range [ − , h i can be set in the range [ − , A (0) ≫ B (0) and the system starts in thesuperposition of all possible computational states. During the evolution from s = 0 to s = 1, the transverse field is reduced and the Ising field strength is increased such that A (1) ≪ B (1). If t f is large enough, the adiabatic theorem (Born & Fock, 1928) guaranteesthat the final state of the system will be the ground state of H ( s = 1). The device hasbeen used for machine learning (Biamonte, Wittek, Pancotti, Rebentrost, Wiebe, & Lloyd,2017; Adachi & Henderson, 2015; Mott, Job, Vlimant, Lidar, & Spiropulu, 2017), imagerecognition (Neven, Rose, & Macready, 2008), and combinatorial optimization (Ushijima-Mwesigwa, Negre, & Mniszewski, 2017; McGeoch & Wang, 2013; Neukart, Compostella,Seidel, Von Dollen, Yarkoni, & Parney, 2017; Venturelli, Marchand, & Rojo, 2015) to namea few.We use the above platform to compute the diameters for all the 560 instances of randomsatisfiability problems and compared the results to MCMC-FlatSAT. As noted in remark 1,each instance of an N -dimensional 3-SAT problem gives rise to 2 N optimizations for D i .We chose the smallest possible annealing time t f = 5 µ s. For each QUBO instance of thisstudy, we did 1000 readouts with 10 gauge transforms each (Boixo, Rønnow, Isakov, Wang,12 iameter by HFS
69 70 71 72 73 74 D i a m e t e r b y D - W a v e (a) Instances for N = 100 , α = 4 . Diameter by HFS
95 96 97 98 99 100 101 102 D i a m e t e r b y D - W a v e (b) Instances for N = 125 , α = 5 . Diameter by HFS
93 93.5 94 94.5 95 95.5 96 96.5 97 D i a m e t e r b y D - W a v e (c) Instances for N = 150 , α = 4 . Diameter by HFS
103 103.5 104 104.5 105 105.5 106 106.5 107 D i a m e t e r b y D - W a v e (d) Instances for N = 200 , α = 4 . Figure 5: Comparison of diameters computed by the D-Wave quantum annealer and HFSalgorithm. 13 iameter by SMT
94 95 96 97 98 99 D i a m e t e r b y D - W a v e (a) Instances for N = 150 α = 4 . Diameter by SMT
104 105 106 107 108 D i a m e t e r b y D - W a v e (b) Instances for N = 200 α = 4 . Figure 6: Comparison of diameters computed by the D-Wave quantum annealer and SMTsolver.Wecker, Lidar, Martinis, & Troyer, 2014). Additional details of this particular processcan be found in, for example, Ref. (Mishra, Albash, & Lidar, 2018). Note that the totalwall clock time to optimize each instance, which includes overheads such as initializing thequbits and measurements, was ≈ . ≈ ,
000 instances)coupled with limited affiliate time on the DW2X annealer.We map each diameter computation to a QUBO. Q i ( ~x ) = X p ∈ S + i [1 − x p − x p + x p x p ] − X p ∈ S − i [1 − x p − x p + x p x p ] . (15)The size of this QUBO problem depends on the size of the sets S + i and S − i (see Section 4of the paper). If the SAT problem has N variables and α clause density, the number ofclauses M = N α is a loose upper bound on the size of these QUBO problems. Since theclauses can contain arbitrary variables, for which DW2X has a finite connectivity graph,we need to find a minor embedding of the D-Wave graph that can fit this QUBO prob-lem (Choi, 2008, 2011). In such embedding, each x p i in Eqn. 15 is represented by a chain ofphysical qubits connected via an ferromagnetic couplings. We used the sapiFindEmbedding function provided by the D-Wave application program interfaces (API) to find such embed-dings. We used the sapiEmbedProblem function to submit the jobs to the processor and14 apiUnembedAnswer function with the minimize_energy option to optimally decode theembeddings back to the variables ~x . We used the heuristic ferromagnetic chain couplingprovided by the API. To find higher-quality solutions, one can optimize this ferromagneticcoupling value such that the chain of physical qubits representing each variable is consistentat the end of the anneal. Thus, our results provide a lower bound on the diameter. Po-tentially, one may be able to obtain improved results by performing the actions suggestedabove and optimizing the annealing process.After computing all the D i ’s for a given instance, we can plot the concentration of mea-sure bounds for the DOS. Note that in McDiarmid’s inequality (Eqn. 9), two constantsappear that can be used to make the bounds on the DOS tight. In particular, we find C = 1 and c = 56 . − .
08 exp( − . N − . . α − .
46) give rise to very closeapproximations of the density of states in the range of 30 ≤ N ≤
200 (as shown in Fig. 3).These parameters were computed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) al-gorithm (Liu & Nocedal, 1989). The functional form for c was obtained empirically byfinding the best c for each of the 560 instances and performing regression with respect to N and α (see Fig. 4). We repeat the D i computations for each of the 560 instances by forcing the classical device tocompute the best possible solution using the HFS algorithm. We again run the computationfor 0 . D i . The best D i is saved and the rest arediscarded. We then compare the D i values obtained using quantum annealing with thosecomputed classically. Note that higher diameter values correspond to “better” solutionsas they correspond to higher quality solutions of the QUBO. Note we intend to conducta comprehensive benchmarking study for the D-Wave quantum annealer (Albash & Lidar,2018) (using our DOS framework) in future work.Out of the 560 random satisfiability instances, the D-Wave quantum computer computeshigher quality solutions (higher diameter values) in 306 instances. However, as shown inFig. 5, the D-Wave provided a marginal improvement on the diameter values. In partic-ular, we found that the average solution computed by the D-Wave machine was around0 .
7% higher than the HFS algorithm. The most favorable result for the D-Wave was thecomputation of a solution that was 7% better than the HFS algorithm. Whether this im-provement holds for larger instances remains to be seen and will be tested in higher qubitsettings. We would like to point out, however, that in no instance did the HFS algorithmfind a higher quality solution when compared to the D-Wave machine. As shown in Fig. 6,the SMT solver does find significantly better solutions than the D-Wave machine. Notethat the computational cost of the solver is significantly higher (taking hours to days tocompute the diameter of some instances).
6. Conclusion
In this work, we have developed a novel concentration of measure inequalities-based ap-proach for estimating the density of states of the k -SAT problem. Existing state-of-the-artMarkov Chain Monte Carlo methods for computing density of states are stymied by compu-tational intractability. Our approach provides estimates for the density of states histogram15y converting the problem into a set of optimization problems that bound the maximumvariability of the cost function (called “diameters”). For the 3-SAT, these optimizationproblems elegantly reduce to quadratic unconstrained binary optimizations, thereby makingthem amenable for commercially available quantum annealers such as the D-Wave machine.To test our approach, we computed diameters for a range of random instances that spanthe phase transition of the 3-SAT. We used the D-Wave quantum annealer and computedbounds on the density of states for those problem instances. We find that by tuning a singleparameter (exponent), one can accurately estimate the entire density of states histogram ,orders of magnitude faster than state-of-the-art Markov Chain Monte Carlo methods. Addi-tionally, we compared the D-Wave diameter computations to equivalent classical techniquessuch as SMT solvers and the HFS algorithm. We find that the D-Wave machine providesa marginal improvement in the diameter computations over the HFS algorithm. In otherwords, the D-Wave annealer computes solutions that correspond to larger values of thediameter. The SMT solver works particularly well for small problem instances. However,we find this approach does not scale well and is not competitive (from a computationalstandpoint) for larger instances of the 3-SAT.In summary, we propose a new approach for estimating density of states of the k -SATproblem that can be implemented on both classical and quantum platforms. Moreover, theproblem is a particularly interesting test for comparing quantum platforms (annealers andother noisy intermediate-scale quantum devices) to classical computation. This is becausethe diameter values provide a real number metric for comparison. In other words, the qualityof solution is a real number as opposed to the standard satisfiability tests for benchmarkingquantum devices that yield inconclusive results of the form “no satisfiable assignmentsfound” for most instances. We hope this problem and the outlined approach can be used toanalyze complex aerospace system requirements from a satisfiability standpoint as well asto test emerging quantum platforms against their classical counterparts. DOS estimationcan be used for probabilistic inference and, thus, an efficient quantum algorithm for thedensity of states estimation will enable development of quantum artificial intelligence.In future work, we plan to expand the use of concentration of measure inequalities fora wider set of problems and compare the above approaches with simulated annealing-basedmethods (Aarts & Korst, 1988).
7. Acknowledgements
The authors thank Federico Spedalieri and Daniel Lidar of the University of SouthernCalifornia-Information Sciences Institute (USC-ISI) for discussions and suggestions. Theauthors thank Lockheed Martin’s quantum computing team (Christopher Elliott, Greg Tal-lant, and Kristen Pudenz) for discussions related to the approach and generously providingaffiliate access to their D-Wave quantum annealer. Dr. Jha acknowledges support from theUS Army Research Laboratory Cooperative Research Agreement W911NF-17-2-0196, andNational Science Foundation (NSF) grants eferences
Aarts, E., & Korst, J. (1988).
Simulated annealing and Boltzmann machines . New York,NY; John Wiley and Sons Inc.Abu-Mostafa, Y. S., Magdon-Ismail, M., & Lin, H.-T. (2012).
Learning from data , Vol. 4.AMLBook New York, NY, USA.Adachi, S. H., & Henderson, M. P. (2015). Application of quantum annealing to trainingof deep neural networks. arXiv preprint arXiv:1510.06356 .Albash, T., & Lidar, D. A. (2018). Demonstration of a scaling advantage for a quantumannealer over simulated annealing.
Physical Review X , (3), 031016.Biamonte, J., Wittek, P., Pancotti, N., Rebentrost, P., Wiebe, N., & Lloyd, S. (2017).Quantum machine learning. Nature , (7671), 195.Biere, A., Heule, M., & van Maaren, H. (2009). Handbook of Satisfiability , Vol. 185. IOSpress.Birnbaum, E., & Lozinskii, E. L. (1999). The good old Davis-Putnam procedure helpscounting models.
Journal of Artificial Intelligence Research , , 457–477.Bjørner, N., Phan, A.-D., & Fleckenstein, L. (2015). ν z-an optimizing SMT solver. In International Conference on Tools and Algorithms for the Construction and Analysisof Systems , pp. 194–199. Springer.Boixo, S., Rønnow, T. F., Isakov, S. V., Wang, Z., Wecker, D., Lidar, D. A., Martinis, J. M.,& Troyer, M. (2014). Evidence for quantum annealing with more than one hundredqubits.
Nat. Phys. , (3), 218–224.Born, M., & Fock, V. (1928). Beweis des Adiabatensatzes. Zeitschrift für Physik , (3-4),165–180.Boros, E., Hammer, P. L., & Tavares, G. (2007). Local search heuristics for quadraticunconstrained binary optimization (qubo). Journal of Heuristics , (2), 99–132.Boucheron, S., Lugosi, G., & Massart, P. (2013). Concentration inequalities: A nonasymp-totic theory of independence . Oxford University Press.Bunyk, P. I., Hoskinson, E. M., Johnson, M. W., Tolkacheva, E., Altomare, F., Berkley,A. J., Harris, R., Hilton, J. P., Lanting, T., Przybysz, A. J., & Whittaker, J. (2014).Architectural considerations in the design of a superconducting quantum annealingprocessor.
IEEE Trans. Appl. Supercond. , (4), 1–10.Chieu, H. L., & Lee, W. S. (2009). Relaxed survey propagation for the weighted maximumsatisfiability problem. Journal of Artificial Intelligence Research , , 229–266.Choi, V. (2008). Minor-embedding in adiabatic quantum computation: I. The parametersetting problem. Quantum Inf. Process. , (5), 193–209.Choi, V. (2011). Minor-embedding in adiabatic quantum computation: II. Minor-universalgraph design. Quantum Inf. Process. , (3), 343–353.Cook, S. A. (1971). The complexity of theorem-proving procedures. In ACM Symposiumon Theory of Computing , pp. 151–158. ACM.17ostello, R. J., & Liu, D.-B. (1995). Metrics for requirements engineering.
Journal ofSystems and Software , (1), 39–63.Dutertre, B. (2014). Yices 2.2. In International Conference on Computer Aided Verification ,pp. 737–744. Springer.Ermon, S., Gomes, C., & Selman, B. (2011). A flat histogram method for computing thedensity of states of combinatorial problems. In
International Joint Conference onArtificial Intelligence .Ermon, S., Gomes, C. P., & Selman, B. (2010). Computing the density of states of Booleanformulas. In
International Conference on Principles and Practice of Constraint Pro-gramming , pp. 38–52. Springer.Ferrante, O., Scholte, E., Pinello, C., Ferrari, A., Mangeruca, L., & Liu, C. (2016). Amethodology for increasing the efficiency and coverage of model checking and its ap-plication to aerospace systems.
International Journal of Aerospace , (2016-01-2053),140–150.Gourgoulias, K., Katsoulakis, M. A., Rey-Bellet, L., & Wang, J. (2017). How biased isyour model? Concentration inequalities, information and model bias. arXiv preprintarXiv:1706.10260 .Hamze, F., & de Freitas, N. (2004). From fields to trees. In Uncertainty in ArtificialIntelligence , pp. 243–250. AUAI Press.Johnson, L. A., et al. (1998). DO-178B, software considerations in airborne systems andequipment certification.
Crosstalk , .Johnson, M. W., Amin, M. H., Gildert, S., Lanting, T., Hamze, F., Dickson, N., Harris,R., Berkley, A. J., Johansson, J., Bunyk, P., et al. (2011). Quantum annealing withmanufactured spins. Nature , (7346), 194.Kirkman, D. (1998). Requirement decomposition and traceability. Requirements Engineer-ing , (2), 107.Kochenberger, G., Hao, J.-K., Glover, F., Lewis, M., Lü, Z., Wang, H., & Wang, Y. (2014).The unconstrained binary quadratic programming problem: a survey. Journal ofCombinatorial Optimization , (1), 58–81.Krentel, M. W. (1988). The complexity of optimization problems. Journal of Computerand System Sciences , (3), 490–509.Leveson, N. G., & Weiss, K. A. (2009). Software system safety. In Safety Design for SpaceSystems , pp. 475–505. Elsevier.Levin, D. A., & Peres, Y. (2017).
Markov chains and mixing times , Vol. 107. AmericanMathematical Soc.Leyendecker, S., Lucas, L. J., Owhadi, H., & Ortiz, M. (2010). Optimal control strategiesfor robust certification.
Journal of Computational and Nonlinear Dynamics , (3),031008.Liu, D. C., & Nocedal, J. (1989). On the limited memory BFGS method for large scaleoptimization. Mathematical Programming , (1-3), 503–528.18cDiarmid, C. (1989). On the method of bounded differences. Surveys in Combinatorics , (1), 148–188.McGeoch, C. C., & Wang, C. (2013). Experimental evaluation of an adiabiatic quantum sys-tem for combinatorial optimization. In ACM International Conference on ComputingFrontiers , p. 23. ACM.Mishra, A., Albash, T., & Lidar, D. A. (2018). Finite temperature quantum annealingsolving exponentially small gap problem with non-monotonic success probability.
Nat.Commun. , (1), 2917.Monasson, R., Zecchina, R., Kirkpatrick, S., Selman, B., & Troyansky, L. (1999). De-termining computational complexity from characteristic ‘phase transitions’. Nature , (6740), 133.Mott, A., Job, J., Vlimant, J.-R., Lidar, D., & Spiropulu, M. (2017). Solving a Higgs opti-mization problem with quantum annealing for machine learning. Nature , (7676),375.Neukart, F., Compostella, G., Seidel, C., Von Dollen, D., Yarkoni, S., & Parney, B. (2017).Traffic flow optimization using a quantum annealer. Frontiers in ICT , , 29.Neven, H., Rose, G., & Macready, W. G. (2008). Image recognition with an adiabaticquantum computer I. Mapping to quadratic unconstrained binary optimization. arXivpreprint arXiv:0804.4457 .Owhadi, H., Scovel, C., Sullivan, T. J., McKerns, M., & Ortiz, M. (2013). Optimal uncer-tainty quantification. SIAM Review , (2), 271–345.Preskill, J. (2018). Quantum computing in the NISQ era and beyond. Quantum , , 79.Rieffel, E. G., & Polak, W. H. (2011). Quantum Computing: A Gentle Introduction . MITPress.Sebastiani, R., & Trentin, P. (2015). OptiMathSAT: A tool for optimization modulo theories.In
International Conference on Computer Aided Verification , pp. 447–454. Springer.Selby, A. (2013). Qubo-chimera.
Online: https://github. com/alex1770/QUBO-Chimera[Accessed 28 October 2019] .Selby, A. (2014). Efficient subgraph-based sampling of Ising-type models with frustration. arXiv preprint arXiv:1409.3934 .Sommerville, I. (2005). Integrated requirements engineering: A tutorial.
IEEE Software , (1), 16–23.Tamura, N., Taga, A., Kitagawa, S., & Banbara, M. (2009). Compiling finite linear CSPinto SAT. Constraints , (2), 254–272.Tao, T. (2010). 254A, Notes 1: Concentration of measure. What’s New (Blog) .Tao, T. (2012).
Topics in random matrix theory , Vol. 132. American Mathematical Soc.Ushijima-Mwesigwa, H., Negre, C. F., & Mniszewski, S. M. (2017). Graph partitioningusing quantum annealing on the D-Wave system. In
International Workshop on PostMoores Era Supercomputing , pp. 22–29. ACM.19enturelli, D., Marchand, D. J., & Rojo, G. (2015). Quantum annealing implementation ofjob-shop scheduling. arXiv preprint arXiv:1506.08479 .Walsh, T. (2000). SAT v. CSP. In
International Conference on Principles and Practice ofConstraint Programming , pp. 441–456. Springer.Wang, F., & Landau, D. (2001). Efficient, multiple-range random walk algorithm to calcu-late the density of states.
Physical Review Letters , (10), 2050.Xu, K., & Li, W. (2000). Exact phase transitions in random constraint satisfaction problems. Journal of Artificial Intelligence Research ,12