Model selection for dynamical systems via sparse regression and information criteria

TL;DR

This paper introduces a novel model selection algorithm for dynamical systems that efficiently evaluates a large number of candidate models using sparse regression and information criteria, enabling automatic identification of the most supported model.

Contribution

It presents a new approach combining sparse identification with Pareto frontier sub-selection to compute information criteria for many models efficiently.

Findings

Successfully identified canonical dynamical systems like Lorenz and SEIR models.

Accurately ranked models using AIC scores into support categories.

Demonstrated the method's effectiveness on complex dynamical systems.

Abstract

We develop an algorithm for model selection which allows for the consideration of a combinatorially large number of candidate models governing a dynamical system. The innovation circumvents a disadvantage of standard model selection which typically limits the number candidate models considered due to the intractability of computing information criteria. Using a recently developed sparse identification of nonlinear dynamics algorithm, the sub-selection of candidate models near the Pareto frontier allows for a tractable computation of AIC (Akaike information criteria) or BIC (Bayes information criteria) scores for the remaining candidate models. The information criteria hierarchically ranks the most informative models, enabling the automatic and principled selection of the model with the strongest support in relation to the time series data. Specifically, we show that AIC scores place each candidate model in the {\em strong support}, {\em weak support} or {\em no support} category. The method correctly identifies several canonical dynamical systems, including an SEIR (susceptible-exposed-infectious-recovered) disease model and the Lorenz equations, giving the correct dynamical system as the only candidate model with strong support.

Equations

$$AI{C}_j=2k-2\ln (L(\mathbf{x},\hat{\mu })),$$

$$AI{C}_c=AIC+2(k+1)(k+2)/(m-k-2),$$

$$\frac{d}{dt}\mathbf{x}(t)=\mathbf{f}(\mathbf{x}(t)),$$

$$\frac{d}{dt}\mathbf{X}=\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }.$$

$$\dot{x}=x-0.2x^3-0.1x^4.$$

Full text

1 Introduction

Nonlinear dynamical systems theory has provided a fundamental characterization and understanding of phenomenon across the physical, engineering and biological sciences. Traditionally, simplified models are posited by domain experts, and simulations and analysis are used to explore the underlying dynamical behavior which may include chaotic dynamics (e.g. Lorenz equations), nonlinear oscillations (e.g. van der Pol, Duffing), and/or bifurcations. The emergence of data-driven modeling methods provides an alternative framework for the discovery and/or inference of governing nonlinear dynamical equations. From this perspective, governing models are posited from time-series measurement data alone. The recent sparse identification of nonlinear dynamics (SINDy) method [1] uses sparse regression and a Pareto analysis to correctly discover parsimonious governing equations from a combinatorially large set of potential dynamical models. This methodology can be generalized to spatio-temporal systems [2, 3] and dynamical systems characterized with rational function nonlinearities which often occur in biological networks [4]. Although previously suggested [4], no explicit connection between the sparse selection process and information theoretic criteria has been established. Information criteria is the standard statistical method established for the model selection process. In this manuscript, we demonstrate that the Akaike information criteria (AIC) and/or Bayesian information criteria (BIC) can be connected with the SINDy architecture to hierarchically rank models on the Pareto front for automatic selection of the most informative model. As outlined in Fig. 1, the AIC/BIC metrics can be used to correctly infer dynamical systems for a given time-series data set from a combinatorially large set of models. To our knowledge, this is the first explicit demonstration of how information theory can be exploited for the identification of dynamical systems.

Model selection is a well established statistical framework for selecting a model from a set of candidate models given time series data, with information theory providing a rigorous criteria for such a selection process. As early as the 1950s, a measure of information loss between empirically collected data and model generated data was proposed to be computed using the Kullback-Leibler (KL) divergence [5, 6]. Akaike built upon this notion to establish a relative estimate of information loss across models that balances model complexity, and goodness-of-fit [7, 8]. This allowed for a principled model selection criteria through the Akaike information criteria (AIC). The AIC was later modified by G. Schwarz to define the more commonly used Bayes information criteria (BIC) [9]. Both AIC and BIC compute the maximum log likelihood of the model and impose a penalty: AIC adds the number of free parameters $k$ of the posited model, while BIC adds half of $k$ multiplied by the log of the number of data points $m$. Much of the popularity of BIC stems from the fact that it can be rigorously proved to be a consistent metric [9]. Thus if a number of models $q$ are proposed, with one of them being the true model, then as $m\rightarrow\infty$ the true model is selected as the correct model with probability approaching unity. Regardless of the selection criterion, AIC or BIC, they both provide a relative estimate of information loss across a selection of $m$ models, balancing model complexity and goodness-of-fit [10].

