1 Introduction
Tropical geometry [brugall2013bit] has been used to find the order of time scales of variables in chemical reaction networks [SamalGrigorievFroehlichWeberRadulescu2015] and for model reduction. It has applications in economics and optimizations like network flows and scheduling.
Satisfiability Modulo Theories (SMT) checking [monniaux2016survey] is usually built on top of SAT (Boolean satisfiability), which is the first problem that was proved, in the form of 3SAT, to be NPcomplete. SMT checking allows one to test a logical formula with unknowns and relations for satisfiability and, if it is so, for an assignment of the unknowns that leads to the formula’s satisfiability. SMT checking is used today in verification of computer hardware and software and has advanced much in recent years due to advances in technology and industrial applications [de2011satisfiability].
We present a novel approach to use SMT checking to compute the tropical equilibrium (resp. prevariety). We believe this to be of use, since SMT is a very active field of research, yet problems of tropical geometry have not been solved with SMT until now.
In the following, we describe the idea of tropical geometry and SMT in the remainder of this section. In Sect. 2 we describe the exact problem and Sect. 3 describes the proposed solution and several possible improvements. Section 4 has results of speed tests of our implementation SMTcut for various SMT solvers and possible optimizations.
1.1 Some Tropical Geometry
Given a system of polynomial equations with zero r.h.s., the basic idea is to express its indeterminates and parameters as powers of some times a value resp. that is roughly 1. That is, indeterminates become and parameters , with .
Let and , be multiindices with . Furthermore, let and be index sets of the terms with positive resp. negative sign. We write the polynomial equation system as
(1) 
Next, we sort all terms with positive sign to one side and all terms with negative sign to the other. Hence, the equation becomes
(2) 
The critical observation is now that on each side almost all of the time one term dominates all others [Viro2000]. In tropical geometry we only look at this dominating term. A tropical root is found when the dominating terms on both sides cancel each other out. Since , the inequality is approximately equivalent to . “Domination” thus means that is minimal.
With the above sketched idea of tropicalization, we transform equation (2) into its tropical counterpart to look for tropical roots:
In order for this equation to hold, the minimum has to be attained at least twice, one time on each side. Observe that the equation consists now only of minima of linear functions and can thus readily be expressed as a set of polyhedra.
One polyhedron is defined by each combination of and
that yields a hyperplane via
(3) 
while for all halfspaces are defined by
(4) 
A set of polyhedra is defined by cycling over all possible choices for and .
Since a polynomial system may consist of multiple equations, we get multiple sets of polyhedra. Because we are looking for solutions where all polynomials are zero at the same time, all constraints (3) and (4) have to hold at the same time and hence the sets of polyhedra have to be intersected. The resulting set of polyhedra is called a tropical equilibrium.
For computations, the parameter in the above construction has to be provided. Furthermore, the are rounded to rationals. Tropical geometry focuses on the dominant terms and hence yields only approximate results. The resulting coarsegraining can be helpful to find broad areas where the dominant terms cancel each other out.
Equations like (1) can be used to model chemical reaction networks, where all quantities (concentrations and reaction speeds) are positive values. The resulting matching of positively and negatively signed terms in (2) is called opposite sign condition.
If no such distinction is made and all terms—regardless of their sign—can cancel each other out, then a tropical prevariety is computed. In the above formalism this can be achieved by defining both and to contain the indices of all terms. The tropical prevariety is a superset of the tropical equilibrium.
We will continue to speak of equilibria, but the presented algorithm works for prevarieties as well, since the difference is only in the input. In Sect. 4 we show runtimes only for computations of equilibria, since the inputs for these benchmarks come from chemical reaction networks.
1.2 Bringing SMT into the Picture
The problem that we are solving is the intersection of several unions of polyhedra. That is, given polyhedra and unions of (convex) polyhedra , we are interested in the intersection , where the are again (convex) polyhedra. In this article we show how to use SMT checking to solve this problem.
Satisfiability Modulo Theories (SMT) checking allows us to decide if a logical formula, with atoms that are themselves equations or inequalities, is satisfiable or not. For example, is such an SMT formula. One has to specify a theory of numbers that unknowns in the formula can assume. In the above example, the problem is satisfiable in the theory of real numbers, but not in the theory of integers. If an SMT problem is satisfiable, SMT can return a model, which is an assignment for all unknowns in the formula.
SMT solvers may be used in incremental mode, where one can add additional assertions, i.e., clauses that are combined with AND, and continue to look for further models after one has found a solution. This can save a lot of time and we will make use of it later.
2 The Problem
A (convex) polyhedron is defined as the intersection of finitely many hyperplanes and halfspaces. Furthermore, each hyperplane can be expressed as two (closed) halfspaces, thus a polyhedron can be described as a finite number of halfspaces [Ziegler1995].
Given the ambient space , a (closed) halfspace is a set
(5) 
Given halfspaces , we define a polyhedron as
(6) 
A bag is what we call a union of polyhedra , i.e.,
Finally, we are looking for the intersection of said bags , that is
(7) 
The naive solution to the problem of computing the intersection is to cycle successively through all combinations. To do that, pick two bags and , , and intersect all polyhedra from one with all polyhedra from the other to form a new bag . Then, remove and from the set of bags and insert instead. Continue this procedure until there is only one bag left, which will then consist of the sought polyhedra . This is the solution that was used in [SamalGrigorievFroehlichWeberRadulescu2015] and, with some refinements, in PtCut [LuedersPtCut].
The problem with this solution is that the complexity is exponential in the number of bags. In practice, it often happens that the number of intermediate results, i.e., the number of polyhedra in some , is very high, even if in the end there are only a few solution polyhedra. This intermediate expression swell makes computing the intersection infeasible for some models.
3 The Procedure
First, we have to formulate our problem as an SMT problem. Fortunately, it is easy to convert a polyhedron as defined in (6) into a logical formula. Set theory maps easily to logical formulas with union mapping to logical OR and intersection to logical AND. In the following, denotes the logical formula that defines the set . Thus, (7) expands to
and definition (5) of employs a linear function that can be used as a formula in SMT. Call the resulting SMT formula . We can use an SMT solver to decide the satisfiability of and, what’s more important, get an that satisfies the constraints, if there is one.
Next, we look for a matching polyhedron that includes and is included in the solution . Since is contained in the intersection of the , it must be contained in at least one polyhedron per bag . Thus, we cycle through all to find a containing , call it . (There may be more than one , but any will do.)
Obviously, the intersection includes , but most likely has higher dimension than that. Furthermore, since is the intersection of exactly one polyhedron per bag it is included in as well. Hence, we have found a polyhedron that includes .
In the next step, we modify our initial formula to exclude the polyhedron , like this:
Notice that we are only adding another assertion to the formula, so we can utilize the incremental mode of SMT solvers to save (a lot of!) time for its next computation.
The important observation here is that we are expanding the original formula —which describes all solution points—to exclude what we already know to be a solution and continue the search. Thus the procedure generates an ever growing subset of the solution , making it an anytime algorithm [zilberstein1996using]. We can iterate this process until formula is unsatisfiable.
This is the algorithm in Pythonstyle pseudocode:
The result of this function is a list of polyhedra. Mathematically, this union of polyhedra describes the equilibrium (resp. prevariety) . Yet, there are some problems that we address in the next section.
The logic used for SMT formulas is QF_LRA, that is, quantifierfree linear real arithmetic (here “real” means rational). This allows Boolean propositional logic of equations/inequalities consisting of linear functions over elements of [barrett2010smt].
3.1 Improvements to the Procedure
Nonmaximal Polyhedra
The main issue we experience with the procedure compute_polyhedron_dnf is that the polyhedron computed from point is often not maximal. That is, is only a lower dimensional face of a higherdimensional polyhedron. The full highdimensional polyhedron will eventually be found by the procedure, but earlierfound lowerdimensional faces would still remain in the result list rr, albeit superfluous.
To avoid this, we test if each newly found polyhedron is included in one of the already computed polyhedra of result list rr. Unfortunately, this causes quadratic runtime in the number of resulting polyhedra. But there is an observation that can reduce the constant.
If the newly found polyhedron is included in some already found polyhedron , then obviously, point is included in as well. Testing if a point is included in a polyhedron is simple and fast, so one can test this first. Only if this test succeeds one must perform the full polyhedron inclusion test. Measurements show that with this heuristic, almost all polyhedron inclusion tests can be avoided. See Sect. 4 for details.
In our procedure, we would have to modify function contains in line 22 and function append in line 23 according to these observations.
Superfluous Constraints
Another issue is the redundancy of the constraints that are collected in line 23. Efficiency can be increased by minimizing the set of constraints: the larger the number of constraints, the larger the memory demand and, of course, SMT checking times.
One can simply cycle through all constraints, test if each of them is really required and if not, drop it. The remaining set is not necessarily a minimal set, though.
Here’s how this can be done: Let be the constraint in consideration, the formula before and after the addition. If is more restrictive than (i.e., makes a difference), then the following is unsatisfiable:
We use SMT checking and apply this formula to all constraints to drop superfluous ones.
Preprocessing
We explored the possibility to improve the speed of the procedure by preprocessing the input, i.e., the sets of polyhedra.
For one, one can collect all constraints from all bags with only one polyhedron each. Call the resulting polyhedron . Because of distributivity these constraints hold for all polyhedra of the solution. Hence, we can intersect all polyhedra in their bags with to test if the intersection is empty, in which case we drop the polyhedron from its bag to reduce computation time.
A more powerful version of this technique can be used to test polyhedra in all bags on if they are required for the definition of the solution in (7). Let be the polyhedron in question, the union of all other polyhedra in ’s bag and the intersection of all other bags. Then the solution is . If is required, then and in particular . Thus, the following set is nonempty:
This can easily be tested with SMT checking for each polyhedron per bag and the superfluous polyhedra are dropped, again reducing the computation time.
4 Benchmarks
To benchmark the procedure compute_polyhedron_dnf, we created input sets of polyhedra of chemical reaction networks. The reaction networks were taken from the BioModels database [BioModels2006], which contains models formulated in SBML [Hucka2003]. The SBML models were converted to ODE systems with polynomial or rational functions using ODEparse [LuedersODEparse]. The resulting ODE systems can be downloaded from ODEbase [LuedersErramiNeidhardtSamalWeber2019] at http://odebase.cs.unibonn.de. If the ODE systems contained rational functions, we multiplied each equation with its common denominator to get a polynomial.
These polynomials were then tropicalized as sketched in Sect. 1.1. The parameter was set to and the logarithms in (3) and (4) were rounded to integers. The sets of polyhedra created by tropicalization were saved with polyhedra expressed as sets of equalities and inequalities.
The software used for polyhedral computation was PtCut 3.5.1 [LuedersPtCut] and was run under Windows 10 64bit using Python 3.7.7. Polyhedral computation were done with the help of the Parma Polyhedra Library (PPL) [BagnaraHillZaffanella2008], version 1.2.
The procedure compute_polyhedron_dnf was implemented in Python as SMTcut and SMTcut 4.6.4 was used for benchmarking. SMTcut is available under a free software license from https://gitlab.com/cxxl/smtcut.
Neither PtCut nor SMTcut make active use of multithreading. All input and output files are available from https://gitlab.com/cxxl/smtcut_data_2 in one large ( 42 MiB) repository.
The machine for tests was an Intel Core i75820K CPU with 48 GB of memory. We used the solveragnostic framework PySMT 0.9.0 [pysmt2015] to access SMT solvers. Of the seven supported solvers, the solvers CUDD, PicoSAT and Boolector could not be used, since they do not support the QF_LRA logic. This left us with the SMT solvers MathSAT 5.6.1 [mathsat5], Z3 4.8.7 [z3smt], CVC4 1.7prerelease [CVC4] and Yices 2.6.2 [yicessmt]. SMTcut was run on Ubuntu Linux 18.04 64bit with Python 3.7.6.
Table 1 shows a comparison of runtimes for all BioModels that completed within two hours, sorted by ascending runtime.
Runtime [s]  Rounds  Runtime / round [s]  

