key: cord-0205170-is46i34e
authors: Bessonov, Mariya; Ilmer, Ilia; Konstantinova, Tatiana; Ovchinnikov, Alexey; Pogudin, Gleb; Soto, Pedro
title: Obtaining weights for Gr"obner basis computation in parameter identifiability problems
date: 2022-02-13
journal: nan
DOI: nan
sha: 2f7a8e62c866db9de4f60f4ff1aebae0273668b7
doc_id: 205170
cord_uid: is46i34e

We consider a specific class of polynomial systems that arise in parameter identifiability problems of models of ordinary differential equations (ODE) and discover a method for speeding up the Gr"obner basis computation by using a weighted ordering. Our method explores the structure of the ODE model to generate weight assignments for each variable. We provide empirical results that show improvement across different symbolic computing frameworks

Structural parameter identifiability is a property crucial for designing high-quality mathematical models of real-world phenomena using ordinary differential equations. The question of identifiability arises when one seeks a value for a particular parameter of the model. There can be at most finitely many such values (local structural identifiability) or the value is unique (global structural identifiability). For each case, there exist algorithms that provide solutions with high degree of correctness. We address computational challenges of computing Gröbner bases [4] for polynomial ideals that arise in parameter identifiability. This approach, for instance, lies at the core of global identifiability tool SIAN [19] . The computation itself is well-known to be heavy for some polynomial systems, despite such algorithms as F4 [11] and F5 [12] . We developed an approach that uses additional domain-specific information about the polynomial ideal to significantly speed up the process.

The Gröbner basis (itself and its computation) of a polynomial system can vary based on the order of monomials. The most common and empirically reliable in terms of computing time ordering of monomials is the so-called total-degreereverse-lexicographic order, or tdeg in MAPLE notation. Weighted ordering adds a layer of comparison to monomial orderings where we first compare variables by the weight value multiplied by its degree exponent and then ties are broken according to, for example, tdeg. Properly chosen weights will have tremendous impact on the computation time. Consider the following motivating example showing benefits of weights in general:

Computing the Gröbner basis of this system with tdeg-order of x 1 , x 2 , x 3 , x 4 , x 5 , x 6 The main result of this paper is the rule for assigning weights in parameter identifiability problems based on the data from the input ODE models. That is, given a polynomial system obtained by SIAN from an input ODE model, the rule produces a weight assignment to the variables based on available data from the original ODE that improves the Gröbner basis computation in the vast majority of cases. We provide experimental results showing that improvements in runtime and memory consumption for MAPLE and Magma. We will take advantage of this improvement to accelerate parameter identifiability solving. The implementation of this strategy will be included in an upcoming release of SIAN and SIAN-Julia [28] .

The rest of this paper is organized as follows. In Section 2, we provide an overview of works related to identifiability and Gröbner basis computation. Section 3 contains the weight generation algorithm. In Section 4, we show the experimental results and benchmarks with our new weight assignment approach. Finally, in Section 5, we provide concluding remarks of our work.

Gröbner basis is the key part of SIAN [19] . Buchberger [4] presented the first algorithm to solve the challenge of finding the basis of polynomial systems. This solution can be resource heavy [27] and depends on many factors, such as selection strategy for polynomials [17] . The improvement of the algorithm itself came from the works of Faugère, whose F4 [11] and later F5 [12] are a tremendous success. The key idea of algorithms such as F4 is to take advantage of mathematical tools from the realm of linear algebra. The work [2] addresses the termination and complexity properties of the F5 algorithm. In general, these algorithms are much more efficient than the default Buchberger and each creates a family of solutions to the Gröbner basis problem. For an excellent taxonomy of F5-based approaches, we refer to [9] .

The analysis of connection between weights and homogenization of ideals appeared in [13] and later in more detail in [14] . Homogeneous ideals are an intriguing special case of inputs for a Gröbner basis algorithm. In the mentioned works, weights were used as a homogenization tool. We have observed in the motivating example (1) above that a weighted ordering can break homogenization, offering large benefits. In fact, polynomial systems in SIAN contain non-homogeneous polynomials in most cases by the nature of the input data. For example, an output function (see Definition 1) of the form

