Algorithmically generating new algebraic
features of polynomial systems
for machine learning

There are difficulties in applying standard ML techniques to problems in NRA. One is the lack of sufficiently large datasets, which is addressed only partially by the SMT-LIB. The experiment in [26] found that the [QF]NRA sections of the SMT-LIB too uniform, and had to resort to random generated examples (although the state of benchmarking in computer algebra is far worse [22]). There have been improvements since then, with the benchmarks increasing both in number and diversity of underlying application. For example, there are now problems arising from biology [4], [23] and economics [37], [38].

Another difficulty is the identification of suitable features from the input with which to train the ML models. There are some obvious candidates concerning the size and degrees of polynomials, and the distribution of variables. However, this provides a starting set (i.e. before any feature selection takes place) that is small in comparison to other machine learning applications. The main focus of this paper is to introduce a method to automatically (and cheaply) generate further features for ML from polynomial systems.

Our main contributions are the new feature generation approach described in Section 3 and the validation of its use in the experiments described in Sections 4$-$5. The experiments are for the choice of variable ordering for cylindrical algebraic decomposition, a topic whose background we first present in Section 2, but we emphasise that the techniques may be applicable more broadly.

A Cylindrical Algebraic Decomposition (CAD) is a decomposition of ordered $\mathbb{R}^{n}$ space into cells arranged cylindrically: the projections of any pair of cells with respect to the variable ordering are either equal or disjoint. The projections form an induced CAD of the lower dimensional space. The cells are (semi)-algebraic meaning each can be described with a finite sequence of polynomial constraints.

A CAD is produced to be truth-invariant for a logical formula (so the formula is either true or false on each cell). Such a decomposition can then be used to perform Quantifier Elimination (QE) over the reals, i.e. given a quantified Tarski formula find an equivalent quantifier free formula over the reals. For example, QE would transform $\exists x,ax^{2}+bx+c=0\land a\neq 0$ to the equivalent unquantified statement $b^{2}-4ac\geq 0$. A CAD over the $(x,a,b,c)$-space could be used to ascertain this, so long as the variable ordering ensured that there was an induced CAD of $(a,b,c)$-space. We test one sample point per cell and construct a quantifier free formula from the relevant semi-algebraic cell descriptions.

CAD was introduced by Collins in 1975 [15] and works relative to a set of polynomials. Collins’ CAD produces a decomposition so that each polynomial has constant sign on each cell (thus truth-invariant for any formula built with those polynomials). The algorithm first projects the polynomials into smaller and smaller dimensions; and then uses these to lift $-$ to incrementally build decompositions of larger and larger spaces according to the polynomials at that level. For further details on CAD see for example the collection [12].

QE has numerous applications throughout science and engineering [41]. Our work also speeds up independent applications of CAD, such as reasoning with multi-valued functions [17] or motion planning [44].

The definition of cylindricity and both stages of the algorithm are relative to an ordering of the variables. For example, given polynomials in variables ordered as $x_{n}\succ x_{n-1}\succ\dots,\succ x_{2}\succ x_{1}$ we first project away $x_{n}$ and so on until we are left with polynomials univariate in $x_{1}$. We then start lifting by decomposing the $x_{1}-$axis, and then the $(x_{1},x_{2})-$plane and so so on. The cylindricity condition refers to projections of cells in $\mathbb{R}^{n}$ onto a space $(x_{1},\dots,x_{m}$) where $m<n$. There have been numerous advances to CAD since its inception: new projection schemes [34], [36]; partial construction [16], [43]; symbolic-numeric lifting [40], [29]; adapting to the Boolean structure [5], [20]; and adaptations for SMT [30], [9]. However, in all cases, the need for a fixed variable ordering remains.

Depending on the application, the variable ordering may be determined, constrained, or free. QE, requires that quantified variables are eliminated first and that variables are eliminated in the order in which they are quantified. However, variables in blocks of the same quantifier (and the free variables) can be swapped, so there is partial freedom. Of course, in the SMT context there is only a single existentially quantified block and so there is a free choice of ordering. So the discriminant in the example above could have been found with any CAD which eliminates $x$ first. A CAD for the quadratic polynomial under ordering $a\prec b\prec c$ has only 27 cells, but needs 115 for the reverse ordering.