BM  MSAT  Z3  CVC4  Yices  MSAT  Z3  CVC4  Yices  MSAT  Z3  CVC4  Yices 
397  0.1  0.1  0.1  0.1  1  1  1  1  0.128  0.108  0.147  0.092 
26  0.1  0.1  0.1  0.1  7  8  8  8  0.016  0.017  0.016  0.012 
41  0.2  0.2  0.2  0.1  7  10  7  8  0.025  0.022  0.030  0.017 
30  0.2  0.2  0.2  0.2  8  8  7  8  0.026  0.028  0.029  0.020 
492  0.3  0.2  0.3  0.2  2  2  2  2  0.141  0.120  0.158  0.083 
491  0.3  0.2  0.3  0.2  2  2  2  2  0.146  0.118  0.130  0.084 
48  0.8  0.7  0.8  0.4  12  13  8  8  0.068  0.054  0.096  0.044 
482  0.5  0.8  0.6  0.4  19  25  21  20  0.028  0.032  0.027  0.020 
28  0.6  0.9  0.4  0.5  25  35  20  32  0.022  0.025  0.022  0.017 
637  0.6  0.7  0.6  0.5  21  27  19  24  0.026  0.027  0.031  0.020 
221  0.7  0.8  0.7  0.6  55  58  53  68  0.012  0.014  0.013  0.009 
21  0.7  0.9  0.7  0.6  56  65  61  67  0.012  0.015  0.011  0.009 
315  0.8  1.0  0.8  0.8  25  30  19  36  0.031  0.033  0.041  0.022 
328  2.0  2.9  1.7  1.5  93  105  89  90  0.022  0.027  0.019  0.017 
200  2.0  1.9  2.2  1.6  38  41  34  45  0.053  0.048  0.063  0.035 
599  2.1  2.3  1.6  1.7  48  48  36  52  0.043  0.048  0.045  0.033 
22  2.2  3.5  1.8  2.3  199  235  165  259  0.011  0.015  0.011  0.009 
666  2.7  3.0  2.2  2.6  93  85  85  111  0.029  0.038  0.026  0.023 
222  3.2  4.0  3.1  2.3  246  257  228  230  0.013  0.016  0.014  0.010 
638  2.6  2.9  3.2  2.4  37  46  29  51  0.071  0.062  0.109  0.047 
489  4.3  4.9  3.0  2.5  72  82  49  54  0.059  0.060  0.061  0.046 
365  4.6  5.1  4.3  3.6  80  82  74  82  0.058  0.062  0.059  0.044 
396  3.9  4.8  4.2  5.0  75  82  86  121  0.053  0.059  0.049  0.041 
147  4.1  6.0  4.4  4.0  84  125  90  117  0.049  0.048  0.049  0.035 
230  5.5  7.0  5.1  5.3  115  143  99  139  0.048  0.049  0.052  0.038 
498  6.5  9.1  6.1  6.9  231  285  242  333  0.028  0.032  0.025  0.021 
431  7.7  9.7  6.6  7.4  199  224  183  243  0.039  0.043  0.036  0.030 
105  9.6  9.4  9.7  6.7  141  153  143  152  0.068  0.061  0.068  0.044 
102  14.8  19.1  12.2  11.8  399  538  395  497  0.037  0.036  0.031  0.024 
32  22.0  23.1  12.8  23.5  418  388  282  504  0.053  0.060  0.046  0.047 
407  20.7  30.2  15.7  55.0  307  389  261  881  0.067  0.078  0.060  0.062 
477  42.0  72.0  33.8  49.3  793  940  745  1051  0.053  0.077  0.045  0.047 
576  55.1  78.0  48.6  62.9  838  1064  881  1076  0.066  0.073  0.055  0.058 
93  124.4  121.0  55.3  139.7  1390  1345  870  1865  0.090  0.090  0.064  0.075 
183  69.8  132.0  77.7  115.1  29  194  14  284  2.408  0.681  5.551  0.405 
501  116.1  192.2  102.8  352.5  1018  1516  1045  2819  0.114  0.127  0.098  0.125 
430  209.8  196.2  128.2  405.5  2321  2448  2142  4798  0.090  0.080  0.060  0.085 
103  231.0  449.2  224.9  337.2  2781  3072  2609  2923  0.083  0.146  0.086  0.115 
74  1165.6  1118.5  941.5  1694.2  14227  13352  13352  23245  0.082  0.096  0.071  0.073 
73  1087.1  1647.9  1067.7  1485.4  14977  16294  15522  21118  0.073  0.101  0.069  0.070 
61  3045.7  4084.1  1661.3  4187.1  12564  16326  12915  24078  0.242  0.250  0.129  0.174 
The lowest computation times were achieved with CVC4, if—with the exception of BioModel 183—computation took more than 12 seconds. Below that time Yices was often faster. If we look at runtime per round we see that Yices had always the lowest runtimes per round for models that took less than 12 seconds to compute. For unknown reasons, Yices often needs more rounds than other solvers, especially CVC4.
The geometric means of runtimes of completed computations for CVC4, Yices, MathSAT and Z3 were 4.93, 5.37, 5.49 and 6.77 seconds, respectively. For that reason, we used CVC4 as solver in further comparisons.
Table 2 shows a comparison of runtimes between PtCut and SMTcut with CVC4 as solver. We can make several observations from the data.
BM  R  Dim  Comb  PH  IntMaxPH  SMTcut  PtCut  Speedup 

