On using floating-point computations to help an exact linear arithmetic decision procedure

On using fl oating-p oin t computations to help an exact linear a r ithmetic decision pro cedure ∗ Da vid Monniaux CNRS / VERIMA G † Octob er 29 , 20 18 Abstract W e consider the decision problem for q uantifier-free form ulas whose atoms are linear inequalities interpreted o ver the reals or rationals. This problem ma y b e decided using satisfiabilit y mo dulo theory (SMT), using a mixture of a SA T solver and a simplex-based decision pro cedure for conjunctions. State-of-th e-art SMT solvers use simplex implementations o ver rational num b ers, whic h p erform well for typical p roblems arising from model-chec king and program analysis (sparse i neq u alities, small co- efficien ts) but are slo w for other a pp lications (denser problems, larger coefficients). W e prop ose a simple prepro cessing phase th at can b e adapted to ex- isting SMT solvers and that may b e optionally triggered. Despite us- ing floating-point computations, our method is sound and complete — it merely affects efficiency . W e implemen ted th e method a nd provide bench- marks sho wing that this change brings a naiv e and slow decision p rocedu re (“textb o ok simplex” with rational num b ers) up to the efficiency of recent SMT solv ers, o ver test cases arising from mo del-chec king, and makes it definitely faster th an state-of-the-art S MT solvers on dense ex amples. 1 In tro duction Decision pro ce dur es for ar ithmetic theories are widely us ed for c o mputer-aided verification. A decision pro ce dur e for a theor y T takes as input a formula of T and outputs a Bo olean: whether the fo r mu la is sa tisfiable. F or many decidable and potentially useful theories, ho wev er, decision pro c edures a r e sometimes to o slow to pro cess problems b eyond small exa mples . This is for insta nc e the ca se of the theo ry o f r e a l clos e d fields (polynomial arithmetic over the real num b ers). Excessive computation times ar ise fr om tw o so urces: the Bo olean structure of the formulas to b e decided (prop ositional satisfiability is cur rently so lved in ex- po nent ial time in the worst cas e), and the in trinsic har dness o f the theo ry . In recent y ear s, SA T mo dulo theory (SMT) techniques hav e a ddressed the former source of inefficiency , by leveraging the power of efficient SA T (Bo olean sa tisfia- bilit y) solvers to deal with the Bo o lean structure. SMT solvers com bine a SA T ∗ This work was partiall y funded by the ANR AR P EGE pro ject “ASOPT ”. † VERIMAG is a joi n t lab oratory of CNRS, Unive rs it ´ e Joseph F ourier and Grenoble-INP. 1 solver with a de c is ion pr o cedure for co njunctions o f a toms in T . If T is linear real arithmetic (LRA), then this decision pro cedure must decide whether a set of linea r inequalities with rational or integer co efficient s has rational so lutions. The pro ble m of testing whether a set of linea r inequa lities has a solution and, if it has, to find a solution tha t maximizes s o me linear co m bination of the v ar iables is kno wn as linear progr amming and has been considera bly studied in op erational re search. V ery efficient implement atio ns exist, whether co mmercial or not, and are able to so lve very larg e pr oblems. They a re not dir ectly appli- cable to our pro blems, how ever, if o nly b ecause they op erate over floating-p oint nu mbers and provide in gener a l no a s surance that their r esult is truthful, de- spite elab orate precautio ns ta ken ag a inst n umerical instabilities. As a result, the decision pro cedures for LRA in SMT solvers are implement ed with rationa l arithmetic, which is slow er than floating-p oint, especia lly if co efficients beco me large, as often happ ens with dense linear pr o blems: large c o e fficien ts fo r ce the use of costly extended precision arithmetic. It thus would seem desirable to leverage the sp eed a nd ma turit y of floating- po in t linear pro gramming s ystems to enhance exac t decision pr o cedures. This article des c rib e s a simple pr epro cessing phase that can b e a dded, with minimal change, to existing r ational simplex implement atio ns used as dec ision pro cedures inside SMT solvers. The pro cedure was implemen ted on top o f a naive and inefficient rational simplex implemen tation; the resulting pro cedure riv a ls recent SMT so lvers. A similar metho d ha s been prop osed in the op e rational re search field [4 ], 1 but there are rea sons why it may perfor m less well fo r the typical optimization tasks of ope rational resear ch than for decision tasks. The nov elty of this article is the application of this technique as a simple mo dification of exis ting SMT algorithms. 2 Simplex SMT so lvers ne e d a decisio n pr o c edure capable o f: • b eing used incr ementally: adding new constraints to the problem, and r e- moving blo cks of constraints, prefera bly witho ut recomputing ev erything; • telling whether the pro blem is sa tisfiable or not; • if the problem is unsatisfiable, outputting a (preferably small o r even min- imal) unsatisfiable subset of the constra in ts; • propa gating theory lemmas, if p ossible at reasonable costs (fro m a con- junction C 1 ∧ . . . ∧ C n , o btain literals L 1 , . . . , L m that a re co ns equences of that conjunction: C 1 ∧ . . . ∧ C n ⇒ L 1 ∧ . . . ∧ L m ). All cur rent SMT so lvers seem to decide general linear real arithmetic (as opp osed to sy ntactic restrictions thereof suc h as difference log ic) us ing the sim- plex algorithm . This algor ithm is exp onential in the worst ca se, but tends to per form w ell in practice; none of the current s olvers seem to use a (p olyno mial- time) in terior p oint metho d. Our metho d is a v ariant of the simplex a lgorithm; we shall thus first descr ibe the “conven tio nal” simplex. 1 W e were not aw are of this article when w e started our work, and we thank Bernd G¨ artner for p oint ing i t to us. 2 2.1 Basic simplex W e shall first give a brief summary on the dual simplex algor ithm on which the LRA decis ion pro cedures in Yices 2 [6, 5] and Z3 3 [3] are ba sed. There o therwise exist many exce llen t textb o o ks on the simplex algorithm [2 , 1 5], though these seldom discuss the sp ecifics of implementations in exact pr ecision or incremental use. T ake a system of linea r equatio ns, e.g.    x − 2 y ≤ 1 − y + 3 z ≥ − 1 x − 6 z ≥ 4 (1) The s y stem is firs t made cano nical. Inequalities ar e scaled so that each left hand side o nly has integer co efficients w ith no co mmon factors. Then, each inequality is optionally negated so that the first coefficient app earing (using some arbitrary or dering of the v aria bles) is pos itive. This ensures that tw o inequalities c o nstraining the same direction in space (e.g. − y + 3 z ≥ − 1 and 2 y − 6 z ≥ 3) app ear with the exa ct sa me left-ha nd side. F or each left-hand side that is not a v a riable, a new v ariable is in tro duced; the system is then conv erted int o a num b er of linea r e q ualities and b ound constraints on the v ariables. F or instance, the ab ov e s y stem gets c o nv erted int o:    α = x − 2 y β = y − 3 z γ = x − 6 z    α ≤ 1 β ≤ 1 γ ≥ 4 (2) The pr o blem is thus formulated as deciding whether a pro duct of int er v als in- tersects a linear subspace given b y a basis. The se t of v aria bles is partitioned in to b asic and nonb asic v ar ia bles; the nu mber of basic v ariables s tays co nstant thr o ughout the algo r ithm. Ba sic v ar i- ables ar e expres sed as linear co m binatio ns of the no n bas ic v aria bles. The main op eration o f the algo rithm is pivoting : a basic v ar iable is made no n bas ic and a nonbasic v a riable is made basic, without c hang ing the linear subspace defined by the system o f equations. F or instance, in the ab ov e example, α = x − 2 y defines the basic v aria ble in ter m o f the nonbasic v a riables x and y . If one wants instead x to b e made basic, one o bta ins x = α + 2 y . The v ariable x then has to be r eplaced b y α + 2 y in all the other equa lities, so that the right-hand sides of these equalities only refer to nonbasic v aria ble s . This replac e men t pro cedure (essentially , replacing a vector u by v + k u ) is the costly part of the alg orithm. A mor e formal description is given in Alg . 1. Let us insist that pivoting do es not change a n ything to the v alidity o f the problem: b oth the bo unds a nd the linear subspace stay the sa me. The idea behind the simplex a lgorithm is to pivot until a p osition is found where it is obvious that the problem has a s olution, or that it ha s not. The algo rithm also maintains a vector of “ current” v alues for the v ariables . This vector is fully defined b y its pro jection on the nonbasic v ariables , since the basic v ar iables ca n b e obtained from them. The curre nt v alues of the no n basic v ar iables always stay within their resp ective b ounds. If all cur rent v alues of the 2 http://y ices.csl.sri.c om/ 3 http://r esearch.micros oft.com/ en- us/um/redmond/projects/z3/ 3 basic v ariables a lso fit within the b ounds, the current v alue is a solution to the problem and the algo rithm terminates. If there are bas ic v aria bles that fall outside the prescr ibe d b ounds, one of them (say , α ) is selected and the corresp onding r ow (sa y , x − 2 y ) is exa mined. Suppo se for the sak e o f the explanation that the cur rent v alue f or α , c α , is strictly grea ter than the max imal prescrib ed v alue M α . One can tr y making x smaller or y larger to comp ensate for the difference. If x is alrea dy at its low er b ound and y at its uppe r b ound, it is impos sible to make α smalle r and the system is unsatisfiable; the algo rithm then t er mina tes. In other words, by p erfo rming interv al arithmetic ov er the equa tion defining α in terms of the nonbasic v ariables , one shows this equa tion to b e unsatisfia ble (replacing the nonbasic v ar iables by their in terv al b ounds, o ne obtains an in terv al that do es not intersect the interv al for α ). Let us now s uppos e that x is no t at its low er b ound; we can try making it smaller. α a nd x a re pivoted: x bec omes basic and α no n basic. α is set to its low er b ound a nd x is a djusted acco rdingly . It can b e sho wn that if ther e is a s olution to the problem, then there is one co nfiguration where any nonbasic v a riable is either at its low er o r upp er bo und. In tuitively , if some nonbasic v aria bles are “in the middle”, then this means we hav e some “ slack marg in” av a ila ble, so we should as well use it. The simplex algorithm appea r s to mov e betw een nu meric a l p oints taken in an infinite contin uous space, but in fa ct, its cur rent config uration is fully defined by stating which v ariables are ba sic and nonbasic, and, for the nonbasic v ariables, whether they a re at their lower or upp er b o und — thus the nu mber of configuratio ns is finite. Remark that we left aside how we cho ose which ba sic v ar iable and which basic v aria ble to pivot. It can b e shown that certain piv oting strategies (say , choose the suitable non basic a nd basic v ariables of least index) necess arily lead, maybe after a gr eat num b er o f pivots, to a co nfiguration where either a solution is obtained, or the constraints of a nonbasic v ariable are clear ly unsatisfiable . The idea of our a rticle is based o n the f ollowing r emark: the simplex al- gorithm, with a s uitable piv oting strategy , a lwa ys terminates with the right answer, but its executio n time ca n v ar y consider ably dep ending o n the initial configuratio n. If it is star ted from a configur ation where it is obvious tha t the system has a solution or do es not hav e one, then it terminates immediately . Otherwise, it may do a great deal of pivoting. O ur idea is to use the output of some untrusted floating - po in t simplex algor ithm to direct the rational simplex to a hop efully go o d star ting po in t. 2.2 Mo difications and extensions The a bove a lgorithm is a quic k description o f “textb o ok simplex”. It is not suf- ficient for obtaining a numerically sta ble implementation if implemen ted over floating-p oint; this is why , after initial trials with a naive floating-p oint imple- men tation, w e mov ed to a b etter off-the-she lf implementation, namely Gl pk (GNU Linea r Prog ramming Kit) [9]. So far , we hav e only considered wide inequa lities. A classic a l metho d to con- vert problems with str ict inequalities in to problems with only wide inequa lities is to intro duce infinitesimals : a c o o rdinate is no long e r one single rationa l num- ber , but a pair o f rational num b ers a a nd b , denoted a + b ε , with lexicogr aphic 4 ordering. ε is th us a sp ecial num ber , g reater than zero but less than all pos itive rationals. x < y is rewritten into x + ε ≤ y . The “ current” v alues in the simplex algorithm are thus pairs ( a, b ) of ratio nals, noted a + bε , the upp er b ounds ar e either + ∞ or a + bε , the low er b ounds either −∞ e ither a + bε . The “ current” vector, a vector of pairs , can b e equiv alently represented a s a pa ir of vectors of rationals ( u , v ), written u + ε v — if v 6 = 0 , it designates a p oint infinitely c lo se to u in the dire c tio n of v , if v = 0 it only means the p oint u . In a ddition to a yes/no answer, decision pro cedures used in SMT need to provide: • In case o f a p ositive a nswer, a solution p oint. If only wide inequalities ar e used, this is given s traightforw ar dly by the “ current” v a lues. F or p oints of the form u + ε v , one considers the half-line u + t v , t > 0 , inject it into the system of inequa lities, and so lve it for t — the solution set will b e a n int erv al o f the fo r m (0 , t 0 ). A suitable s olution p oint is thus u + t 0 2 v . • In c ase o f a negative answer, a contradiction witness: nonnegative co ef- ficient s such that by multiplying the or ig inal inequalities by thos e co effi- cients, one g ets a tr iv ially unsatisfia ble inequality (0 < c where c ≤ 0, or 0 ≤ c where c < 0 ). This contradiction witness is obtained b y using an auxiliary tableau tracking how the equalities b − P n t b,n n = 0 defining the ba sic v ariables were o btained as linear combinations of the or iginal equalities de fined at the initialization of the simplex, as des c r ibe d in [6, Sec. 3.2.2 ]. 3 Mixed floating-p oin t / rational strategy Our pro cedure tak es as input a rational simplex problem in th e forma t de- scrib ed a t Sec. 2.1: a ta ble a u of linear equalities and b ounds o n the v ariables. It initializes a floating-p oint simplex by converting the ratio nal problem: the infinitesimals are dis carded and the ratio na ls rounded to near est. It then ca lls, as a subroutine, a floating - po in t s implex algorithm which, on exit, indicates whether the pr oblem is satisfiable (a t least acco r ding to floating-p oint compu- tation), a nd a co rresp onding config uration (a partition into basic and nonbasic v ar iables, and, for each no nbasic v ariable , whether it is set to its upper or lo wer bo und). There are several suitable pack ages av ailable; in o rder to p erform ex- per iment s, w e implemented our metho d us ing the GNU Linear pr ogramming to olkit ( G lpk ) [9]. In or der for the resulting mixed procedure to b e used incr emen tally , the floating-p oint solver s hould supp o rt incre mental use. Commercial linear pro- gramming so lvers are desig ned for large pr oblems sp ecified as a whole; there may be a lar ge ov erhea d for lo ading the problem into the s olver, even if the problem is sma ll. Instea d, w e nee d solvers capa ble of incrementally adding and withdrawing constraint bounds at minima l co st. Glpk suppo rts incr emen tal use, s ince it keeps the facto rization of the bas is in memor y b etw een calls [9, p. 2 0 ]; this factorizatio n is costly to co mpute but needs to b e computed only if the bas is matrix changes: in our ca se, this basis stays the sa me. A t this p oint, if successful, and unless there has b een so me fatal numeric degeneracy , the floating-p oint simplex outputs a floating-p oint approximation to a solution p oint. How ever, in genera l, this approximation, converted back 5 int o a ratio nal p oint, is no t neces s arily a true solution point. The reaso n is that simplex pro duces so lution p oints a t a vertex of the solution poly hedron, and, nu meric a lly sp eaking, it is in g e neral imp ossible to b e exactly on that p oint; in general, the solution p oint obtained is found to b e v ery slightly outside of the po lyhedron when membership is tested in ex a ct arithmetic. It is therefore not feasible to simply take this solution p oint as a witness. The ra tional simplex tableau is then “forc e fully piv oted” , using Algorithm 2, un til it ha s the same basic/no n bas ic partition as the floating-p oint output. This amounts to a pass o f Gaussia n elimination for changing the basis of the linear subspace. This phase c an partially fail if the final basic /nonbasic partitio n re- quested is infeasible in exa ct precis ion arithmetic — maybe be c ause of bug s in the floating- po in t simplex algorithm, or simply b ecaus e of floating -po int inac- curacies. The “cur r ent” v a lues o f the nonbasic v ariables are then initialized according to the output of the floating -po int simplex: if the floating- po in t simplex se le cted the upper bound for no n basic v ariable n then its curre n t v alue in the r ational simplex is set to its upp er bo und, a nd similarly for lower bounds. If the final ba- sic/nonbasic partition pr o duced by the floating-p oint simplex is infea sible, then there ar e nonb as ic v aria bles of the rational simplex for which no information is known: these a re left untouc hed o r set to arbitrary v a lues within their b ounds (this do es not affect co rrectness). The cur rent v alues of the bas ic v ariable s are computed using the rational table a u. The ratio na l simplex is then started. If thing s hav e gone well, it ter minates immediately by noticing that it has fo und either a so lution, or a co nfig uration showing that the system is unsa tisfiable. If things have go ne mo dera tely w ell, the ratio nal simplex do es a few additiona l r ounds of pivoting. If things hav e gone ba dly , the ra tional simplex pe rforms a full s earch. The rational simplex alg orithms are well known, and we have already pre- sented them in Sect. 2.1. The corr ectness of o ur mixed alg orithm relies on the correctnes s of the r ational simplex and the “forced piv oting” pha se maintaining the inv a riant that the linear equa lities define the same solutions as those in the initial s ystem. Algorithm 1 Pivot ( table au , b, n ): pivot the basic v ariable b a nd the no n bas ic v ar iable n . t v is the line defining basic v ariable v , t v, w is the co efficient of t v corres p onding to the non basic v a riable w . The a i are the optional auxilia ry tableau descr ibed in Sec. 2 .2. Require: b ba sic, n no n bas ic , t b,n 6 = 0 p := − 1 /t b,n ; t n := p.t b ; t n,n := 0; t n,b := − p ; a n = p.a b B := B ∪ { n } \ { b } for all b ′ ∈ B do p := t b ′ ,n ; t b ′ ,n := 0; t b ′ := t b ′ + p.t b ; a b ′ := a b ′ + p.a b end fo r t b := 0 The “for all” loo p is the most exp ensive pa rt of the whole simplex algor ithm. Note that, dep ending on the way the sparse arrays a nd auxiliary structures are implemen ted, this lo o p may b e parallelized, each iteratio n b eing independent of the other s. This g ives a p erformance b o ost on dense matrices. 6 Algorithm 2 For cedPivot ( table au , B f ): force pivoting un til the se t of ba sic v ar iables is B f B := B i rep eat hasPivote dSomething := false for all b ∈ B \ B f do if ∃ n ∈ B f \ B t b,n 6 = 0 then Cho ose n in B f \ B such that t b,n 6 = 0 Pivot ( table au , b, n ) { This inv olves B := B ∪ { n } \ { b }} hasPivote dSomething := true end i f end fo r un til ¬ hasPivo te dSomething return B = B f W e s hall now descr ibe in more detail the “forced pivoting” a lgorithm (Alg. 2). This algor ithm takes as input a simplex tableau, with ass o ciated partition of the set V of v aria bles into basic ( B i ) and non basic v aria bles ( ¯ B i ), a nd a final partition of basic ( B f ) and nonbasic v ariables ( ¯ B f ). F or each bas ic v ar iable b , the tableau contains a line b = P n ∈ ¯ B t b,n n . Additional constr aints are that the table a u is well-formed (basic v a r iables are combination o f only the nonbasic v ar iables) and tha t | B i | = | B f | (since | B | is a constant). Assuming that all arithmetic o per ations ta ke unit time (whic h is not true when one co nsiders dense problems, since coefficient sizes quickly gr ow, but is almost true for sparse problems with small co efficient s), the r unning time of the forced pivoting algo r ithm is at most cubic in the size of the problem. This mo tiv ates our s uggestion: instead of per forming an exp ensive and p oten- tially expo nen tial- time search dir ectly with rational arithmetic, we per form it in floa ting-p oint, with all the p ossibilities of optimization of floating -po int lin- ear arithmetic offered b y mo dern libr aries and compilers , a nd then p erform a cubic-time pa ss with rationa l a rithmetic. Not a ll final config urations are feas ible: it is p ossible to ask For cedPivot to p erfo r m an imp ossible transfor ma tion, in which case it returns “ false”. F or instance, if the input system o f equations is x = a + b ∧ y = a + b , thus with B i = { x, y } , then it is imp ossible to mov e to B f = { a, b } , for ther e is no way to express a and b a s linea r functions of x and y . More precisely , we say that a configur ation is f ea s ible if it is p ossible to write the basic v ariables of the configuratio n as linear function o f the nonbasic v ariables and obta in the same solutions a s the initial system. Lemma 1. For cedPivot suc c e e ds ( and r eturns “true”) if and only if the final p artition define d by B f is fe asible, otherwise it r eturn s “false”. Pr o of. Let S denote the space of so lutions of the input system of equations. At all iterations, S is e x actly defined by the sy stem of equations, and dim S = | ¯ B | . The only w ay the pro cedure fails is when B 6 = B f and y et, for all basic v aria ble b ∈ B \ B f and nonbasic v ariable n ∈ B f \ B , it is imp ossible to pivot b and n becaus e t b,n = 0. In other words, all s uch b are linear combinations of the nonbasic v a riables in ¯ B ∩ ¯ B f . All v a riables in ¯ B f are thus linear combinations 7 0 5000 10000 15000 20000 25000 0 1000 2000 3000 4000 5000 6000 cumulated time (s) number of benchmarks naive rational Simplex mixed Simplex Yices 1.0.9 Z3 SMT-COMP 2008 Figure 1: Benchmarks on unsatisfiable conjunctions extracte d from vSMT verification prob- lems. Even though our i m plemen tation of spars e ari thmetic and the rational si mplex ar e not up to those i n state-of-the-art s ol v ers (as shown by the lack of p erform ance on the “easy” examples on the left), and our pro cedure is not geared tow ards the sparse problems typical of v erification applications, it still p erforms faster than Yi ces 1. In 4416 cases out of 5736 (77%), no additional simplex pivots are needed after Forc ePivot . of v aria bles in ¯ B ∩ ¯ B f , and since we hav e supp osed that B 6 = B f , ¯ B ∩ ¯ B f ( ¯ B th us | ¯ B ∩ ¯ B f | < | ¯ B | = dim S . But then, | ¯ B f | < dim S and B f cannot b e a feasible config uration. One can still add a few ea s y improv ements: • Befor e embarking into an y simplex, we first test that the o riginal prob- lem do es not contain a trivially unsatisfiable tableau row: one where the b ounds obtained by interv a l arithmetic on the r ight-hand side of b = P n t b,n n the equality hav e an e mpty intersection with thos e for b . • During the forced pivoting pr o cedure, we detect whether the new equality obtained is trivia lly unsatisfiable, in which case we terminate immediately . • F o rced pivots c an b e done in any order. At the b eginning o f the pr o cedure, we sor t the v ar iables in B i \ B f according to the num b er of nonzer o co ef- ficient s in the co rresp onding equation. When choosing the basic v a riable to be pivoted, w e take the lea st one in that order ing . The idea is that the pivoting steps are faster when the equation defining the basic v ariable to be pivoted has few nonze ro co efficients. • Similarly , one can pre-sort the v a riables in B f \ B i according to their nu mber of o ccurr ences in e quations. 8 0 5000 10000 15000 20000 25000 30000 0 10 20 30 40 50 60 70 80 cumulated time (s) number of benchmarks naive rational Simplex mixed Simplex Yices 1.0.9 Z3 SMT-COMP 2008 Figure 2: Benc hmarks on systems of 100 inequalities with 50 v ariables, where all coefficients are taken unifor mly and independently dis tributed integers in [ − 100 , 100 ]. 31 are unsatisfiable, 46 are satisfiable. F or each solv er, b ench marks sorted by increasing time spent. In 58 cases out of 82 (71%), no additional simplex pivots are need ed after For cePivot . The SMT pr o cedure may hav e at its disp osa l, in addition to the “exact” theory test a pa rtial a nd “fast” theory test, whic h may err on the side of sat- isfiability: first test sa tisfiability in flo a ting-p oint, and test in exact ar ithmetic only if negative. The “fast” theory test ma y b e used to test whether it seems a go o d ide a to inv estigate a branch o f the search tree or to backtrack, while the “ exact” results should b e used for theory propa gation or w hen it s eems a go o d idea to chec k whether the cur rent branch truly is satisfia ble . V ar ious sys- tems are p ossible, depending on how the SA T and theory parts interact in the solver [8]. Note that after an “ex act” theory chec k, we hav e a n exa ct r a tional simplex tableau corr esp onding to the constra in ts, from which it is p ossible to extract theory propag ation information. F or instance, if interv al analysis on a r ow x = 3 y + 5 z , using interv als from y and z , sho ws that x < 3, th en one can immediately conclude that in the c ur rent SA T sea rch branch, the liter al x ≥ 4 is false. 4 Implemen tation and b enc hmarks W e implemented the above algo r ithm into a toy SMT solver. 4 The SA T par t is handled by Minisa t [7]. Implementing a full SMT so lver for LRA was useful 4 Benc hmarks and implement ation are a v ailable from http://w ww- verimag .imag.fr/ ~ monniaux /simplexe/ . 9 for testing the so ft ware, aga inst the SMT-LIB examples . 5 W e how ever did not use SMT- LIB for b enchmarking: the p erforma nce o f a complete SMT solver is influenced b y ma n y factors, including the p erfor ma nce of the SA T s olver, the ease of adding new clauses on the fly , e tc., outside of the pure s pee d of the decision pro cedure. The flo ating-p oint simplex used is the dual s implex implemented by option GLP_DU ALP in Glpk [9], with default parameter s. The ratio nal simplex tableau is implemen ted using spar se vectors. 6 Rational num ber s ar e implemented a s quotients of tw o 32-bit nu mbers ; in case of overflo w, ex tended precision rationals from the GMP library [10] a re used. 7 The rea son behind a tw o-tier system is that GMP ra tionals inherently inv olve some inefficiencies, including dynamic memory allo cation, ev en for small nu mera tor and denomina to r. In man y exa mples arising from v er ific a tion problems, one never needs to call GMP . The simplex algo rithm, including the pivot strategy , is implemented straight from [6, 5]. It is therefor e likely that our system ca n b e implement ed as a prepr o cessor into any of the current SMT pla tforms for linear real ar ithmetic, o r even those for linear integer arithmetic, since these a re based on relaxa tions to ratio nals with additional constraints (branch-and-bound o r Gomory cuts). W e b enchm ar ked four to ols: • Our “ naive” implementation o f ra tional simplex. • The same, but with the floa ting-p oint prepro ces sing and for c ed pivoting phase descr ibed in this a r ticle (“mixed s implex”). • Bruno Dutertre and Leona rdo de Mour a’s (SRI International) Yices 1.0.9 • Nikola j Bjørner and Leonardo de Moura’s ( Micro soft Res earch) Z3, a s presented at SMT-COMP ’08. W e used tw o g roups of b enchmarks: • Benchmarks k indly provided to us b y Leo nardo de Moura, extr a cted from SMT-LIB pr oblems. Each b enchmark is an unsatisfiable conjunction, used by Z3 as a theor y lemma. Thes e ex amples are fair ly sparse, with small co efficients, and rationa l simplex seldom needs to use ex tended precis io n rational arithmetic. Despite this, the p erforma nce o f the “mixed” imple- men tation is b etter than that of Yices 1 (Fig . 1). Since the sourc e co de of neither Yices nor Z3 are av aila ble, the reasons why Z3 p erfor ms b etter than Yices 1 and b oth perfor m muc h b etter than our own rational sim- plex implementation a re somewhat a matter of conjecture. The differences probably arise from b oth a b etter s pa rse ma trix implementation, and a better pivoting strategy . • Random, dense b enchmarks. On these, o ur mixed implementation p er- forms faster tha n all other s, including the state-of-the-art SMT solvers (Fig. 2). 5 http://g oedel.cs.uiowa .edu/smt lib/benchmarks/QF_LRA.tar.gz 6 More precisely , usi ng boost::numeri c::ublas::compressed vector fr om the Boost li- brary , av ailable at http://www.bo ost.org/ . 7 http://g mplib.org/ 10 On a few ex a mples, Glpk cras hed due to an internal err or (failed a ssertion). W e are unsur e whether this is due to a bug inside this library or our misusing it — in either cas e, this is rather unsur prising g iven the complexity of current nu meric a l pa ck ag es. It is also p ossible that the numerical phas e o utputs incor - rect results in s ome c a ses, b ecause of bugs or lo ss o f precision in floating- p oint computations. Y et, this has no imp ortance — the output of the numerical phase do es not affect the cor rection of the final res ult, but only the time that it ta kes to rea ch this result. 5 Related w ork The high co st of exact-ar ithmetic in linea r progr amming has lo ng b een recog- nized, and line a r progr amming pack ages geared tow ards op erational resea rch applications seldom feature the option to p erfor m computations in ex a ct arith- metic. In co n tra s t, most SMT so lvers (e.g. Yices , Z3 ) or computational geo m- etry pack ag e s implement exa ct arithmetic. F aure et a l. ha ve exp er imented with using commercial, floa ting-p oint SMT solvers such as CPLEX 8 inside an SMT solver [8]. Their appr o ach is different from ours in that they simply s ought to reuse the yes/no re tur n v a lue pr o duced by the inexact s olver, while we also reuse the basis structure that it pro duces. Many of the difficulties t hey repo rt — for insta nce, not being able to reuse the output of the inexact so lver for theory propag ation — dis app e ar with our system. Some of their remarks still a pply: for instance, the floating-p oint solver should allow incremental use, which means that, for instance, it should not per form LU matrix factorization every time a b ound is c hanged but only at basis initializa tion. The idea o f co mbining exact and floa ting-p oint a r ithmetic for the simplex algorithm is no t new. G¨ artner prop osed a n algo rithm where mo st o f the compu- tations inside the simplex a lg orithm a re p erformed using or dina ry floating- po in t, but some are p erformed using extended pr e c ision arithmetic (floating-p oint with an arbitrary n umber of digits in the mant issa ), and repo rts improvemen ts in pre - cision compared to the industrial solver CPLE X at mo derate costs [11]. It is how ever difficult to compare his algorithm to ours , beca us e they are gea red to- wards different kinds of problems. Our alg orithm is gear ed towards the decisio n problem, a nd is mean t to be potentially incremental, and to output b oth sa t- isfiability and unsatisfiability witnesses in exact r ational arithmetic. G¨ artner’s is g eared tow ar ds o ptimization problems in computational g eometry , and do es not pr ovide o ptimality or unsatisfiability witnesses. The work that seems clo sest to ours see ms to b e LPex [4 ], a linear pro- gramming solver that fir st obtains a floating-p oint solution, then recrea tes the solution basis inside an exact precision solver and then verifies the s olution and p ossibly do es additional iterations in exac t precisio n in o rder to repair a “wrong” basis. The origina l inten t of the a uthors w as to apply this technique to optimization, not decision problems, and there ar e in fact arguments a gainst using this technique for optimizatio n that do not apply to using it for decision. A simplex- based optimizatio n solver actually runs tw o optimizatio n phases: 8 CPLEX is a professional optimization pack age geared to wards large op erational r esearc h problem, published by ILOG. http:/ /www.ilog.com /products/cplex/ 11 1. A sear ch for a fea s ible so lution (e.g. the algorithm from Sec. 2.1), which can b e equiv alently presented as an optimization problem. A ne w v aria ble δ is added, and all inequalities P a i,j x j ≤ b are repla ced by P a i,j x j − δ ≤ b . By taking δ sufficien tly large, o ne ca n find a solution p oint of the mo dified pr oblem. Optimization iter a tions then reduce δ un til it is nonnegative. Intuitiv ely , δ measures how muc h afar we ar e from a solutio n to the origina l problem. 2. O nce a so lution p oint to the orig inal pro blem is found, the original ob jec- tive function f is optimized. There ar e ther e fo re t wo ob jective functions at work: o ne for finding a solution po in t, and the “true” ob jective function. The floating-p oint simplex pr o cess, optimizes δ , then f . Then, the rational s implex, seeking to “ r epair” the resulting basis, star ts optimizing δ a gain; then it ha s to mov e back to optimizing f . In general, changing ob jective functions in the middle of an o ptimiza tion pro cess is ba d for efficiency . How ever, since we are interested in dec ision, we optimize a single function and the ob jectio n do es not hold. This appro ach w as la ter improved to co mputing exact so lutions for a ll prob- lems in NETLIB-LP , a p opular b enchmark library fo r linear pr o gramming prob- lems [14]. One improvemen t was to r e-run the floating- p oint simplex in higher precision rather than attempt to repair the “wrong” output basis using ex- act arithmetic — th us the suggestion to use extended-precision floating-p oint arithmetic 9 with increas ing pr ecisions until the solution found can b e check ed in exact precision. This last algor ithm was implemented inside the QSopt ex solver [1]. 10 These prop osa ls a r e how ever different from o ur s in that they inv olve more mo difications to the under lying ex a ct ar ithmetic solver. F or once, they c o mpute exact precision factor izations o f the input matrices, which pro bably saves time for the large op era tional resear ch pr oblems that they consider but may have a to o high overhead for the s ma ller proble ms that arise in SMT applications. In contrast, our a lgorithm can b e a dapted as a prepross essing step to the s implex pro cedure used in existing SMT solvers with hardly any mo dification to that pro cedure. 6 Conclusion and future w ork Our work leverages the maturity and efficiency of flo a ting-p oint simplex solver (inherent efficiency o f hardware floating- po in t versus r ationals, adv a nced pricing and pivoting stra tegies...) in order to sp eed up exact decision pro cedures in ca ses where these p erfor m p o orly . The main application of SMT solvers so far has b een progr am or specifica tion verification. Such pr oblems are typically sparse, with small co efficients. On s uc h problems, r ecent SMT solvers such as Yices and Z3 t ypically perfor m well using the simplex algor ithm ov er ratio na l num ber s. P erfor mance, how ever, decreases considerably if they a r e used ov er dense pro blems, since the size of numerators and denominators inv olved c a n be c ome prohibitively large. In this cas e, running our prepro cessing phase b e fo re embarking o n costly extended precis ion simplex 9 Av ail able through GNU MP’ s mpf t yp e or the M PFR library , for instance. 10 http://w ww2.isye.gatec h.edu/ ~ wcook/qs opt/index.html 12 can sa ve s ignificant time. W e suggest th at our pr o cedure b e added to such implemen tations and activ ated, as a heur is tic, when the ra tional co efficients bec ome to o larg e. Allowing SMT so lvers to sca le be yond program a nalysis examples may prov e useful if they are used for so me other a pplications than progr am pro ofs, for in- stance, forma l mathematica l pro ofs. As an example of the use of for mal a rith- metic pro ofs outside of pr ogram verification, Ha le s prov ed Kepler’s conjecture using many lemmas obtained by optimiza tio n techniques [13], but mathemati- cians o b jected tha t it was unclear whether these le mmas truly held. As a result, Hales launched a pro ject to formally prov e his theorem, including all le mma s obtained using n umerical optimization. He prop os ed transfor mations o f his or ig- inal linear progra mming pr oblems in to problems for which it is pos sible to turn the b ounds obtained by n umerical techniques in to numerical b ounds [12]. This suggests that there ar e p otential mathematical applications of efficient decision pro cedures for linear arithmetic. The applica tio ns of our improv ement s are no t limited to the linea r theor y of the rea ls or rationals . They also a pply to the linear theory of integers, or the mixed linear theory of rationals/ r eals and in teger s. In most SMT solvers, decision for in teger and mixed linear problems is implemented by relaxing the problem to the real cas e. If there is no real solution, then there is no integer so- lution; if there is a solution where the v ariables that are supp o sed to b e in teger s are integers, then the in teger o r mix e d problem has a solution. O therwise, the search is r estarted in parts of the state space, delimited by “branch and b ound” or Gomory cuts. Any efficiency impr ov ement in the ra tional simplex ca n thus translate into an improvemen t to the integer or mixed linear decision pro cedure. The rea son why we s till hav e to p er fo rm exp ensive r ational pivots even after computing a floating-p oint solution is that the floating-p oint solutions pr o duced by the simplex algo rithm a lmost a lways lie o utside of the solution po ly hedron when tes ted ov er exa ct arithmetic, as ex pla ined in Sec. 3. W e therefore think of inv e stigating interior p oint metho ds, searching bo th for a satisfiability witness and for an unsatisfiability witness. References [1] David Appleg a te, Willia m Co ok, Sanjeeb Dash, and Daniel Espinoza. Exact solutions to linear progr amming problems. Op er. R es. L ett . , 35 (6):693– 6 99, 2007. [2] Geor g e Dantzig. Line ar Pr o gr amming and Extensions . Pr inceton University Press, 19 98. [3] Leona rdo de Moura and Nik ola j Bjørner. Z3: An efficien t s m t solver. In T ACAS , pages 3 37–34 0, 2008. [4] Marc el Dhiflaoui, Stefan F unke, Car sten Kwappik, Kurt Mehlhorn, Michael Seel, Elma r Sch¨ omer, Ra lph Sch ulte, and Dennis W eb er. Cer tifying and repairing so lutions to la rge LP s: how go o d ar e LP-so lvers? In SO DA ’03: Pr o c e e dings of the fourte enth annual ACM-SIAM symp osium on Discr ete algorithms , pag es 255–2 56. SIAM, 200 3. 13 [5] Bruno Dutertre and Leonar do de Mour a. A fas t linear- arithmetic so lver for DPLL(T). In Thomas Ball and Rob ert B. Jones, editors, Computer-aide d verific ation (CA V) , volume 4144 of LNCS , pages 81–9 4. Springer, 20 06. [6] Bruno Dutertre and Leo nardo de Moura . Int egr ating simplex with DPLL(T). T echnical Rep or t SRI-CSL- 06-01 , SRI In ternatio nal, May 200 6. [7] Niklas E´ en and Niklas S¨ orensso n. An ex tensible SA T-solver. In The ory and Applic ations of Satisfiabili ty T esting (SA T ’03) , volume 2919 o f LN CS , pages 333 –336. Springer, 2 004. [8] Germain F aure, Rob e rt Nieu wenh uis, Albert Oliveras, and Enric Ro dr ´ ıguez-Carb onell. SA T mo dulo the theory o f linear a rithmetic: Ex- act, inexact and commercial solvers. In 11th Int. Conf. on T he ory and Applic ations of Satisfiability T esting (SA T) , volume 499 6 of LNCS , pag es 77–90 . Spr inger, May 2008. [9] F r ee Soft ware F oundation. GNU Line ar Pr o gr amming Kit R efer enc e Man- ual, version 4.34 , December 20 08. [10] F ree Softw a re F oundation. GNU MP The GN U Multiple Pr e cision Arith- metic Libr ary , 4.2.4 edition, September 200 8. [11] Ber nd G¨ a rtner. Exact ar ithmetic at low cost — a case study in linear progra mming. In S O DA ’98: Pr o c e e dings of t he ninth annual ACM-SIAM symp osium on D iscr ete algorithms , pa ges 157– 166. SIAM, 19 98. [12] Thoma s C. Hales. Some algorithms ar ising in the pro of of the Kepler conjecture. Av aila ble as arXiv:math/02 05209 v1. [13] Thoma s C. Hales. A pro o f of the Ke ple r co njectur e. Ann. Math. , 1 62:106 5– 185, 20 0 5. [14] Thor sten Ko ch . The final NETLIB-LP results. Op. Res. L et t ers , 32(2 ):138– 142, Mar ch 2004. [15] Alexa nder Schrijv er. The ory of Line ar and Int e ger Pr o gr amming . Wiley , April 199 8. 14