Efficient Local Search for Nonlinear Real Arithmetic

Satisfiability Modulo Theories (SMT) is the problem of determining the satisfiability of a formula containing both logical operators and functions interpreted in one or more custom theories [BarrettT18]. Commonly considered theories include equality, arithmetic, bit-vectors, arrays, and strings. After about two decades of development, SMT has gained widespread applications in program verification, model checking, planning, and many other areas.

The arithmetic theories can be divided according to the type of numbers involved into integer and real theories, and according to the operations allowed into difference logic, linear, and nonlinear theories. The case of (quantifier-free) nonlinear real arithmetic (NRA) considers satisfiability of equalities and inequalities involving polynomials of degree greater than one, and where the arithmetic variables take on real values. It has applications in the analysis of nonlinear hybrid automata [CimattiMT12], generating ranking functions for termination analysis [HeizmannHLP13, LeikeH15], constraint answer set programming [SusmanL16, ShenL18], and even analysis of biological networks [AkutsuHT08]. Problem instances from many of these applications are collected in the SMT-LIB benchmarks [BarFT-SMTLIB].

Similar to the SAT case, methods for solving SMT problems can be roughly divided into complete and incomplete methods. Complete methods are usually based on DPLL($T$) or close variants. They are able to both find solutions and prove unsatisfiability. Incomplete methods, such as those based on local search [HoosS2004], explores the solution space heuristically, usually by changing the assignment of one variable at a time, in an attempt to find a satisfying solution. Local search methods are not able to prove unsatisfiability, but can have an advantage over complete methods on some satisfiable instances.

Complete methods for the theory of nonlinear real arithmetic make crucial use of cylindrical algebraic decomposition (CAD) [Caviness2004QuantifierEA], which permits deciding the satisfiability of a conjunction of polynomial (in)equalities over real numbers. This can then serve as the theory solver in DPLL($T$) [NieuwenhuisOT06]. An innovation over DPLL($T$) for arithmetic theories is the MCSAT algorithm [JovanovicM12, MouraJ13], which constructs models involving both boolean and arithmetic variables at the same time. A nice overview of DPLL($T$) and MCSAT for nonlinear real arithmetic can be found in Kremer’s thesis [Kremer20].

Exploration of applying local search to solve SMT problems over arithmetic theories began only recently. The work [CaiLZ22] applied local search to the theory of linear integer arithmetic. It introduced the concept of critical moves, a change in one arithmetic variable that satisfies a previously unsatisfied clause. The algorithm iteratively applies critical moves that most improves the score, a weighted count of unsatisfied clauses. The presence of both boolean and arithmetic variables are dealt with by alternately working in the integer mode and the boolean mode, when assignments to only integer or boolean variables are changed respectively. Switching between modes are performed after the number of non-improving steps reaches a certain threshold.

Compared to linear integer arithmetic, the problem of nonlinear real arithmetic held additional challenges for local search methods. Unlike integer theories, there is an uncountable number of possible assignments to choose from, including infinite number of choices in any finite interval. Unlike the linear case, there is potential need to solve polynomials of degree greater than one, which is both costly in time and may result in variable assignments that are irrational (e.g. algebraic) numbers, causing a slow-down of the ensuing search process. It is also possible that a nonlinear constraint cannot be satisfied by changing the value of one variable alone, resulting in the lack of critical moves, so that other heuristics are needed in such scenarios.

There has already been some work exploring local search for nonlinear real arithmetic. However, largely due to the challenges listed above, none of the existing work covers the entirety of the theory. The work [abs-2303-06676] considers the multilinear case, where each variable appears with degree at most one in each polynomial constraint. The work [LiXZ23] considers problems where all equality constraints contain at least one variable that is linear. In both of these works, the problem of variable assignments to irrational values is avoided by either assuming linearity in each variable, or by limiting higher-degree constraints to strict inequalities only.