Since we can switch the order of quantified variables in a statement when the quantifier is the same, we also have some choice on the ordering of quantified variables. For example, a QE problem of the form $\exists x\exists y\forall a\,\phi(x,y,a)$ could be solved by a CAD under either ordering $x\succ y\succ a$ or ordering $y\succ x\succ a$.

The choice of variable ordering can have a great effect on the time and memory use of CAD, and the number of cells in the output. Further, Brown and Davenport presented a class of problems in which one variable ordering gave output of double exponential complexity in the number of variables and another output of a constant size [10].

Heuristics have been developed to choose a variable ordering, with Dolzmann et al. [18] giving the best known study. After analysing a variety of metrics they proposed a heuristic, sotd, which constructs the full set of projection polynomials for each permitted ordering and selects the ordering whose corresponding set has the lowest sum of total degrees for each of the monomials in each of the polynomials. The second author demonstrated examples for which that heuristic could be misled in [6]; and then later showed that tailoring to an implementation could improve performance [21]. These heuristics all involved potentially costly projection operations on the input polynomials.

In [28] the second author of the present paper collaborated to use a support vector machine to choose which of three human made heuristics to believe when picking the variable ordering, based only on simple features of the input polynomials. The experiments identified substantial subclasses on which each of the three heuristics made the best decision, and demonstrated that the machine learned choice did significantly better than any one heuristic overall. This work was picked up again in [24] by the present authors, where ML was used to predict directly the variable ordering for CAD, leading to the shortest computing time, with experiments conducted for four different ML models.

Both [28] and [24] used a set of $11$ human identified features. These did lead to good performance of the models, with ML outperforming the prior human created heuristics, but a starting set of 11 features is relatively small for ML and so we hypothesise that identifying more would improve the results.

An early heuristic for the choice of CAD variable ordering is that of Brown [8], which chooses a variable ordering according to the following criteria, starting with the first and breaking ties with successive ones.

1. (1)

   Eliminate a variable first if it appears with the lowest overall degree in the input.

2. (2)

   For each variable calculate the maximum total degree for the set of terms in the input in which it occurs. Eliminate first the variable for which this is lowest.

3. (3)

   Eliminate a variable first if there is a smaller number of terms in the input which contain the variable.

Despite being computationally cheaper than the sotd heuristic (because the latter performs projections before measuring degrees) experiments in [28] suggested this simpler measure actually performs slightly better, although the key message from those experiments is that there were substantial subsets of problems for which each heuristic made a better choice than the others.

The Brown heuristic inspired almost all the features used by the authors of [28], [24] to perform ML for CAD variable ordering, with the full set of 11 features listed in Table 1 (column 3 will be explained later).

Our new feature generation procedure is based on the observation that all the measurements taken by the Brown heauristic, and all those features used in [28], [24] can be formalised mathematically using a small number of functions. For simplicity, the following discussion will be restricted to polynomials of $3$ variables as these were used in the following experiments, but everything generalises in an obvious way to $n$ variables.

Let a problem instance $\bm{Pr}$ be defined by a set of $P$ polynomials

This is the case for producing a sign-invariant CAD. Of course, any problem instance consisting of a logical formula whose atoms are polynomial sign conditions can also have such a set extracted.

In the following we define the notation for polynomial variables and coefficients that will be used throughout the manuscript. Each polynomial with index $p$, for $p=1,\dots,P$ contains a different number of monomials, which will be labelled with index $m$, where $m=1,\dots,M_{p}$ and $M_{p}$ denotes the number of monomials in polynomial $p$. We note that these are just labels and are not setting an ordering themselves. The degrees corresponding to each of the variables $x_{1},x_{2},x_{3}$ are a function of $m$ and $p$. These need to be explicitly labelled in order to allow a rigorous definition of our proposed procedure of feature generation.

We next write each polynomial as