Although successful and statistically rigorous, model selection in its standard implementation is typically performed on $q$ predetermined candidate models, where $q$ is often ten or less[10, 11, 12, 13, 14, 15]. For modern applications to dynamical systems where rich, high-fidelity time series data can be acquired, the restriction on the number of models limits the potential impact of AIC/BIC metrics for discovering the correct nonlinear dynamics. Instead, it is desired to consider a combinatorially large set of potential dynamical models as candidates, thus enforcing that $q\gg 1$. This is computationally intractable with standard model selection, as each of the models from the combinatorially large set would have to be simulated and then evaluated for a BIC/AIC score.

As an alternative, sparse regression techniques, embodied by the lasso (least absolute shrinkage and selection operator) method of Tibshirani [16], have enabled variable selection algorithms capable of optimally choosing among a combinatorially large set of potential predictors. Specifically, a lasso regression analysis, or one of its many generalizations and variants, performs both variable selection and regularization in order to enhance the prediction accuracy and interpretability of the statistical model it produces. Such a mathematical tool provides a critically enabling framework for model selection, in particular for identifying dynamical systems.

In this work we demonstrate a new mathematical framework that leverages information criteria for model selection with sparse regression for evaluating a combinatorially large set of candidate models. Specifically, we circumvent a direct computation of information criteria for the combinatorially large set of models by first sub-selecting candidate functional forms from which are most consistent with the time series data. Thus we integrate two maturing fields of statistical analysis: (i) sparse regression for nonlinear systems identification via SINDy and (ii) model selection via information criteria. Our algorithm is demonstrated to produce a robust procedure for discovering parsimonious, nonlinear dynamical systems from time series measurement data alone. We demonstrate the methodology on a number of important examples, including the SEIR (susceptible-exposed-infectious-recovered) disease model and the Lorenz equations, and demonstrate its efficacy as a function of noise, length of time series and other key regression factors. Our sparse selection of dynamical models from information theory criteria ranks the candidate models and further shows that the correct model is strongly supported by the AIC/BIC metrics. Ultimately, the method provides a cross-validated and ranked set of candidate nonlinear dynamical models for a given time-series of measurement data, thus enabling data-driven discovery of the underlying governing equations.

2 Background

2.1 Model selection via information criteria

The process of model selection fundamentally enables the connection of observations or data to a mathematical model. Further, a well-selected model, which describes a governing law or physical principle underlying the system process, can be utilized for prediction outside of the sampled data and parameter configuration [10]. The substantial challenge facing the selection process is discovering the best predictive model from a combinatorially large space of available models. To emphasize the enormity of this task, consider the number of possible polynomial models up to degree $4$ with $5$ state variables. Approximately $10^{38}$ models would need to be constructed, fit to the data, and compared according to a goodness-of-fit metric [4]. Thus, model selection quickly becomes computationally intractable for a modest number of variables and polynomial degree.

Typically, a sub-selection of models occurs based on prior scientific knowledge of the process to produce a subset, $\mathcal{O}(10)$, of heuristically defined candidate models [10, 11, 12, 13, 14, 15]. Recent research has focused on automatically expanding the number of candidate models [17, 18, 19]. Once a subset of models is chosen, the model selection procedure balances the goodness-of-fit with model complexity, i.e. the number of free parameters. A wide-variety of rigorous statistical metrics have been developed to balance model parsimony and predictive power including popular methods such as the Akaike information criterion (AIC) [7, 8], Bayesian information criterion (BIC) [9], cross-validation (CV) [20], deviance information criterion (DIC) [21], and minimum description length (MDL) [22]. Methods such as AIC explicitly balance parsimony and relative information loss across models, penalizing the number of parameters in the model to avoid overfitting.