26  10  6  32  0.1  0.1  1.0  
397  *  10  0  —  0.1  14.9  101.6  
30  16  6  48  0.2  0.3  1.7  
41  11  4  924  0.2  1.4  6.8  
491  36  1  1  0.3  17.2  66.6  
492  18  1  1  0.3  2.9  9.1  
28  22  17  119  0.4  0.4  1.0  
482  *  23  17  495  0.6  2.3  4.1  
637  16  12  1140  0.6  4.2  7.1  
221  *  19  50  2573  0.7  2.1  3.1  
21  *  33  46  3408  0.7  3.0  4.3  
48  *  10  5  1160  0.8  29.2  38.1  
315  13  13  432  0.8  2.5  3.3  
599  24  24  456  1.6  5.8  3.6  
328  17  86  140  1.7  1.4  0.8  
22  *  67  147  1170  1.8  2.0  1.1  
200  27  20  4704  2.2  64.0  29.7  
666  *  22  64  464  2.2  5.6  2.6  
489  24  42  4824  3.0  57.5  19.4  
222  *  8  192  12516  3.1  11.0  3.5  
638  8  13  2124  3.2  32.3  10.2  
396  *  18  54  972  4.2  21.7  5.2  
365  19  70  15030  4.3  583.7  134.2  
147  30  54  5069  4.4  18.7  4.2  
230  36  68  3330  5.1  24.6  4.8  
498  *  50  214  1750  6.1  3.3  0.5  
431  47  155  984  6.6  13.2  2.0  
105  27  130  21088  9.7  644.2  66.5  
102  27  322  4784  12.2  7.2  0.6  
32  43  244  1092  12.8  110.0  8.6  
407  23  212  15010  15.7  968.2  61.7  
477  35  467  23460  33.8  571.0  16.9  
576  57  756  10752  48.6  55.8  1.1  
93  52  596  47772  55.3  6113.4  110.5  
183  19  1  —  77.7  > 72000.0  > 926.4  
501  35  916  —  102.8  > 72000.0  > 700.3  
430  34  1676  4683  128.2  58.0  0.5  
103  30  1938  111402  224.9  763.5  3.4  
74  12  9685  —  941.5  > 72000.0  > 76.5  
73  21  13449  —  1067.7  > 72000.0  > 67.4  
61  34  10084  1784868  1661.3  26408.6  15.9  
14  86  > 2663  —  > 7200.0  > 72000.0  —  
151  66  > 7734  —  > 7200.0  > 72000.0  —  
410  53  > 10612  —  > 7200.0  > 72000.0  —  
560  59  > 12577  —  > 7200.0  > 72000.0  —  
730  *  45  > 15310  —  > 7200.0  > 72000.0  — 
runtime. Minimal runtimes are set in bold. Columns “BM”, “Dim”, “Comb”, “PH”, and “IntMaxPH” contain the BioModel number, dimension, theoretic number of combinations, number of polyhedra in the solution, and maximal intermediate polyhedra, respectively. A star in column “R” signifies a model with a rational vector field.