In this section, we formally define SMT problems over nonlinear real arithmetic, followed by a review of existing local search algorithms for arithmetic theories.

The syntax of a general SMT formula over nonlinear real arithmetic is as follows:

Here $x$ is an arithmetic variable, $c$ is a (rational) constant, $b$ is a boolean variable. A literal is either an atom or its negation. A clause is a disjunction of literals. In practice, we assume that input problem instances are given in conjunctive normal form (CNF), that is as a collection of clauses to be satisfied. Note that strict inequalities $p\neq 0$, $p<0$ and $p>0$ can be represented as $\neg(p=0)$, $\neg(p\geq 0)$ and $\neg(p\leq 0)$, respectively. We allow problem instances to contain boolean and arithmetic variables at the same time. We define boolean literal and arithmetic literal to mean literal whose atom is a boolean variable and a polynomial inequality, respectively.

A polynomial $p$ is linear in some variable $x$ if all terms of the polynomial have degree at most one in $x$. Alternatively, $p$ can be written in the form $p_{1}\cdot x+p_{2}$, where $p_{1},p_{2}$ do not contain $x$. A polynomial is multilinear if it is linear in each of its variables.

Given a problem instance containing boolean variables $b_{i}\leavevmode\nobreak\ (1\leq i\leq m)$ and arithmetic variables $x_{j}\leavevmode\nobreak\ (1\leq j\leq n)$, a complete assignment is a mapping from each $b_{i}$ to $\{\top,\bot\}$ and each $x_{j}$ to $\mathbb{R}$. We will only deal with complete assignments in this paper, and hence sometimes call it assignment for short. A formula is satisfied under an assignment if it evaluates to true under the standard interpretation of boolean and arithmetic operators. An assignment is a solution to a problem instance if it satisfies all its clauses. The SMT problem for nonlinear real arithmetic is to determine whether a given problem instance is satisfiable by some assignment.

Local search algorithms attempt to determine satisfiability of a problem instance by searching in the space of complete assignments, usually by changing the value of one variable at a time. In the SAT case, each move in local search flips the assignment of one boolean variable. Which flip to make is determined by factors such as the number of clauses that become satisfied/unsatisfied as result of the flip, weight of the clauses, which variables are flipped recently, and so on. For SMT over arithmetic theories, the analogous move changes the value of one arithmetic variable, usually in order to make some clause become satisfied. Such moves, called critical moves, are introduced in [CaiLZ22] for linear integer arithmetic.

For nonlinear real arithmetic, the basic procedure used to determine possible moves is root isolation for polynomials. Given a polynomial $p$ in a single variable $x$ (where here we allow the coefficients of $p$ to be algebraic numbers), the procedure computes the roots of the equation $p(x)=0$, together with the sign of the polynomial in each interval separated by the roots. Algebraic numbers in the coefficients of $p$ and in the output of root isolation are represented by their minimal polynomials, together with intervals with rational endpoints that bracket a root of that polynomial.

Given a complete assignment, and an arithmetic literal $l$ involving variable $x$, we can use root isolation to compute the set of values that assignment of $x$ may be moved to in order for $l$ to be satisfied. This is done by substituting in assignments to other variables in $l$, resulting in a polynomial containing only variable $x$, then perform root isolation and check the value of $l$ in each resulting interval. The answer is given in terms of a set of intervals (which may contain $\pm\infty$ as one of the endpoints, and may be either open or closed at each endpoint). We state the definitions precisely as follows.

A critical move is defined to be a change in the assignment of some variable $x$ to a value that satisfies some previously unsatisfied clause. The basic local search algorithm then performs critical moves at each iteration, using various scoring metrics to determine which move is chosen next. Such moves can also be interpreted as jumping between CAD cells as in [LiXZ23]. An innovation in [LiXZ23] is that when no critical moves are available for some literal, moves that change multiple variables at once along some straight line are also explored.