Here, $x_{v}$ represents the polynomial variables ($v=1,2,3$). Thus for each monomial in each polynomial there is a tuple $(m,p)$ of positive integers that label it. Then in turn we denote by $d_{v}^{m,p}$ the degree of variable $x_{v}$ in that monomial, and by $c^{m,p}$ the constant coefficient, i.e., tuple superscripts are giving a label for a monomial in a problem. The original indices are simply a labelling and not an ordering of the variables $x_{1},x_{2},x_{3}.$

Therefore, any one of our problem instances $\bm{Pr}$ is uniquely represented by a set of sets

Observe now that each of Brown’s measures can be formalised as a vector of features for choosing a variable as follows.

1. (1)

   Overall degree in the input of a variable: $\max_{m,p}d_{v}^{m,p}$.

2. (2)

   Maximum total degree of those terms in the input in which a variable occurs:
   $\max_{m,p}\text{sgn}(d_{v}^{m,p})\cdot(d_{1}^{m,p}+d_{2}^{m,p}+d_{3}^{m,p})$

3. (3)

   Number of terms in the input which contain the variable: $\sum_{m,p}\text{sgn}(d_{v}^{m,p})$

In the latter two we use the sign function to discriminate between monomials which contain a variable (sign of degree is positive) and those which do not (sign of degree is zero). Of course the sign of the degree is never negative.

Define now also the averaging functions

Then the features in Table 1 can be formalised similarly to Brown’s metrics, as shown in the third column of Table 1.

We can place all of these scalars into a single framework:

where

and $g_{1},g_{2},g_{3}$, and $g_{4}$ are all taken from the set

In the above set $\max_{0}$, $\sum_{0}$, $\text{av}_{0}$ and $\text{sgn}_{0}$ are all equal to the identity function.

For example, let $\bm{Pr}=\{x_{1}^{2}x_{2}-x_{3},x_{1}x_{2}^{4}x_{3}^{2}+x_{1}x_{3}\}$. If $m=1,p=2$, then $h^{1,2}\left(\bm{Pr}\right)\in\big{\{}d_{v}^{1,2},\,\text{sgn}\left(d_{v}^{1,2}\right)\cdot\left(\textstyle\sum_{v^{\prime}}d_{v^{\prime}}^{1,2}\right)\,|\,v=1,2,3\big{\}}=\big{\{}1,1\cdot 7,4,4\cdot 7,2,2\cdot 7\big{\}}$.

We will thus consider deriving all of the other features which fall into this framework, but to do so we must first impose a number of rules.

1. 1.

   The functions $g_{1},g_{2},g_{3},g_{4}$ must all belong to distinct categories of function, i.e. one each of $\max$, $\sum$, av, and sgn.

2. 2.

   Exactly one of the functions $g_{1},g_{2},g_{3},g_{4}$ is computed over $p$ and exactly one is computed over $m$ (it may be the same one).

3. 3.

   The computation over $p$ is always performed by a function $g_{i}$ with an index $i$ greater or equal to that of the function computing over $m$.

The expression of $f\left(\bm{Pr}\right)$ can be interpreted as follows. The values $h^{m,p}\left(\bm{Pr}\right)$ are functions of variables $m$ and $p$. Each of the functions $g_{1},g_{2},g_{3},g_{4}$ either leave the function unchanged, or they turn it into a function of fewer variables (first into a function of $p$, and then into a scalar value, representing the ML feature).

The rules above are justified as follows. Rule $1$ reduces the redundancy in the feature set. Rules $2$ and $3$ guarantee that the feature $f_{v}\left(\bm{Pr}\right)$ is well defined and is a scalar number. In particular, Rule $3$ is necessary because the computation over the terms in a polynomial is dependent on their number, which is not the same for all polynomials.

The final set $\{f^{(1)}(\bm{Pr}),\dots,f^{(N_{f})}(\bm{Pr})\}$ has size $N_{f}=1728$ for a problem with $3$ variables. This number is attained as follows: we have 12 possible distributions of indexes to the functions $g_{1},\dots,g_{4}$ as shown in Table 2; then $4!$ possible orderings of those functions; and 6 possible choices for $h$. $4!\cdot 6\cdot 12=1728$.