BioModels 14, 151, 410, 560 and 730: each of them could not be computed with PtCut and it is likely infeasible to do so. SMTcut was at least able to compute part of the solution, even though it is unknown how large a part of the full solution this constitutes.

BioModels 183 and 491, 492: Here SMTcut was able to play out its full potential: with only one polyhedron in the solution, it took only some rounds until the maximal polyhedron was found. On the other hand, the computation of BioModel 183 with PtCut was terminated after 20 hours of work with an intermediate number of 15000 polyhedra and still only 5 of 65 iterations done.

BioModel 397 does not have a solution at all. SMTcut find this in 0.1 s, whereas PtCut takes almost 15 s.

BioModels 73, 74, 93, 105, 183, 365, 397, 407, 491, and 501: computation time of SMTcut was at least 60 times lower than with PtCut.

BioModels 102, 328, 430 and 498: here PtCut was up to 2 times faster than SMTcut. The reason is that for models of dimension less than 20 and with many polyhedra, PtCut can be faster if the intermediate expression swell is not too large.

In this overview, SMTcut was the better choice in all cases where PtCut needed more than 58 s of computation time, and SMTcut never took more than twice the time as PtCut.
4.1 Benchmarks of Preprocessing
Table 3 shows a comparison of times with and without preprocessing the input. We make the following observations:
BM  Dim  Comb  PH  No PP  With PP  PP time  After PP  Speedup 