Scoring of critical moves is usually based on a weighted count of unsatisfied clauses. Adaptive weighting schemes assign a weight to each clause, reflecting its importance during the current search. Existing work on local search for arithmetic theories mostly use a probabilistic version of the PAWS weighting scheme [ThorntonPBF04]. This scheme is parameterized by a smoothing probability sp. Whenever there is no moves available that improves the score, with probability $1-\mathit{sp}$ the weight of each unsatisfied clause is increased by 1, and with probability $\mathit{sp}$ the weight of each satisfied clause with weight greater than 1 is decreased by 1. Then, the make-break score of each critical move equals the total weight of clauses that become satisfied by the move, minus the total weight of clauses that become unsatisfied by the move. A move is improving if its score is greater than zero.

One key contribution of [abs-2303-06676] is the introduction of make-break intervals. The idea is that instead of considering only the (in)feasible intervals of a variable $x$ with respect to some clause, we combine the (in)feasible interval information of $x$ with respect to all clauses. This results in a partition of the real line into intervals, with each interval associated to the make-break score for moving the value of $x$ into that interval. We illustrate this idea with the following example.

If boolean variables are present, the make-break score of flipping each boolean variable is defined in a similar way, as the total weight of clauses that become satisfied by the flip, minus the total weight of clauses that become unsatisfied.

Algorithm 1 shows the structure of the basic local search procedure. Begin by initializing the assignments to all variables (line 1). At each iteration, first try to find a move with the largest make-break score. If the score is greater than 0 (the variable necessarily comes from an unsatisfied clause), then perform the corresponding move (line 9). If no move has score greater than 0, it indicates that we reached a local minimum. Update the clause weights according to the PAWS scheme (line 11), and then try to make a move that makes a randomly chosen clause satisfied (line 16). If that is also not possible after several tries, randomly change the assignment of some variable in some unsatisfied clause according to some heuristic (line 19). This continues until all clauses are satisfied (line 4) or the time or step limit is reached (line 6).

There are possible variations in the choice between boolean and arithmetic variables on line 7. In [CaiLZ22, abs-2303-06676], the search is separated into modes where only boolean or arithmetic variables are considered. Alternatively, we can combine the lists of moves and decide between them based purely on make-break scores. We take the latter approach in this paper. There are other aspects of the algorithm that are left unspecified, including how the best move is computed on line 7, and the heuristic choice of moves on line 19. These will be specified in more detail in the next sections.

One key step in Algorithm 1 is computing the move with the best make-break score. The computation for boolean variables is standard (and is in any case not the bottleneck here), hence we focus on critical moves for arithmetic variables. The default approach is to loop over all variables in all unsatisfied clauses. For each variable, compute its score with respect to each clause and then combine the results (as demonstrated in Example 2.2). However, it is clear that this may result in repeated computations across iterations. For example, the feasible set of some variable with respect to a clause may be recomputed, even though none of the variables in that clause are changed in the previous step. Following the idea of caching and updating scores in GSAT [SelmanLM92], we propose data structures for caching and updating score information for arithmetic variables.

We define a boundary to be a quadruple $\langle\mathit{val},\mathit{is\_open},\mathit{is\_make},\mathit{cid}\rangle$, where $\mathit{val}$ is a real number, $\mathit{is\_open}$ and $\mathit{is\_make}$ are boolean values, and $\mathit{cid}$ is a clause identifier. It indicates that there is a change in make-break score when moving from less than to greater than $\mathit{val}$ due to clause $\mathit{cid}$. If $\mathit{is\_make}$ is $\top$, the score is increased by the weight of clause $\mathit{cid}$, otherwise it is decreased by the weight. If $\mathit{is\_open}$ is $\top$, the change is not active at $\mathit{val}$, otherwise it is already active at $\mathit{val}$. There is a natural ordering among boundaries, first order by $\mathit{val}$ and then by $\mathit{is\_open}$ (with $\bot<\top$). The make-break score information of each variable with respect to each clause can be characterized by a starting score (indicating the make-break score of large negative values), together with a set of boundaries. The make-break information of each variable with respect to all clauses is formed by summing the starting score and taking the union of the sets of boundaries. We illustrate the computations in the following example.