where g i is a polynomial, is inevitably transformed into an equation with a free term of degree 0. Therefore, polynomial systems of SIAN always have a non-homogeneous polynomial. By design of the algorithm, the weight of any variable in g i will be 1. For other polynomials that do not have a free term and may be homogeneous, the maximum possible degree in the system will either increase or remain the same. In this sense, we do not necessarily make polynomials "more homogeneous".

Machine learning offers other ways of improving computer algebra algorithms (e.g., an algorithm for Gröbner basis computation). In particular, [30] shows how reinforcement learning improves some aspects of Gröbner basis computation, such as the selection strategy. Machine learning has found applications in other areas of computer algebra, such as choosing order for cylindrical algebraic decomposition [10, 21, 15, 16] or choosing pivots in Gaussian elimination [23] .

The problem of parameter identifiability has been of great interest in the modeling community, and with recent computational advances it gained a lot of attention in the broader scientific computing community. SciML organization recently included an implementation of the global identifiability tool from [8] . The identifiability software SIAN [19] that we consider in this paper has been recently used, for instance, in [6, 32, 33, 31] . Together with the algorithm from [29] , it has been turned into a web-based analyzer package [22] which was applied in, for instance, [25] .

In this section, we present the method for obtaining weighted ordering that improves the runtime of the Gröbner basis computation for polynomial systems formed out of ordinary differential equations. Let us define a model in a state-space form which represents the input to the identifiability software.

Definition 1 (Model in the state-space form). A model in the state-space form is a system

where f = (f 1 , . . . , f n ) and g = (g 1 , . . . , g m ) with f i = f i (x, µ, u), g i = g i (x, µ, u) are rational functions over the complex numbers C. The vector x = (x 1 , . . . , x n ) represents the time-dependent state variables andẋ represents the derivative.

The vectors u = (u 1 , . . . , u s ), y = (y 1 , . . . , y m ), µ = (µ 1 , . . . , µ λ ), and x * = (x * 1 , . . . , x * n ) represent the input variables, output variables, parameters, and initial conditions, respectively. SIAN [19] transforms the input system from the form of Equation (2) to a collection of differential polynomials before assessing global identifiability. To this end, SIAN computes truncated Taylor polynomials of the output functions at time t = 0. For more details about the truncation bound, we refer to [20, Theorems 3.16, 4.12] . The final global identifiability assessment relies on: sampling x * , µ in the polynomial system, calculating y i and its derivatives based on the sample, and using Gröbner basis computation to check uniqueness of x * , µ that provide y i values. This procedure makes SIAN a Monte-Carlo algorithm with user specified probability of correctness. Our goal is to find a weight assignment to each variable of the resulting polynomial system.

3.1.1 The weight assignment rule.

Given a system (2), one can define the Lie derivative L(h) of a function h ∈ C(x, µ, u, u ′ , . . .) with respect to the system by

By applying this formula to each output function y i , we can define, for each state variable or a parameter a ∈ {x, µ} a level as Level(a) := min i ∃ y j ∈ y : a appears in L i (y j ) .

Using that value, we assign weight as follows:

• for a state variable x i ∈ x (and all its derivatives)

Weight(x i ) := Level(x i ) + 1;

• for a parameter µ i ∈ µ:

1, otherwise.

Remark 1 (How we found it?). Empirically, we observe that lower occurrence frequency of a state variable in the right-hand side of the ODE system (excluding output functions) leads to a more successful weight assignment. By success here we mean the assignment's ability to reduce the Gröbner basis algorithm runtime. In our original bruteforce search for weights we found that if a state variable is present in the right-hand side of the output functions y then assigning it a weight w > 1 would decrease performance. Using these observations, we apply a Lie derivative operator L to each output y i . If the differentiation results in no new state variables, we skip the output function y i in the next iteration. We repeat this until all states (and parameters) are found. At the same time, we assign each variable a "level" value corresponding to a derivative order in which it appears. To assign weights, we follow the principle "deeper level means higher weight".