397  10  0  0.1  0.5  0.5  0.0  0.29  
183  19  1  77.7  2988.7  2988.7  0.0  0.03  
491  36  1  0.3  3.4  3.4  0.0  0.08  
492  18  1  0.3  3.4  3.4  0.0  0.09  
41  11  4  0.2  2.6  2.6  0.1  0.08  
48  10  5  0.8  29.8  29.6  0.2  0.03  
26  10  6  0.1  0.9  0.8  0.1  0.14  
30  16  6  0.2  2.7  2.6  0.1  0.07  
637  16  12  0.6  5.3  5.0  0.3  0.11  
315  13  13  0.8  9.3  8.9  0.4  0.08  
638  8  13  3.2  54.9  54.2  0.8  0.06  
28  22  17  0.4  2.5  2.0  0.4  0.18  
482  23  17  0.6  2.8  2.5  0.4  0.20  
200  27  20  2.2  29.2  28.3  0.9  0.07  
599  24  24  1.6  9.7  8.4  1.2  0.17  
489  24  42  3.0  59.9  57.6  2.3  0.05  
21  33  46  0.7  2.2  1.6  0.6  0.31  
221  19  50  0.7  2.9  2.3  0.6  0.24  
147  30  54  4.4  20.8  18.1  2.7  0.21  
396  18  54  4.2  20.6  18.3  2.3  0.20  
666  22  64  2.2  4.8  2.7  2.1  0.45  
230  36  68  5.1  49.6  44.4  5.2  0.10  
365  19  70  4.3  13.6  10.8  2.8  0.32  
328  17  86  1.7  3.5  1.8  1.7  0.49  
105  27  130  9.7  93.1  85.6  7.5  0.10  
22  67  147  1.8  3.7  1.8  2.0  0.48  
431  47  155  6.6  18.2  11.8  6.4  0.36  
222  8  192  3.1  6.7  3.7  3.0  0.47  
407  23  212  15.7  27.4  10.9  16.5  0.57  
498  50  214  6.1  8.6  3.2  5.4  0.71  
32  43  244  12.8  17.4  5.8  11.6  0.74  
102  27  322  12.2  25.7  15.7  10.0  0.48  
477  35  467  33.8  26.5  5.6  20.8  1.28  
93  52  596  55.3  86.2  39.0  47.2  0.64  
576  57  756  48.6  62.0  20.0  42.0  0.78  
501  35  916  102.8  87.0  13.2  73.9  1.18  
430  34  1676  128.2  119.1  7.7  111.4  1.08  
103  30  1938  224.9  589.2  422.2  167.0  0.38  
74  12  9685  941.5  861.7  4.6  857.1  1.09  
61  34  10084  1661.3  1579.5  139.5  1439.9  1.05  
73  21  13449  1067.7  1024.0  5.0  1018.9  1.04 