During local search, after each move of variable $v$ to a new value, only those variables $v^{\prime}$ that shares a clause with $v$ need to have their make-break score information updated (this is analogous to the concept of dependent variables in the SAT case), and moreover boundary information need to be updated for the shared clauses only. This is summarized in Algorithm 2. The set $S$ collects the set of variables sharing a clause with $v$. Line 5 recomputes starting score and boundary information. Line 7 recomputes best critical move and score for each updated variable.

Equality constraints with degree greater than one pose special difficulty for local search, since it may force assignments of variables to irrational (e.g. algebraic) numbers. For example, for the constraint $x^{2}+y^{2}+z^{2}=1$, with most rational assignments to $x$ and $y$, the assignment to $z$ would be forced to irrational in order for the constraint to be satisfied. While it is possible to represent and compute with algebraic numbers during local search, the time-cost of such computation is significantly increased. Even without considering algebraic numbers, numbers with increasingly large denominators are also problematic for slowing down the search process.

Both [abs-2303-06676] and [LiXZ23] avoid algebraic numbers by either limiting to the multilinear case, or considering only equality constraints with at least one linear variable (and solving only for those linear variables). The work [abs-2303-06676] further incorporated comparison of size of denominators (as well as absolute value) of potential assignments into the scoring heuristic, in order to keep the complexity of assigned values as low as possible.

We propose a novel approach to address the problem of assignments to irrational values and values with large denominators, that allow the algorithm to be applied efficiently to the full set of nonlinear arithmetic problem instances. The approach still relies on comparing complexity of assigned values, hence we first define it below.

The relaxation mechanism can be described simply as follows: whenever some equality (or inequality) constraints force an assignment of some variable to a comparatively complex value, those constraints are relaxed before continuing the local search process, so that such assignments never actually occur. The implementation is parameterized by two thresholds. The parameter $\epsilon_{v}$ (for variable threshold) specifies the complexity of assigned values (according to Definition 4.1) beyond which relaxation of constraints should be applied. The parameter $\epsilon_{p}$ (for polynomial threshold) specifies the amount of relaxation of polynomial constraints. Both $\epsilon_{v}$ and $\epsilon_{p}$ are chosen to be $10^{-4}$ in the implementation.

It should be noted that the constraints $p\geq 0$ and $p\leq 0$ together can also force the assignment of variables to irrational values. These constraints may appear as part of clauses with more than one literal, and hence are not equivalent to $p=0$. This means in general we consider relaxation of non-strict inequalities, although we will still use the slightly imprecise (but more intuitive) description of relaxing equalities throughout the paper.

The detailed method for determining which constraints to relax is as follows. When computing the best make-break score of a variable $v$, if that score comes from a one-point interval, the set of clause identifiers in the boundaries contributing to that interval are recorded. If the variable $v$ is chosen, and the new value of $v$ is more complex than both $\epsilon_{v}$ and any other previously assigned value (according to Definition 4.1), all equalities and non-strict inequalities in the recorded clauses that contribute to the boundary are relaxed. The result of relaxation is as follows.

• •

  If the constraint is of the form $p=0$, it is relaxed into the pair of inequalities $p<\epsilon_{p}$ and $p>-\epsilon_{p}$.

• •

  If the constraint is of the form $p\geq 0$, it is relaxed into $p>-\epsilon_{p}$. Likewise, if the constraint is of the form $p\leq 0$, it is relaxed into $p<\epsilon_{p}$.

Note that strict inequality constraints cannot force a variable to a particular value. After relaxation, the local search process proceeds as before, but with all evaluation of literals and computation of make-break scores according to the relaxed interpretation of literals.