However, many of these features will be identical (e.g. to a different placement of the identify function). We do not identify these manually now: the task that is trivial for a given dataset, but substantial to do in generality.

We use the nlsat dataset¹¹1Freely available from http://cs.nyu.edu/˜dejan/nonlinear/ produced to evaluate the work in [30], thus the problems are all fully existentially quantified. Although there are CAD algorithms that reduce what is being computed based on the quantifiers in the input (most notably via Partial CAD [16]), the conclusions drawn are likely to be applicable outside of the SAT context.

We use the $6117$ problems with $3$ variables from this database, so each has a choice of six different variable orderings. We extracted only the polynomials involved, and randomly divided into two datasets for training ($4612$) and testing ($1505$). Only the former was used to tune the parameters of the ML models.

We used the CAD routine CylindricalAlgebraicDecompose which is part of the RegularChains Library for Maple. This algorithm builds decompositions first of $n$-dimensional complex space before refining to a CAD of $\mathbb{R}^{n}$ [14], [13], [3]. We ran the code in Maple $2018$ but used an updated version of the RegularChains Library (http://www.regularchains.org). Training and evaluation of the ML models was done using the scikit-learn package [39] v0.20.2 for Python 2.7. The features for ML were extracted using code written in the sympy package v1.3 for Python 2.7, as was Brown’s heuristic. The sotd heuristic was implemented in Maple as part of the ProjectionCAD package [25].

CAD construction was timed in a Maple script that was called separately from Python for each CAD (to avoid Maple’s caching of results). The target variable ordering for ML was defined as the one that minimises the computing time for a given problem. All CAD function calls included a time limit. For the training dataset an initial time limit of $4$ seconds was used, doubled incrementally if all orderings timed out, until CAD completed for at least one ordering (a target variable ordering could be assigned for all problems using time limits no bigger than $64$ seconds). The problems in the testing dataset were processed with a larger time limit of $128$ seconds for all orderings (time outs set as $128$s).

When computed on a set of problems $\{\bm{Pr}_{1},\dots,\bm{Pr}_{N}\}$, some of the features $f^{(i)}$ turn out to be constant, i.e. $f^{(i)}(\bm{Pr}_{1})=f^{(i)}(\bm{Pr}_{2})=\dots=f^{(i)}(\bm{Pr}_{N}).$ Such features will have no benefit for ML and are removed. Further, other features may be repetitive, i.e. $f^{(i)}(\bm{Pr}_{n})=f^{(j)}(\bm{Pr}_{n}),\forall n=1,\dots,N.$ This repetition may represent a mathematical equality, or just be the case of the given dataset. Either way, they are merged into a single feature for the experiment. After this step, we are left with $78$ features: so while a large majority were redundant, we still have seven times those available in [28], [24].

Feature selection was performed with the training dataset to see if any features were redundant for the ML. We chose the Analysis of Variance (ANOVA) F-value to determine the importance of each feature for the classification task. Other choices we considered were unsuitable for our problem, e.g. the mutual information based selection requires very large amounts of data.

The training dataset consists of $N=6117$ problems with $3$ variables, and each problem is assigned a target ordering, or class $c=1,\dots,C$, where $C=6$. Let $\bm{Pr}_{c,n}$ denote problem number $n$ from the training dataset that is assigned class number $c$, $c=1,\dots,C$ and $n=1,\dots,N_{c}$, where $N_{c}$ denotes the number of problems that are assigned class $c.$ Thus $\sum_{c=1}^{C}N_{c}=N.$

The F-value for feature number $i$ is computed as follows [35].

where $\bar{f}_{c}^{(i)}$ is the sample mean in class $c$, and $\bar{f}^{(i)}$ the overall mean of the data:

The numerator in (5) represents the between-class variability or explained variance and the denominator the within-class variability or unexplained variance.

Of the 78 features the three with the highest F-values were the following

The new features may be translated back into natural language. For example, feature $65$ is the proportion of monomials containing variable $x_{2}$, averaged across all polynomials; feature $46$ the sum of the degrees of the variables in all monomials containing variable $x_{2}$, averaged across all monomials and summed across all polynomials;and feature $76$ the maximum sum of the degrees of the variables in all monomials containing variable $x_{2}$, summed across all polynomials.

Feature selection did not suggest to remove any features (they all contributed meaningful information), so we proceed with our experiment using all 78.

Four of the most commonly used deterministic ML models were tuned on the training data (for details on the methods see e.g. the textbook [2]).

• •

  The K$-$Nearest Neighbours (KNN) classifier [2, §2.5].

• •

  The Multi-Layer Perceptron (MLP) classifier [2, §2.5].

• •

  The Decision Tree (DT) classifier [2, §14.4].

• •

  The Support Vector Machine (SVM) classifier with RBF kernel [2, §6.3].

Each model was trained using grid search $5$-fold cross-validation, i.e. the set was randomly divided into $5$ and each possible combination of $4$ parts was used to tune the model parameters, leaving the last part for fitting the hyperparameters with cross-validation, by optimising the average F-score. Grid searches were performed for an initially large range for each hyperparameter; then gradually decreased to home into optimal values. This lasted from a few seconds for simpler models like KNN to a few minutes for more complex models like MLP. The optimal hyperparameters selected during cross-validation are in Table 3.

The ML approaches were compared in terms of prediction accuracy and resulting CAD computing time against the two best known human constructed heuristics [8], [18] as discussed earlier. Unlike the ML, these can end up predicting several variable orderings (i.e. when they cannot discriminate). In practice if this were to happen the heuristic would select one randomly (or perhaps lexicographically), however that final pick is not meaningful. To accommodate this, for each problem, the prediction accuracy of such a heuristic is judged to be the the percentage of its predicted variable orderings that are also target orderings. The average of this percentage over all problems in the testing dataset represents the prediction accuracy. Similarly, the computing time for such methods was assessed as the average computing time over all predicted orderings, and it is this that is summed up for all problems in the testing dataset.

The minimum total computing time, achieved if we select an optimal ordering for every problem, is $8\,623$s. Choosing at random would take $30\,235$s, almost 4 times as much. The maximum time, if we selected the worst ordering for every problem, is $64\,534$s. The K-Nearest Neighbours model achieved the shortest time of our models and heuristics, with $9\,178$s, only $6\%$ more than the minimal possible.

Since they are not affected by the new feature framework of the present paper the findings on the human made heuristics are the same as in [24]. Of the two human-made heuristics, Brown performed the best, surprising since the sotd heuristic has access to additional information (not just the input polynomials but also their projections). Obtaining an ordering for a problem instance with sotd hence takes longer than for Brown or any ML model $-$ generating an ordering with sotd for all problems in the testing dataset took over $30$min. Using Brown we can solve all problems in 10,951s, 27% more than the minimum. While sotd is only 0.7% less accurate than Brown in identifying the best ordering, it is much slower at $11\,938$s or 38% more than the minimum. So, while Brown is not much better at identifying the best, it is much better at discarding the worst!

The results show that all ML approaches outperform the human constructed heuristics in terms of both accuracy and timings. Moreover, the results show that the new algorithm for generating features leads to a clear improvement in ML performance compared to using only a small number of human generated features in [24]. For all four modules both accuracy has increased and computation time decreased. The best achieved time was 14% above the minimum using the original 11 features but now only 6% above with the new features.

The computing time for all the methods lies between the best ($8\,623$s) and the worst ($64\,534$s). Therefore, if we scale this time to $[0,100]$ so that the shortest time corresponds to $0$ and the slowest to $100$, then the best human-made heuristic (Brown) lies at $4.16$, and the best ML method (KNN) lies at $0.99$. So using ML allows us to be $4$ times closer to the minimum possible computing time.

Figure 1 shows that the human-made heuristics result in computing times that are often significantly larger than $1\%$ of the corresponding minimum time for each problem. The ML methods, on the other hand, all result in over $1000$ problems ($\sim 75\%$ of the testing dataset) within $1\%$ of the minimum time.