Remark 2. In our search for weights, we aim at minimizing the difference in treatment of state variables and constant parameters of the ODE input. Indeed, any constant parameter can be considered a state with zero derivative. The only exception for constant parameters is that we assign weight to those from the maximal possible level. We only consider locally identifiable constant parameters.

Consider the following ODE system

Differentiating once:ẏ

We see that A, B, x 2 all occur after the first differentiation and hence will have weight of 2. At the same time, state x 1 was already at level 0 and will not be considered further. If we differentiate once more, we geẗ

bringing out C. Differentiating further leads to no new information. The weight assignment then is as follows:

In this section, we present several examples of ODE systems, for which we observe reduction in both the runtime and memory. All simulations were run on a cluster with 64 Intel Xeon CPU with 2.30GHz clock frequency and 755 GB RAM. We ran the computation using MAPLE and Magma computer algebra systems.

Normally, MAPLE's implementation of the F4 algorithm picks multiple different large prime numbers to be the characteristic of the underlying field. In our experiments, we specify a single positive characteristic p = 11863279 for better highlights of any time and memory consumption difference between weighted and non-weighted cases. MAPLE does not directly support the use of weighted orderings with a compiled F4 implementation that is sufficiently fast. To avoid any potential slowdowns, we substitute any variable v in the polynomial system that has weight w greater than 1 with v w . To illustrate this method, if we have a polynomial system

and we wish to use the weight of 2 for variable x, our approach is to compute the basis for a new polynomial system

keeping the variable order as total degree reverse lexicographic. Empirically, there may be a difference observed between computing Gröbner basis with x > y and y > x. In our computations, we order the variables by the degree of the derivative. For example, consider a simple input ODE model SIAN produces the following polynomial system where the double index in x i;j shows that the variable is the j-th derivative of x i in jet-notation.

the order of variables that will provide the best speed without weights is

That is, we order variables from higher to lower derivative degree grouping the same degree together (all order-4 derivatives, then all order-3, etc.)

In what follows, we apply the weights on top of the default variable ordering eq. (8) that has proven itself to be empirically faster. We will consider several ODE models and provide Gröbner basis results over a field of integers with positive prime characteristic p = 11863279. Each example will be summarized by the following metrics in Tables 1 and 2: 1. Number of polynomials in the polynomial system. Once the Gröbner basis computation is finished, the weights are removed by a back substitution to answer the identifiability query. 

We presented an algorithm for generating weighted orderings of variables that we used in Gröbner basis computation for parameter identifiability problems. SIAN [19, 20] generates these polynomial ideals by successive differentiation of output functions. We observed significant improvements for multiple models that vary in complexity, number of polynomials, and number of variables.

Our main idea for weight generation lies in the observation that occurrence of parameters and states in the ODE makes a difference for the effect of a weighted ordering. These empirical observations translated into a sequential Lie differentiation of output functions, mimicking a Breadth-First Search algorithm in graph theory. Effectively, this differentiation produces Taylor coefficients of output functions y in terms of states at a fixed time t = 0 and parameters. We assign weights depending on the depth of these Taylor coefficients, thus, effectively, leveraging the outputs "sensitivity".

In cases where the systems were already relatively quick to return the answer, the weights did not have a negative impact. In fact, in examples where computation slowed down, the memory usage still showed a positive effect, decreasing by around 80%. There has also been a case where the program ran around 44% faster but consumed 30% It is necessary to note that the weight assignment provided here is not unique. In fact, we can even find an alternative assignment given the same weight generation procedure as described earlier. Instead of simply using the rule of higher level Equation (4), we can generate the following assignment:

where M = max Table 3 . Machine learning can offer an automated approach for determining weights. Given a feature set of both the ODE model and the polynomial system, a supervised learning algorithm can output a weight vector or a classification of "Slow" and "Fast" with respect to Gröbner basis. The difficulty of this method lies in determining meaningful features and their representation and, for weights output, knowing the target weights in advance.