After local search finds a “solution” under the relaxation of some constraints, it is only an approximate solution. In fact, there is no guarantee that there is an exact solution nearby. Instead, we use local search itself to try to move the approximate solution to an exact solution. First, restore the relaxed constraints to their original forms, and then proceed with local search, until either an exact solution is found, or no improvement is made for a certain number of steps. In the latter case, the search continues with a minor restart (see Section 5.2) and with relaxation of constraints allowed as before.

The above description is summarized as a modification of the overall algorithm, shown in Algorithm 3. The main change is to add a boolean flag called $\mathit{use\_slack}$, that specifies whether relaxation of constraints is active. If all clauses are satisfied and $\mathit{use\_slack}=\top$, then we may only have an approximate solution, so we restore the relaxed constraints and continue local search (line 4-5). If all clauses are satisfied and $\mathit{use\_slack}=\bot$, then we obtained an exact solution (line 7). If no improvement is made for $T_{1}$ steps in the $\mathit{use\_slack}=\bot$ mode, we assume that no exact solution is found nearby, so we assign $\mathit{use\_slack}=\top$ and continue after minor restart (line 11-12).

In this section, we describe the implementation in more detail. First, we explain our choice of heuristic move selection when encountering a literal without critical moves. Then, we describe some further details on preprocessing, restart mechanism, and other efficiency improvements.

In this section, we compare the implementation with those of complete procedures in existing SMT solvers Z3 [MouraB08], CVC5 [BarbosaBBKLMMMN22] and Yices [Dutertre14], as well as previous work on local search for (fragments of) nonlinear arithmetic. We also perform ablation study on the two improvements described in Section 3 and 4.

The benchmark used in the evaluation comes from SMT-LIB’s QF_NRA theory. The benchmark consists mostly of industrial problems from various applications of constraint solving in nonlinear real arithmetic. Among the benchmarks, we took those problem instances whose status is either satisfiable or unknown. This results in a total of 6604 instances. The experiments are run on a cluster of machines with Intel Xeon Platinum 8153 processor at 2.00GHz. Each experiment is run with a time limit of 20 minutes (as in the SMT competition) and memory limit of 30GB.

In this paper, we presented improvements to the local search algorithm for solving SMT problems in nonlinear real arithmetic. Building upon the basic structure of local search, we presented incremental computation of variable scores and temporary relaxation of constraints. We also described heuristic move selection with look-ahead for dealing with literals without critical moves, and implementation details to improve efficiency. The resulting implementation is competitive against complete algorithms based on DPLL($T$) and MCSAT on satisfiable problem instances, as implemented in other SMT solvers. It is the first local search algorithm designed for the entirety of nonlinear real arithmetic, covering a wider range of problems than existing work [abs-2303-06676, LiXZ23].

While the methods proposed in this paper made progress in addressing challenges of local search for nonlinear real arithmetic, there are remaining problems that represent interesting directions of future work.

• •

  Look-ahead for critical moves presents another way to improve upon random search in cases when no critical move is available. On the other hand, methods based on CAD or polynomial optimization would give a more complete way for determining assignments to satisfy a certain literal. A major challenge is how to incorporate such algorithms into local search in an efficient way.

• •

  In the current work, after an approximate solution is found that satisfies the relaxed version of equalities, we use local search itself to attempt to find an exact solution nearby. While there may be more straightforward ways to find exact solutions for limited cases, e.g. when each relaxed clause has a unique variable that can be adjusted, the situations we found in actual benchmarks are usually more complicated. Designing a more general algorithm for finding exact solutions near approximate solutions (or determining that none exist) is an interesting problem that we leave to future work.

• •

  Finally, local search can be incorporated into complete methods such as DPLL, improving its performance even on unsatisfiable instances, as shown by the works [CaiZ21, CaiZFB22]. It is interesting to investigate this possibility for SMT problems over nonlinear real arithmetic. The improvements in efficiency in our work would be very helpful for such combination, as in those cases local search is only given very short running times.