The computation for BioModel 183 took almost one hour to preprocess. This is due to the extraordinary number of polyhedra in some input bags.

Taking the time for preprocessing into account, models with more than 800 polyhedra as solution—with the exception of BioModel 103— gained a moderate speedup (4%–18%) through preprocessing.

Preprocessing increased the overall runtime for all models with less than 400 polyhedra as solution, sometimes by a factor of 10 or more.
4.2 Profiling of the Different Parts
The relative time needed in different parts of the procedure to compute the whole solution varies with the number of resulting polyhedra. Table 4 shows the relative time spent in three different parts of the procedure:

S: Searching for another point outside the already known polyhedra,

M: Minimizing the found polyhedron,

I: Inserting that polyhedron into the list of already known ones and testing for inclusion.
BM  PH  Time  % in S  % in M  % in I 

183  1  77.7  91.7  4.3  0.1 
638  13  3.2  63.6  27.1  1.3 
200  20  2.2  52.8  35.9  1.5 
599  24  1.6  28.0  56.2  3.0 
489  42  3.0  31.6  55.2  2.9 
147  54  4.4  34.2  53.3  2.9 
396  54  4.2  19.1  62.5  6.3 
666  64  2.2  16.9  60.2  6.8 
230  68  5.1  33.5  53.3  4.4 
365  70  4.3  31.3  55.0  3.5 
328  86  1.7  16.9  56.8  5.7 
105  130  9.7  43.7  45.8  3.1 
22  147  1.8  19.1  45.4  5.8 
431  155  6.6  22.3  60.2  6.0 
222  192  3.1  27.6  40.6  7.2 
407  212  15.7  13.8  67.2  11.9 
498  214  6.1  19.4  53.2  12.1 
32  244  12.8  11.6  62.2  17.2 
102  322  12.2  26.9  51.9  8.9 
477  467  33.8  9.0  51.9  30.5 
93  596  55.3  15.7  64.3  13.3 
576  756  48.6  12.7  60.1  19.7 
501  916  102.8  11.5  68.2  15.5 
430  1676  128.2  10.2  50.8  32.1 
103  1938  224.9  27.5  54.5  12.6 
74  9685  941.5  4.5  17.8  71.6 
61  10084  1661.3  7.9  31.5  56.4 
73  13449  1067.7  4.2  15.9  73.4 
14  > 2663  > 7200.0  11.9  73.1  14.5 
151  > 7734  > 7200.0  1.9  22.5  74.8 
410  > 10612  > 7200.0  5.6  40.5  52.6 
560  > 12577  > 7200.0  1.9  29.3  67.5 
730  > 15310  > 7200.0  4.3  36.2  57.5 
From the numbers it is obvious that the inclusion check time takes relatively more time as the number of polyhedra grows.
In contrast, SMT checking for every next point needs relatively less time as the number of polyhedra grows and the time needed for minimization is slowly getting less prominent as well.
It looks like the quadratic time inclusion check is the limiting factor for models with many polyhedra.
Some Profiling of Polyhedral Inclusion Testing
As was described in Sect. 3, the test for inclusion of a newly found polyhedron in the set of already found ones can be sped up by testing containment of the included point which was found by the SMT solver.
Some cursory investigation shows that the number of full checks that still have to be done over the course of the whole computation is about 0.26–0.68 times the number of polyhedra in the solution.
Yet, even though the constant is low, the runtime of the inclusion test is quadratic and it becomes the most dominant part of the computation for large result lists.
5 Conclusion
We presented a novel method to compute tropical equilibria (resp. prevarieties) from an input of sets of polyhedra. We sketched an algorithm for that purpose and discussed several possible improvements. Furthermore, we ran extensive benchmarks with different SMT solvers to compute equilibria of tropicalizations of 46 different BioModels.
The conclusion is that the novel method is working and its computation times compare favorably with a known algorithm using purely polyhedral methods. The runtimes were always smaller for problems that would otherwise take more than 58 seconds to compute, sometimes by a factor of 60 or more. The novel method has also the advantage to be an anytime algorithm, hence it computes more parts of the solution given more time or computation power. This is of importance if computation of the entirety of the solution is infeasible.
The CVC4 SMT solver was overall the fastest solver in this application, yet Yices outperformed CVC4 for models that could be solved in less than 12 seconds.
Preprocessing the input yielded only moderate speedups and only on models that had more than 800 polyhedra as solution. Conversely, preprocessing was always more costly for models with less than 400 polyhedra.
5.1 Future Work
There is obvious potential for a parallel implementation of the procedure. We should expect an almost linear speedup.
The rising percentage of time spent checking for inclusion of already known polyhedra should be addressed. If one could assign onedimensional properties (like dimension) to polyhedra, only parts of the list of already known polyhedra would have to be checked for inclusion.
The computation of nonmaximal polyhedra should be avoided. Hence a good and fast heuristic for choosing the to construct a polyhedron of high dimension would lower the number of rounds needed to compute the solution.
Another avenue is to find better and faster preprocessing to minimize the problem.
5.2 Acknowledgments
This paper is written in grateful memory of the author’s advisor Andreas Weber, who unexpectedly passed away recently.
The author thanks Thomas Sturm, who kickstarted his interest in SMT. Furthermore, he thanks Jörg Zimmermann and Ovidiu Radulescu for discussions and valuable input, and the three anonymous referees for their helpful comments.
This work has been supported by the bilateral project ANR17CE400036 / DFG391322026 SYMBIONT.
Comments
There are no comments yet.