Reinforcement learning, in which an algorithm is trained via an interactive reward system, is another possible option. This was explored in [30] for selection strategies in Buchberger's algorithm. For our case, to find best weights, a reinforcement learning agent would sample weight values, compute the basis, aggregate reward if the computation is improved or lose the reward otherwise. So far in our experiments, this approach has not yielded an improving weight combination that would outperform default order. Another concern is poor generalization properties of this method where one would have to retrain an agent for each new input system. This could be addressed by expanding the state space to include the input ODE systems, so that the agent could make a step by not only choosing weights but also adjusting inputs.

Finally, we discuss the possible effect weights have on the F4 algorithm. The F4 algorithm operates by creating a matrix with rows representing monomials up to a certain degree [11] . From the log outputs of MAPLE's Gröbner basis function, we observed that without weights, the algorithm spends a considerable amount of time on a degree value while with weights there is either no such delay or it is postponed to a much later stage (higher degree). Assume that computation slows down at degree n. After we introduce a weight of w, whatever monomial caused the delay will no longer be present, since degrees have changed (recall that we apply weight by raising variables to the power w). The problematic computation then moves to higher degrees.

In this section, we present details about models considered. Specifically, we will describe the differential equations and the resulting weights of the models used in the analysis of this paper.

The first model we consider is provided in Equation (A.1). The system comes from [18] and describes time periodicity in cell behavior. This example has 4 state variables x 1,2,3,4 and 6 parameters.

The differentiation algorithm applied to this model returns the following weight distribution x 1 ⇒ 1, x 4 ⇒ 2, x 2 , x 3 ⇒ 3, β ⇒ 4. This indicates that, e.g., x 2 and all of its derivatives receive a weight of 3. On the other hand, only constant β receives weight of 4.

SIAN uses an auxiliary variable z aux to account for the presence of denominators in the right-hand side of the original input ODE system. We observe that giving a weight of at most 3 to this variable does not decrease performance.

In this next example, we considered a SEIR-model of epidemics from [26, table 2, ID=14]. The example originally had 3 output functions. We reduced it to 1 to create more of a computational challenge for our program. We also use the term µ i s instead of µ iµ s in the third equation. The state-space form of the model is presented in Equation (A.2).

As a result of our weight ordering pick, we assign these weights:

C ⇒ 1, γ, δ, E ⇒ 4, I ⇒ 3, S ⇒ 2

The following model also comes from [26] 

In this model, the computation without weights has not finished in reasonable time, consuming all available memory. Our algorithm returns the following weight assignment:

A d , A n , S n , β a,i , h 1,2 , γ ai , f, ǫ a,s ⇒ 2, I n , S d ⇒ 1

This is a biomedical model applied to COVID-19 present in [26] . The outputs were changed to make the system more of a computational challenge to SIAN.

The weights for this system were assigned as follows:

ρ ⇒ 3, E ⇒ 2, I, S ⇒ 1

The next two SEIR models were presented in [26] in examples 34 and 16. Example 34 is presented in Equation (A.5) E ⇒ 2, γ ⇒ 4, S, ρ, β ⇒ 3. I, R ⇒ 1.

The following model was presented in [5] . This is a SIR-model with an oscillating forcing term given by equations for

The weight assignment here is as follows:

We considered two HPV models studied in [3] . The model itself is given by eq. Output set 1 weights:

Output set 2 weights:

This model comes from [24] and was used for identifiability analysis in [1] . The ODE system consists of 15 equations, eq. (A.13), and the outputs, eq. (A.14).

We specify values of the following parameters based on the known information from [24] to reduce the number of target identifiability candidates a 1 , a 2 , a 3 , c 1a , c 5a , c 1c , c 3c , c 2c , c 1 , c 2 , c 3 , c 4 , e 1a , k v .

The output functions of eq. (A.13) yields the weights as below (states not listed below get weight of 1)

This model comes from [7] describing pharmacokinetics of glucose-oxidase. We make one modification setting a 1 = a 2 . The model is small but presents a significant computational challenge for global identifiability, that is, it is very difficult to compute Gröbner basis of this model's polynomial system in SIAN.

The weight assignment for this system is

A.10 Example with slowdown: a SIR-model