In this manuscript, we utilize the ubiquitous and well-known AIC as the statistical criterion for comparing candidate models. The AIC value for each candidate model $j$ is:

[TABLE]

where $L(\mathbf{x},{\mu})=P(\mathbf{x}|\mu)$ is the likelihood function (conditional probability) of the observations $\mathbf{x}$ given the parameters $\mu$ of a candidate model, $k$ is the number of free parameters to be estimated, and $\hat{\mu}$ is the best-fit parameter values for the data [7, 8]. In practice, the AIC requires a correction for finite sample sizes given by

[TABLE]

where $m$ is the number of observations. A common likelihood function uses the residual sum of squares (RSS), given by $RSS=\sum_{i=1}^{m}(y_{i}-g(x_{i};\mu))^{2}$, where $y_{i}$ are the observed outcomes, $x_{i}$ are the observed independent variables, and $g$ is the candidate model. The RSS is a well-known objective function for least squares fitting. In this case AIC can be expressed as $AIC=m\ln(RSS/m)+2k$  [10]. Note that (1) penalizes, by increasing the AIC score, the models that have a large number of free parameters and which are unable to capture the characteristics of the observed data.

2.2 SINDy and sparse model selection

Identifying dynamical systems models from data is increasingly possible with access to high-fidelity data from simulations and experiments. With traditional methods, only a small handful of model structures may be posited and fit to data via regression. Indeed, simultaneous identification of both the structure and parameters of a model generally requires an intractable search through combinatorially many candidate models. Genetic programming has been recently used to determine the structure and parameters of dynamical systems [23, 17, 24] and control laws [25], enabling the efficient search of complex function spaces. Sparsity-promoting techniques have also been employed to simultaneously identify the structure and parameters of a dynamical system model. Compressed sensing was first used to determine the active terms in the dynamics [26], although it does not work well with overdetermined systems that arise when measurements are abundant. In contrast, the sparse identification of nonlinear dynamics (SINDy) algorithm [1] uses sparsity-promoting regression, such as the lasso [16] or sequential thresholded least-squares algorithm [1], to identify nonlinear dynamical systems from data in overdetermined situations.

Here we review the SINDy architecture for identifying nonlinear dynamics from data. The general observation underlying SINDy is that most dynamical systems of a state $\mathbf{x}\in\mathbb{R}^{n}$,

[TABLE]

have only a few active terms in the dynamics, making them sparse in a suitable function space. To identify the structure and parameters of the model, a set of candidate symbolic functions are first concatenated into a library $\bm{\Theta}(\mathbf{x})=\begin{bmatrix}\theta_{1}(\mathbf{x})&\cdots&\theta_{p}(\mathbf{x})\end{bmatrix}$. With time-series data ${\bf X}\in\mathbb{R}^{m\times n}$, where each row is a measurement of the state $\mathbf{x}^{T}(t_{k})$ in time, it is possible to evaluate the candidate function library $\bm{\Theta}(\mathbf{X})\in\mathbb{R}^{m\times p}$ at the $m$ time points. Finally, with derivative data $\dot{\bf{X}}\in\mathbb{R}^{m\times n}$, either measured or obtained by numeric differentiation, it is possible to pose a regression problem:

[TABLE]

The few active terms in the dynamics, given by the nonzero entries in the columns of $\bm{\Xi}$, may be identified using sparse regression. In particular, the sparsest matrix of coefficients $\bm{\Xi}$ is determined that also provides a good model fit, so that $\|\dot{\bf{X}}-\bm{\Theta}({\bf X})\bm{\Xi}\|_{2}$ is small. Sparse regression has the added benefit of avoiding overfitting, promoting stability and robustness to noise. Since the original SINDy method, there have been numerous innovations and extensions to handle rational function nonlinearities [4], partial differential equations [2], highly corrupted data [27], and to build Galerkin regression models in fluids [28]. The method is also connected to the dynamic mode decomposition (DMD) [29] if only linear functions are used in $\bm{\Theta}$.

In most sparse regression algorithms, there is a parameter that determines how aggressively sparsity is promoted. The successful identification of the model in Eq. (3) hinges on finding a suitable value of this sparsity-promoting parameter. Generally, a the parameter value is swept through, and a Pareto front is used to select the most parsimonious model. However, the Pareto frontier may not have a sharp elbow or may instead have a cluster of models near the elbow, thus compromising the automatic nature of the model selection using SINDy alone.

3 Methods

Our algorithm integrates sparse regression for nonlinear system identification with model selection via information criteria. This approach enables the automatic identification of a single best-fit model from a combinatorially large model space. In the first step of the algorithm, the SINDy method provides an initial sub-selection of models from a combinatorially large number of candidates. The sub-selection of candidate models near the Pareto frontier is critically enabling as it is computationally intractable to simulate and compare against the time series data. Importantly, the sub-selection can take the number of candidate modes from $10^{9}$ (for our 2-D cubic example) to a manageable $10^{2}$. This then allows for a tractable computation of AIC or BIC scores for the remaining candidate modes. The information criteria hierarchically ranks the most informative models, enabling the automatic selection of the model with the strongest support. This is in contrast to a standard Pareto front analysis which looks for a parsimonious model at the elbow of the error versus complexity curve. However, the Pareto frontier may not have a sharp elbow or may instead have a cluster of models near the elbow. Algorithm 1, using AIC as our information criteria, is executed for model selection. Figure 1 illustrates the algorithm in practice.

When evaluating dynamical systems models, there is some ambiguity about what constitutes an ”observation.” We take time series data for a given set of initial conditions to be an observation, rather than taking each measurement at each time point. To obtain a representative error for the $i$th time-series observation, we calculate the average absolute error over the entire time series: $E_{\mbox{avg}}=\sum_{\tau}|y_{i}^{\tau}-g(x_{i}^{\tau};\mu)|$. We then substitute this representative error in for $(y_{i}-g(x_{i};\mu)$, and take the sum of the squares so that $AIC=m\>\ln((\sum_{i=1}^{m}E_{avg}(x_{i},y_{i}))/m)+2k$. We use this value for AIC in (2).

The candidate models with the lowest scores are ranked as the most likely. To be more precise, the AIC scores for each candidate model can have a wide range of values which require a rescaling by the minimum AIC value ($AIC_{min}$) [10, 30]. The rescaled AIC values $\Delta_{j}=AIC_{i}-AIC_{min}$ can be directly interpreted as a strength-of-evidence comparison across models. Models with $\Delta_{j}\leq 2$ have strong support, $4\leq\Delta_{j}\leq 7$ have weak support, and $\Delta_{j}\geq 10$ have no support [30]. These rankings allow for a principled procedure for retaining or rejecting models within the candidate pool of models. For time series data with enough data samples, large enough signal to noise and/or sufficiently large set of candidate models, only one model is typically strongly supported by AIC metrics. In the examples used to demonstrate the method, the only strongly supported candidate model is, in fact, the correct model.

4 Results: Model selection

4.1 1-D polynomial model

We demonstrate the relationship between SINDy model selection and AIC ranking procedure, as illustrated in Fig. 1 d, with a simple single state variable model with three polynomial terms,

[TABLE]

We compute three time-series for this model, as shown in Fig. 2 a.

Assuming a feature library with up to 5th order polynomials, or $d=6$, there are $\sum_{i=1}^{6}\binom{6}{i}=63$ possible models. A representation of the combinatorial set of models is shown in Fig. 2 b, and we perform least-squares fitting for the coefficients of each model. All models are cross-validated using 100 initial conditions, and relative AICc scores are calculated. This combinatorially complete set of models (black dots) populate the AICc Pareto front in Fig. 2 c. The relative AICc metric successfully characterizes the true, 3 term model as the only model with ”strong support.” All other models within the full combinatorial space fall within the ”no support” range.

SINDy enables us to sub-select a set of models from the feature library and, for this simple system, we can compare against the combinatorial set. The models selected by SINDy are highlighted in blue in Fig. 2 c and d. SINDy finds the correct 3 term model, as well as models with 1, 2, 4, 5, and 6 terms each. These models and coefficients exactly match a subset of those found by combinatorial least-squares, which is expected given that this implementation of SINDy uses least-squares to fit the coefficients. Notably, the 1 and 2-term models selected by SINDy are the most meaningful reductions of the true underlying system (with smaller coefficients set to zero), rather than the 1 and 2 term models with the lowest error in the full feature space.

4.2 2-D cubic model

For larger systems, enumerating all possible models represented in a given feature library would be computationally infeasible, but by using SINDy a sub-selection of the most relevant models are generated and selected. For example, consider the relatively simple example of a 2-state variable system ($n=2$) with a 5th order polynomial library ($d=6$). In this case there are $N_{m}={n+d\choose d}=28$, possible monomials, and $N_{p}=\sum_{i=1}^{N_{m}}{N_{m}\choose i}=268,435,455$ potential models. Figure 3 demonstrates the results of performing SINDy and AIC evaluation on a cubic model.

With only a single time-series for each state variable (x and y) as input (Fig. 3 a), the SINDy selected models with varying number of terms (blue circles) include the correct model (magenta circle). Using 100 randomly selected initial conditions for cross-validation, relative AICc to ranks the correct model as strongly supported, and all other models as having no support.

4.3 3-D Disease transmission model

Next we apply our method to the susceptible-exposed-infectious-recovered (SEIR) disease transmission model. Models of this type are often used to determine disease transmission rates, detect outbreaks and develop intervention strategies. However, generating the appropriate model for interactions between different populations is currently done heuristically and then evaluated using information criteria [15]. Using SINDy with AIC evaluation would provide data-driven model generation and selection from a wider library of possible interaction terms. As a first step, we apply SINDy with AIC to a discrete, deterministic SEIR model as shown in Fig. 4.

We input a single time series representing an outbreak for the $S$, $E$, and $I$ state variables (n =3), and provide a library of polynomial terms up to 2nd order (d= 3). For this example, the total number of models represented in the library is $N_{p}=\sum_{i=1}^{N_{m}}{N_{m}\choose i}=1023$ with $N_{m}={n+d\choose d}=10$ possible monomials. A complication of the SEIR system is that $R$ is a redundant state variable; $S$, $E$, and $I$ have no dependence on $R$, and $R$ depends only on a term already represented in the $I_{k+1}$ equation $R_{k+1}=R_{k}+0.04I_{k}$. A result of this redundancy, SINDy cannot find the correct equations with $R$ included in the library. Without $R$, SINDy selects a set of models from 1023 possible in the library (blue circles), and with 100 cross-validation measurements the relative-AICc evaluation ranks only the correct 6-term model as having any support (magenta circle) in Fig. 4 b.

4.4 Lorenz model

As a final example, we demonstrate SINDy with AIC on the chaotic Lorenz system [31]. Using a library of polynomials up to 2nd order (d=3) for the 3-state variable system (n = 3), there are once again $N_{p}=1023$ models represented in the function library. Providing one time-series for each state variable, as shown in Fig. 5 a, SINDy recovers a subset of these models (circles in Fig. 5 b). Using 100 cross-validation measurements on each model, the relative-AICc criteria ranks the correct 7-term model as having strong support (magenta circle). Unlike in previous examples, one other model is ranked as having ”weak support” (blue circle). This model has an additional small constant term in the equation for $x$: $\dot{x}=8.5\times 10^{-6}+10(y-x)$.

A possible explanation for the level of support for this model is the chaotic nature of the Lorenz system. Even when the recovered model is correct, small variations in the recovered coefficients ($\approx 10^{-6}$ for this case), will cause the calculated time-series for the recovered model to diverge from the ”true” model after some length of time ($>1$ unit time for these parameters). For the example in Fig. 5 b, the cross-validation uses time series of length $t=5$ (arb. units). In a true model-selection situation, we would not know this characteristic length scale ahead of time, and a sensitivity analysis would need to be preformed. We discuss this and other challenges to practical implementation in the next section.

5 Practical implementation: noise and number of measurements

SINDy with AIC ranking can successfully select the correct model for a variety of known systems, given low enough measurement noise and a large amount of data for cross-validation. Under practical conditions, the signal-to-noise may be lower than desired and the amount of data available for cross-validation may be restricted. In Fig. 6 we show the effects of increasing noise and number of cross-validation experiments on the selection of the correct model for the Lorenz system. Reading from the top of Fig. 6 a, for low noise, $\epsilon=0.1$, only the correct model (magenta) falls within the supported (strong or weak) range. Increasing the noise to $\epsilon=0.2$ and $0.5$ cause other models to descend into the weak support regime and eventually into the strong support range. Around this level of noise, the relative AICc scores for the incorrect models are very sensitive to the random additive noise. Repeating the computation for different instances of randomly generated measurement noise causes the position of these models to fluctuate (data not shown), although the true model maintains the lowest score (over 10 instances). This suggests that data sub-sampling could be used to test for models with noise-fitted terms.

Above a certain level of noise ($\epsilon=1$) the method is unable to robustly select the correct model, and for even higher levels of noise ($\epsilon=5$) a larger number of incorrect models appear to have support. For the particular case of the Lorenz system used here, SINDy is unable to select the correct model, and therefore the true model is not evaluated by AICc at these noise levels. As we are using a relative AICc score, the score of the lowest model will always be zero, even when the actual error between that model and the data is very high. These examples highlight the importance of examining the error between the model-generated time-series against and the data in addition to the relative AICc metric.

The number of time-series used for cross-validation also has an impact on the relative AICc score between the true model and other models given strong support as shown in Fig. 6 b. For $\epsilon=0.5$ and only 2 time-series used for cross-validation and AICc calculation, the correct model is selected by SINDy, but a simpler (and incorrect) 1-term model is assigned a lower score. Given only 2 more time-series for cross-validation, ($N_{cross}=4$), the correct model has the lowest score. Increasing the number of cross-validation measurements used for evaluation to $N_{cross}=60,80,$ and $100$ increases the relative score of the incorrect 8-term model to the true 7-term model.

The success of the method also relies on the duration and sampling of the time series, especially in the case of chaotic systems like the Lorenz model. Here the cross-validation experiments run from $t=0$ to $t_{end}=1$ (arb. units), as opposed to those in Fig. 5 which ran to $t_{end}=5$. Even with higher noise ($\epsilon=0.1$ compared to $\epsilon=0.001$), only the true model is supported. Again, sensitivity to sub-sampling of data can help differentiate between noise-fitting and mechanistically essential terms.

6 Discussion and conclusions

The integration of mathematical techniques advocated here, (i) sparse regression for nonlinear systems identification via SINDy and (ii) model selection via information criteria, provides a new paradigm for model selection of dynamical systems. Algorithmically, the critical methods combine as follows. In the first step of the algorithm, the SINDy method, which is based upon sparse regression, provides an initial sub-selection of models from a combinatorially large number of candidate models. The selection of candidate models is critically enabling as it reduces the number of potential models to a manageable number, which can each be evaluated through simulation and comparison to the time series data. Indeed, the remaining candidate models, which are now on the order of ten models, are each evaluated using information criteria such as AIC or BIC.

The candidate models with the lowest scores are ranked as the most likely. Specifically, in what follows, we work with the AIC score and show that these scores place each candidate model in the strong support, weak support or no support category. For time series data with enough data samples, large enough signal to noise and/or a sufficiently large set of candidate models, only one model is typically strongly supported by AIC metrics. In the examples used to demonstrate the method, the only strongly supported candidate model is, in fact, the correct model.

The method presented provides an important contribution to standard model selection as well as to the SINDy paradigm. In particular, each of these methods has a significant shortcoming. In model selection, the shortcoming is centered around the inability of the standard AIC/BIC criteria to assess a combinatorially large set of candidate models. For SINDy, the sparse selection process for identifying the underlying dynamical systems lacks a principled method for selecting the correct dynamical model. The algorithm here circumvents both of these shortcomings. Specifically, the sparse regression of SINDy allows for the consideration of a combinatorially large number of candidate models. The sub-selected set of models can then each be evaluated using information criteria to select the correct dynamical system. The connection between information criteria and automatic model selection can also be integrated with genetic algorithms for selecting the structure and parameters of dynamical systems [23, 17, 24, 25]. The process can be automated for data-driven discovery of physical principles and laws of motion, which is now often referred to as the 4th paradigm of science [32].

Acknowledgements

JNK acknowledges support from the Air Force Office of Scientific Research (FA9550-15-1-0385). SLB and JNK acknowledge support from the Defense Advanced Research Projects Agency (DARPA contract HR0011-16-C-0016). JLP and NMM would like to thank Bill and Melinda Gates for their active support of the Institute for Disease Modeling and their sponsorship through the Global Good Fund.