In eq. (A.16), we present an example in which the weight assignment generated by our algorithm increases the running time of F4. This example comes from [26, Table 1 , ID 26] . While running it in MAPLE, we observed an increase in CPU time from around 12 to 50 minutes. The memory usage slightly decreases from 11.5 to 10.8 GB. In Magma, this system shows a larger increase in memory from 5.6 to 18.8 GB with an increase in CPU time from around 7 to 32 minutes.

= b N − S (I λ + λ Q ǫ a ǫ q + λ ǫ a A + λ ǫ j J + d + 1), I = k 1 A − (g 1 + µ 2 + d 2 ) I, R = g 1 In + g 2 J − d 3 R, A = S (I λ + λ Q ǫ a ǫ q + λ ǫ a A + λ ǫ j J) − (k 1 + µ 1 + d 4 ) A, Q = µ 1 A − (k 2 + d 5 ) Q, J = k 2 Q + µ 2 I − (g 2 + d 6 ) J, y 1 = Q, y 2 = J (A. 16) B Additional tables B.1 Inverted weights

In Table 3 , we provide MAPLE benchmarks for the weights given by Equation (9) . We see that for some systems, such as eq. (A.1), there is a considerable speedup, even compared to the main weight assignment of this paper. The eqs. (A.8) and (A.10) model, however, does not finish. Overall, this illustrates that while this may be inferior to the weight assignment by Equation (4), there are multiple ways to create an efficient weight distribution. The difference in improvement factors is presented in Figures 5 Figure 6 : Memory improvement difference between inverted (eq. (9)) and regular weight assignments when running MAPLE's Gröbner basis.

An iterative identification procedure for dynamic modeling of biochemical networks

On the complexity of the F5 Gröbner basis algorithm

Transmission heterogeneity and autoinoculation in a multisite infection model of HPV

A theoretical basis for the reduction of polynomials to canonical forms

Parameter estimation of some epidemic models. The case of recurrent epidemics caused by respiratory syncytial virus

Estimating vaccination threshold and impact in the 2017-2019 hepatitis A virus outbreak among persons experiencing homelessness or who use drugs in

Effect of prosthetic sugar groups on the pharmacokinetics of glucose-oxidase

Differential elimination for dynamical models via projections with applications to structural identifiability

A survey on signature-based algorithms for computing Gröbner bases

Comparing machine learning models to choose the variable ordering for cylindrical algebraic decomposition

A new efficient algorithm for computing Gröbner bases (F4)

A new efficient algorithm for computing Gröbner bases without reduction to zero (F5)

On the complexity of computing Gröbner bases for quasihomogeneous systems

On the complexity of computing Gröbner bases for weighted homogeneous systems

Algorithmically generating new algebraic features of polynomial systems for machine learning

Improved cross-validation for classifiers that make algorithmic choices to minimise runtime without compromising output correctness

One sugar cube, please" or selection strategies in the Buchberger algorithm

Oscillatory behavior in enzymatic control processes

SIAN: software for structural identifiability analysis of ODE models

Global identifiability of differential models

Applying machine learning to the problem of choosing a heuristic to select the variable ordering for cylindrical algebraic decomposition

Web-based Structural Identifiability Analyzer

Good pivots for small sparse matrices

Mathematical model of NF-kB regulatory module

Quantification of Type I Interferon Inhibition by Viral Proteins: Ebola Virus as a Case Study

Structural Identifiability and Observability of Compartmental Models of the COVID-19 Pandemic

The complexity of the word problems for commutative semigroups and polynomial ideals

SIAN: Structural Identifiability Analyzer

Computing all identifiable functions of parameters for ODE models

Learning selection strategies in Buchberger's algorithm

Long-term regulation of prolonged epidemic outbreaks in large populations via adaptive control: a singular perturbation approach

Delicate Balances in Cancer Chemotherapy: Modeling Immune Recruitment and Emergence of Systemic Drug Resistance

An integrated framework for building trustworthy data-driven epidemiological models: Application to the COVID-19 outbreak in New York City

The authors are grateful to CCiS at CUNY Queens College for the computational resources and to Andrew Brouwer for pointing out the HPV example. This work was partially supported by the National Science Foundation under grants CCF-1563942, CCF-1564132, DMS-1760448, DMS-1853650, and DMS-1853482.