How Entropic Regression Beats the Outliers Problem in Nonlinear System Identification

TL;DR

Entropic Regression is a robust nonlinear system identification method that leverages information theory to effectively handle noise and outliers, outperforming existing techniques especially in high-dimensional and complex systems.

Contribution

The paper introduces Entropic Regression, a novel information-theoretic approach for nonlinear system identification that addresses limitations of current methods regarding outliers, parameter sensitivity, and high-dimensional systems.

Findings

Outperforms state-of-the-art methods in robustness and accuracy.

Effectively handles noise and outliers in system data.

Successfully applied to chaotic systems like Lorenz and Kuramoto-Sivashinsky.

Abstract

In this work, we developed a nonlinear System Identification (SID) method that we called Entropic Regression. Our method adopts an information-theoretic measure for the data-driven discovery of the underlying dynamics. Our method shows robustness toward noise and outliers and it outperforms many of the current state-of-the-art methods. Moreover, the method of Entropic Regression overcomes many of the major limitations of the current methods such as sloppy parameters, diverse scale, and SID in high dimensional systems such as complex networks. The use of information-theoretic measures in entropic regression poses unique advantages, due to the Asymptotic Equipartition Property (AEP) of probability distributions, that outliers and other low-occurrence events are conveniently and intrinsically de-emphasized as not-typical, by definition. We provide a numerical comparison with the current state-of-the-art methods in sparse regression, and we apply the methods to different chaotic systems such as the Lorenz System, the Kuramoto-Sivashinsky equations, and the Double Well Potential.

Equations

$$\dot{\mathbf{z}}=\mathbf{F}(\mathbf{z}),$$

$${\dot{z}}_i=F_i(\mathbf{z})=\sum _{k=0}^{\infty }{a}_{ik}{\varphi }_k(\mathbf{z}),$$

$$\mathbf{\varphi }=[{\varphi }_0(\mathbf{z}),{\varphi }_1(\mathbf{z}),{\varphi }_2(\mathbf{z}),\dots ]=[1,z_1,z_2,\dots ,z_N,z_1{z}_2,z_1{z}_3,\dots ,z_{N-1}{z}_N,\dots ]$$

$${\hat{F}}_i(\hat{\mathbf{z}}(t))=\sum _{k=0}^K{a}_{ik}{\varphi }_k(\hat{\mathbf{z}}(t))+{\xi }_i(t),(t=t_0,\dots ,t_{\ell };i=1,\dots ,N).$$

$$\displaystyle\left(\begin{array}[]{ccc}|&|&|\\ \dot{z}_{1}(t_{i})&\dots&\dot{z}_{N}(t_{i})\\ |&|&|\end{array}\right)$$

$${\mathbf{f}}^{(i)}=\Phi {\mathbf{a}}^{(i)}+{\mathbf{\xi }}^{(i)},$$

$$\mathbf{f}=\Phi \mathbf{a}+\mathbf{\xi },$$

$$\{\begin{array}{l}{k}_1={\mathrm{argmax}}_{k}I(\Phi R(\Phi ,\mathbf{f},\{k\});\mathbf{f}),\\ \hat{\mathbf{a}}=R(\Phi ,\mathbf{f},k_1),\\ {\mathbf{z}}_1=\Phi R(\Phi ,\mathbf{f},k_1).\end{array}$$

$$\{\begin{array}{l}{k}_i={\mathrm{argmax}}_{k\notin \{k_1,\dots ,k_{i-1}\}}I(\Phi R(\Phi ,\mathbf{f},\{k\});\mathbf{f}|{\mathbf{z}}_{i-1}),\\ \hat{\mathbf{a}}=R(\Phi ,\mathbf{f},\{k_1,\dots ,k_i\}),\\ {\mathbf{z}}_i=\Phi R(\Phi ,\mathbf{f},\{k_1,\dots ,k_i\})\end{array}$$

$$\{\begin{array}{l}\hat{\mathbf{a}}=R(\Phi ,\mathbf{f},\{k_1,\dots ,k_i\}/\{k_i\}),\\ {\bar{\mathbf{z}}}_j=\Phi R(\Phi ,\mathbf{f},\{k_1,\dots ,k_i\}/\{k_i\}),\end{array}$$

$${\hat{z}}_k(t_i)=z_k(t_i)+{\eta }_{ki},$$

$$\{\begin{array}{l}{\dot{z}}_1=F_1(\mathbf{z})=\sigma (z_2-z_1),\\ {\dot{z}}_2=F_2(\mathbf{z})=z_1(\rho -z_3)-z_2,\\ {\dot{z}}_3=F_3(\mathbf{z})=z_1{z}_2-\beta z_3,\end{array}$$

$$u_t=-\nu u_{xxxx}-u_{xx}+2u{u}_x,(t,x)\in [0,\infty )\times (0,L)$$

$$u(x,t)=\sum _{-\infty }^{\infty }{b}_k(t)e^{ikqx},\text{where}q=\frac{2\pi }{L}.$$

$${\dot{b}}_k=(k^2-\nu k^4)b_k+ik\sum _{m=-\infty }^{\infty }{b}_m{b}_{k-m}$$

$${\dot{a}}_k=(k^2-\nu k^4)a_k-k\sum _{m=-\infty }^{\infty }{a}_m{a}_{k-m}.$$

$$f(x)=x^4-x^2.$$

$$F_i(\mathbf{z})=\sum _{k=0}^K{a}_{ik}{\varphi }_k(\mathbf{z}),$$

$${\hat{F}}_i(\hat{\mathbf{z}}(t))=\sum _{k=0}^K{a}_{ik}{\varphi }_k(\hat{\mathbf{z}}(t))+{\xi }_i(t),(t=1,\dots ,T;i=1,\dots ,n).$$

$${\mathbf{f}}^{(i)}=\Phi {\mathbf{a}}^{(i)}+{\mathbf{\xi }}^{(i)},$$

$$\mathbf{f}=\Phi \mathbf{a}+\mathbf{\xi },$$

$$\underset{\mathbf{a}\in {\mathbb{R}}^K}{\min }\Vert \Phi \mathbf{a}-\mathbf{f}{\Vert }_2.$$

$${\mathbf{a}}_{\text{mbox}(LS)}={\Phi }^{†}\mathbf{f},$$

$$\{\begin{array}{l}(k_1,a_{k_1})={\mathrm{argmin}}_{k,c}\Vert \mathbf{f}-c{\mathbf{\varphi }}_k{\Vert }_2,\\ {\mathbf{r}}_1=\mathbf{f}-{\mathbf{\varphi }}_{k_1}{a}_{k_1}.\end{array}$$

$$\{\begin{array}{l}(k_{\ell +1},a_{k_{\ell +1}})={\mathrm{argmin}}_{k,c}\Vert {\mathbf{r}}_{\ell }-c{\mathbf{\varphi }}_k{\Vert }_2,\\ {\mathbf{r}}_{\ell +1}={\mathbf{r}}_{\ell }-{\mathbf{\varphi }}_{k_{\ell +1}}{a}_{k_{\ell +1}}.\end{array}$$

$$\underset{\mathbf{a}\in {\mathbb{R}}^K}{\min }\left(\Vert \Phi \mathbf{a}-\mathbf{f}{\Vert }_2^2+\lambda \Vert \mathbf{a}{\Vert }_1\right),$$

$$\{\begin{array}{l}{\mathrm{argmin}}_{\mathbf{a}}\Vert \mathbf{a}{\Vert }_1,\\ \text{mbox}subjectto\Vert \Phi \mathbf{a}-\mathbf{f}\Vert \le ϵ,\end{array}$$

$$\displaystyle\left(\begin{array}[]{ccc}6&3&2\\ 2&1&1\\ 1&2&1\end{array}\right)\bm{a}=\left(\begin{array}[]{c}6\\ 2\\ 4\end{array}\right).$$

$$\displaystyle\left(\begin{array}[]{ccc}6&3&2\\ 2&1&1\end{array}\right)\bm{a}=\left(\begin{array}[]{c}6\\ 2\end{array}\right)$$

Full text

Keywords: System Identification, Sparse Regression, Data-Driven Modeling, Entropy, Conditional Mutual Information, Asymptotic Equipartition Property, Nonlinear dynamics, Complex Systems.

System identification (SID) is a central concept in science and engineering applications whereby a general model form is assumed, but active terms and parameters must be inferred from observations. Most methods for SID rely on optimizing some metric-based cost function that describes how a model fits observational data. A commonly used cost function employs a Euclidean metric and leads to a least squares estimate, whereas recently it has become popular to also account for model sparsity such as in compressed sensing and Lasso. While the effectiveness of these methods has been demonstrated in previous studies including in cases where outliers exist in sparse samples, SID remains particularly difficult under more realistic scenarios where each observation is subject to non-negligible noise, and sometimes even contaminated by large noise outliers. Here we report that existing sparsity-focused methods such as compressive sensing, when applied in such scenarios, can result in “over sparse” solutions that are brittle to outliers. In fact, metric-based methods are prone to outliers because outliers by nature have an unproportionally large influence. To mitigate such issues of large noise and outliers, we develop an entropic regression approach for nonlinear SID, whereby true model structures are identified based on an information theoretic criterion describing relevance in terms of reducing information flow uncertainty, versus not necessarily (just) sparsity. The use of information-theoretic measures in entropic regression poses unique advantages, due to the asymptotic equipartition property of probability distributions, that outliers and other low-occurrence events are conveniently and intrinsically de-emphasized as not-typical, by definition.

A basic and fundamental problem in science and engineering is to collect data as observations from an experiment, and then to attempt to explain the experiment by summarizing data in terms of a model. When dealing with a dynamical process, a common scenario is to describe the underlying process as a dynamical system, which may be in the form of a differential equation (DE). Traditionally this means “understanding the underlying physics,” in a manner that allows one to write a DE from first principles, including those terms to capture the delicate but important (physical) effects. Validation of the model may come from comparing outputs from the model to those from experiments, where outputs are typically represented as multivariate time-series. Building a DE model based on fundamental laws and principles requires strong assumptions, which might be evaluated by how the model fits data. Weigenbend and Gershenfeld made a distinction between weak modeling (data rich and theory poor) and strong modeling (data poor and theory rich), and suggest that it is related to “…the distinction between memorization and generalization…” [15].

The problem of learning a (dynamical) system from observational data is known as system identification (SID), and often times involves the underlying assumption that the structural form of the DE is known (which kinds of terms to include in the functional description of the equation), but only the underlying parameters are not known. For example, suppose we observe the dynamics of a simple dissipative linear spring, then we may express the model as $m\ddot{x}+\gamma\dot{x}+kx=0$ based on Hooke’s law. However, the parameters $m,\gamma,$ and $k$ might be unknown and need to be estimated in order to completely specify the model for purposes such as prediction and control. One may directly measure those parameters by static testing (e.g., weighing the mass on a scale). Alternatively, here we are interested in utilizing the observational data generated by the system without having to design and perform additional experiments, to estimate the parameters corresponding to the model that best fits empirical observations, which is a standard viewpoint in SID. In this thought experiment, the SID process is performed with the underlying physics understood (the form of the Hooke spring equation). In general it can be applied in the scenario where very little information is previously known about the system, in a black box manner.

Suppose that observations $\{\bm{z}(t)\}$ come from a general (multidimensional, coupled) DE, represented by

[TABLE]

where $\bm{z}=[z_{1},\dots,z_{N}]^{T}\in\mathbb{R}^{N}$ is the (multivariate) state variable of the system and $\bm{F}=[F_{1},\dots,F_{N}]^{\top}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ is the vector field. Each component function $F_{i}(\bm{z})$ can be represented using a series expansion (for example a power series or a Fourier series), writing generally,

[TABLE]

for a linear combination of basis functions $\{\phi_{k}\}_{k=0}^{\infty}$. The basis functions do not need to be mutually orthogonal, and the series can even include multiple bases, for example to contain both a polynomial basis and a Fourier basis [5]. The coefficients $\{a_{ik}\}$ are to be determined by contrasting simulations to experimental measurements, in an optimization process whose details of how error is measured distinguishes the various methods we discuss here. This was the main theme in previous approaches on nonlinear SID, with different methods differ mainly on how a model’s fit is quantified. The different approaches include using standard squared error measures [53, 9], sparsity-promoting methods [26, 5, 52, 51] as well as using entropy-based cost functions [19]. Among those, sparsity-promoting methods have proven particularly useful because they tend to avoid the issue of overfitting, thus allowing a large number of basis functions to be included to capture possibly rich dynamical behavior [26, 26, 5, 52].

Regardless of the particular method or system, most previous work on nonlinear SID focused on the low-noise regime and demonstrated success only when there is a sufficient amount of clean observational data. In practice, an observation process can be subject to external disturbances in unpredictable ways. Consequently, the effective noise can be quite large and even with frequently occurring “outliers” both of which may contaminate the otherwise perfect data. Can SID still work under the presence of large noise and outliers? At a glance, the answer should be yes, given that several recent SID methods for nonlinear systems are readily deployable in the presence of noise. For example, compressive sensing can handle noise by relaxing the constraint set whereas least squares and Lasso can be applied off the shelf—the important question however is whether the quality of solution is compromised or not, and to what extent. Recently Tran and Ward considered the nonlinear SID problem under the presence of outliers in observational data and showed that so long as there the outliers are “sparse” leaving sufficient amount of “clean” data available, existing techniques such as SINDy can be extended to reconstruct the exact form of a system with high probability [49]. In the current work, we are interested in the more realistic scenario where effective noise is present everywhere and thus all data points are contaminated by non-negligible noise and sometimes outliers. These features effectively creates a “high noise and low data amount” regime, where we found that existing nonlinear SID methods including recent ones that specialize in promoting sparsity, fall short.

In this work we depart from most standard approaches for nonlinear SID. We identify the error quantification via metric-based cost functions as a root cause of existing methods to fail under large noise and outliers because outliers tend to deviate from the rest of sample data as measured by metric distance; thus trying to “fit” the outliers almost inevitably causes the model to put (much) less weights on the “good” data points. To resolve this important issue, we propose to infer the (sparsity) structure of a general model together with its parameters using a novel information theoretic regression approach that we call Entropic Regression (ER). As we will show, while standard metric-based methods emphasize the data in ways as designed by the chosen metric, the proposed ER approach is robust with regards to the presence of noise and outliers in the data. Instead of searching for the sparsest model and thus risk forcing a wrong sparse model, ER is emphasizing “information relevance” according to a model-free, entropic criterion. Basis terms will be included in the model only because they are relevant and not (necessarily) because they together make up the sparsest model. We demonstrate the effectiveness of ER in several examples, including chaotic Lorenz systems, Kuramoto-Sivashinsky equations, and a double well potential, where in each case the observed data contains relatively large noise and outliers. We also remark on the computational complexity and convergence in small-data regime, as well as discuss open problems and future directions.

Results

Nonlinear System Identification: Problem Statement and Formulation

Following the standard routine in nonlinear SID [34], the starting point is to recast the nonlinear SID problem into a computational inverse problem, by considering an appropriate set of basis functions that span the space of functions including the system of interest [53, 51]. A common choice is the standard polynomial basis

[TABLE]

where each term is a monomial. Using a set of basis functions, one can represent the individual component functions of $F$ as a series as in (2). The specification of the location of nonzero parameters are referred to as the structure of the model.

Consider time series data $\{\bm{z}(t)=[z_{1}(t),\dots,z_{m}(t)]^{\top}\}_{t=t_{0},\dots,t_{\ell}}$ and corresponding $\{\bm{F}(\bm{z}(t)\}_{t=t_{0},\dots,t_{\ell}}$ generated from a nonlinear, high-dimensional dynamical system (1), possibly subject to observational noise. From $\bm{z}(t)$, one can estimate the derivatives by any of the standard Newton-Cotes methods, explicit Euler’s method of course being the simplest, giving $F_{i}(\bm{z}(t_{k}))=\frac{z_{i}(t_{k+1})-z_{i}(t_{k})}{\tau_{k}}+\mathcal{O}((t_{k+1}-t_{k}))$, or central difference which has improved accuracy: $F_{i}(\bm{z}(t_{k}))=\frac{z_{i}(t_{k+1})-z_{i}(t_{k-1})}{t_{k+1}-t_{k-1}}+\mathcal{O}((t_{k+1}-t_{k-1})^{2})$. The problem of nonlinear system identification is to reconstruct the functional form as well as parameters of the underlying system, that is, to infer the nonlinear function $\bm{F}$.

Under the basis representation (2), the identification of $\bm{F}$ becomes equivalent to estimating all the parameters $\{a_{ik}\}$. In practice, the empirically observed state variable is subject to noise: $\hat{\bm{z}}(t)=\bm{z}(t)+\bm{\eta}(t)$ with $\bm{\eta}(t)$ representing the (multivariate) noise and $\hat{F}_{i}$ denoting the approximated value of $F_{i}$. For noisy observations $\hat{\bm{z}}(t)$, the difference between $\hat{F}_{i}(\hat{\bm{z}}(t))$ and $F_{i}(\hat{\bm{z}}(t))$ originates from several sources: the infinite series is truncated and the derivatives are estimated numerically and by using approximate states. Nevertheless, we can represent the aggregated error as an effective noise $\bm{\xi}(t)$ term and express the forward model as

[TABLE]

Note that because of the combined and accumulated effects of observational noise, approximation error and truncation, even if the observational noise of the states $\eta_{i}(t)$ are iid, this is not necessarily true for the effective noise $\xi_{i}(t)$. In matrix form, the forward model (4) has the approximate expression

[TABLE]

Figure 1 shows the structure of the Lorenz system under standard polynomial basis up to quadratic terms.

In vector form, under a choice of basis and truncation, the nonlinear system identification problem can be recast into the form of a linear inverse problem

[TABLE]

where $\bm{f}^{(i)}=[\hat{F}_{i}(\hat{\bm{z}}(t_{1})),\dots,\hat{F}_{i}(\hat{\bm{z}}(t_{\ell}))]^{\top}\in\mathbb{R}^{\ell\times 1}$ represents the $i$-th component of the estimated vector field from the observational data, $\Phi=[\bm{\phi}^{(1)},\dots,\bm{\phi}^{(K)}]\in\mathbb{R}^{\ell\times K}$ (with $\bm{\phi}^{(k)}=[\phi_{k}(\hat{\bm{z}}(t_{1})),\dots,\phi_{k}(\hat{\bm{z}}(t_{\ell}))]\in\mathbb{R}^{\ell\times 1}$) represent sampled data for the basis functions, $\bm{\xi}^{(i)}=[\xi_{i}(t_{1}),\dots,\xi_{i}(t_{\ell})]^{\top}\in\mathbb{R}^{\ell\times 1}$ represents noise, and $\bm{a}^{(i)}=[a_{i1},\dots,a_{iK}]^{\top}\in\mathbb{R}^{K\times 1}$ is the vector of parameters which is to be determined. Note that the form of the equation (16) is the same for each $i$, and solving each $\bm{a}^{(i)}$ can be done separately and independently for each $i$. In what follows we omit the index when discussing the general methodology, and consider the following linear inverse problem

[TABLE]

where $\bm{f}\in\mathbb{R}^{\ell\times 1}$ and $\Phi\in\mathbb{R}^{\ell\times K}$ are given, with the goal to estimate $\bm{a}\in\mathbb{R}^{K\times 1}$. This general problem is in the form of an inverse problem and is typically solved under various assumptions of noise by methods such as least squares, orthogonal least squares, lasso, compressed sensing, to name a few. Each of these methods, in addition to the recent approach of SINDy and its generalization, is mentioned in the Results section and reviewed in the Methods section. In what follows we develop a unique information-theoretic approach called entropic regression, which we demonstrate has significant advantages.

Entropic Regression

To overcome the competing challenges of potential overfitting, efficiency when limited data points are given, and robustness with respect to noise and in particular outliers in observations, we propose a novel framework that combines the advantage of information-theoretic measures and iterative regression methods. The framework, which we term entropic regression (ER), is model-free, noise-resilient, and efficient in discovering a “minimally sufficient” model to represent data. The key idea is that, for given set of basis functions, a model should be considered minimally sufficient if no basis function that is not already included in the model can help increase the information relevance between the model outputs and observed data. In other words, the residual between the model fit and observational data is statistically independent from any basis function that is not included in the model—because otherwise the dependence can be harvested to reduce the discrepancy by including such a basis function in the model. We emphasize that, although the idea seems related to classical model selection principles such as AIC [1], ours combines model construction with selection. In addition, even though it is not uncommon for entropy measures to be adopted in system identification [19, 37], the proposed method is unique as it fuses entropy optimization with regression in a principled manner that enables scalable computation and efficient estimation in reconstruction nonlinear dynamics under noisy data. As we shall see below, the proposed ER method is applicable even in the small-sampling regime (by adopting appropriately defined entropy measures and efficient estimators) and naturally allows for a computationally efficient procedure to build up a model from scratch. In particular, we use (conditional) mutual information as an information-theoretic criterion and iteratively select relevant basis functions, analogous to the optimal causation entropy algorithm previously developed for causal network inference [47, 48] but here including an additional regression component in each step. Thus, ER can be thought of as an information-theoretic extension of the orthogonal least squares regression, or as a regression version of optimal causation entropy.

We now present the details of ER. The ER method contains two stages (also see Algorithm 1 for the pseudocode): forward ER and backward ER. In both stages, selection and elimination are based on an entropy criterion and parameters are updated in each iteration using a standard regression (e.g., least squares). Consider the inverse problem (17). For an index set $S\subset\mathbb{N}\cup\{0\}$, the estimated parameters can be thought of as a mapping from the joint space of $\Phi$, $\bm{f}$ and $S$ to a vector denoted as $\hat{\bm{a}}=R(\Phi,\bm{f},S)$. For instance, under a least-squares criterion the mapping is given by $R(\Phi,\bm{f},S)_{S}=\Phi_{S}^{\dagger}\bm{f}$ ($\Phi_{S}$ denotes the columns of matrix $\Phi$ indexed by $S$) and $R(\Phi,\bm{f},S)_{i}=0$ for all $i\notin S$. Using the estimated parameters, the recovered signal can be computed as $\Phi R(\Phi,\bm{f},S)$. In the ER algorithm, we start by selecting a basis function $\bm{\phi}_{k_{1}}$ that maximizes its mutual information with $\bm{f}$, compute the corresponding parameter $a_{k_{1}}$ using the least squares method, and obtain the corresponding regression model output $\bm{z}_{1}$ according to

[TABLE]

Here $I(\bm{x};\bm{y})$ denotes mutual information between $\bm{x}$ and $\bm{y}$, which is a model-free measure of the statistical dependence between two distributions (that is, $\bm{x}$ and $\bm{y}$ are independent if and only if their mutual information equals zero) [12]. Next, in each iteration of the forward stage, we perform the following computations and updates, for $i=2,3,\dots$,

[TABLE]

The process terminates when $\max_{k}I(\Phi R(\Phi,\bm{f},k);\bm{f}|\bm{z}_{i-1})\approx 0$ (or when all basis functions are exhausted), indicating that none of the remaining basis function is relevant given the current model, in an information-theoretic sense. The result of the forward ER is a set of indices $S=\{k_{1},\dots,k_{m}\}$ together with the corresponding parameters $a_{k_{1}},\dots,a_{k_{m}}$ ($a_{j}=0$ for $j\notin S$) and model $f\approx a_{k_{1}}\bm{\phi}_{k_{1}}+\dots+a_{k_{i}}\bm{\phi}_{k_{i}}.$ Finally, we turn to the backward stage, where the terms that had previously been included are re-examined for their information-theoretic relevance and these that are redundant will be removed. In particular, we sequentially check for each $j=k_{i}\in S$ to determine if the basis term $\bm{\phi}_{j}$ is redundant by computing

[TABLE]

and updating $S\rightarrow S/\{j\}$ (that is, remove $j$ from the set $S$) if $I(\Phi R(\Phi,\bm{f},S);\bm{f}|\bm{\bar{z}}_{j})\approx 0$. The result of the backward ER is the reduced set of indices $S=\{\ell_{1},\dots,\ell_{n}\}$ with $n\leq m$, together with the corresponding parameters $a_{\ell_{1}},\dots,a_{\ell_{n}}$ ($a_{j}=0$ for $j\not\in S$) computed as $\bm{a}=R(\Phi,\bm{f},S)$, and accordingly the recovered model $\bm{f}\approx\bm{\phi}\bm{a}=\bm{\phi}_{S}\bm{a}_{S}=a_{\ell_{1}}\bm{\phi}_{\ell_{1}}+\dots+a_{\ell_{n}}\bm{\phi}_{\ell_{n}}$. In practice, mutual information and conditional mutual information need to be estimated from data, and whether or not the estimated values should be regarded as zero is typically done via (approximate) significance testing, the details of which are provided in Methods (also see Sec. SI.3).

Numerical Experiments: Outliers, Expansion Order, and the Paradox of Sparsity

To demonstrate the utility of ER for nonlinear system identification under noisy observations, we benchmark its performance against existing methods including least squares (LS), orthogonal least squares (OLS), Lasso, as well as SINDy and its extension by Tran and Ward (TW). The details of the existing approaches are described in the Methods Section. The examples we consider represent different types of systems and scenarios, including both ODEs and PDEs. In addition, we consider different noise models and especially the presence of outliers in order to evaluate the robustness of the respective methods.

For each example system, we sample the state of each variable at a uniform rate of $\Delta{t}$ to obtain a multivariate time series $\{\bm{z}(t_{i})\}_{k=1,\dots,N;i=1,\dots,\ell}$ where $\bm{z}=[z_{1},\dots,z_{d}]^{\top}\in\mathbb{R}^{d}$; then we add noise to each state variable and obtain the noisy empirical time series denoted by $\{\hat{\bm{z}}(t_{i})\}$, where

[TABLE]

with $\eta_{ki}$ representing state observational noise. The vector field $\bm{F}$ is estimated using central difference on the noisy time series $\{\hat{\bm{z}}(t)\}$.

Example 1. Chaotic Lorenz system. Our first detailed example data set was generated by noisy observations from a chaotic Lorenz system, which is represented by a three-dimensional ODE which is a prototype system as a minimal model for thermal convection obtained by a low-ordered modal truncation of the Saltzman PDE [41], and for many parameter combinations exhibits chaotic behavior [35]. In our standard notation, we have $\bm{z}=[z_{1},z_{2},z_{3}]^{\top}$ and

[TABLE]

with default parameter values $\sigma=10$, $\rho=28$ and $\beta=8/3$ unless otherwise specified. We consider a standard polynomial basis as in Eq. (3). Over recent years, the Lorenz system has become a favorable and standard example for testing SID methods and typically requires tens of thousands of measurements for accurate reconstruction [49, 5].

First, we compare several nonlinear SID methods in reconstructing the Lorenz system when the state observational noise is drawn independently from a Gaussian distribution, $\eta\sim\mathcal{N}(0,\epsilon^{2})$. As we discussed before, this translates into effective noise that is not necessarily Gaussian or even independent. Fig. 2 shows the error in the estimated parameters where, $error=\|\bm{a}_{true}-\bm{a}_{estimated}\|_{2}$. As shown in Fig. 2, even with observational noise as low as $\epsilon=10^{-4}$, ER and OLS outperform all other methods. In this low noise regime, SINDy required more measurements (around 4 times) to reach similar accuracy as ER. In comparison, as noted in [49, 5] and in the implementation provided by the authors, for SINDy and TW methods to yield accurate reconstruction the number of measurements is at the order of $10^{4}$.

Next, to explore the performance of SID methods under the presence of outliers, we conduct additional numerical experiments. The extent to which outliers present is controlled by a single parameter $p$: each observation is subject to an added noise $\eta$, where $\eta\sim\mathcal{N}(0,\epsilon_{1}^{2})$ with probability $1-p$ and $\eta\sim\mathcal{N}(0,\epsilon_{1}^{2}+\epsilon_{2}^{2})$ with probability $p$. Here we use $\epsilon_{1}=10^{-5}$, $\epsilon_{2}=0.2$ and $p=0.2$. The results of SID are shown in Fig. 3. Compared to Fig.(2), we see that with $p>0$ OLS performance drops due to the increasing occurrence of large noise and outliers whereas ER remains its capacity of accurately identifying the underlying system. As an example, in each of the side panels of Fig. 3, we show the trajectory of the identified dynamics using the median solution of each method. It is clear that under such noisy chaotic dynamics and at a relatively under-sampled regime, ER method successfully recovers the system dynamic. As an ample amount of data becomes available, we note that TW method starts to produce excellent reconstruction which is consistent with recent findings reported in Ref. [49].

Given that a major theme of modern SID is to seek for sparse representations, and the Lorenz system under standard polynomial basis is indeed sparse, it is worth asking: what are the respective structure identified by the different methods? In Fig. 4 we compare the structure of the identified model using different methods across a range of parameter values for $\rho$. In this case, under the presence of large noise and outliers ($p=0.2$), none of the methods examined here, including recently proposed sparsity-promoting (CS, SINDy) and outlier-resilient (TW) methods, is able to identify the correct structure. The proposed ER method, however, does identify the correct structure. It is worth pointing out that, often times when expressed in the right basis, a model will appears to be sparse, the converse is not true: just because a method return a sparse solution does not suggest (at all) the such a solution gives a reasonable approximation of the true model structure. Interestingly, as we show in Fig. SI.1 and Fig. SI.2, for the same system and data, as more basis functions are used–that is, when the true dynamics becomes sparser–the reconstructed dynamics using existing methods (such as CS) can become worse.

Example 2. Kuramoto-Sivashinsky equations. To further demonstrate the power of ER, we consider a nonlinear PDE, namely the Kuramoto-Sivashinsky (KS) equation [32, 31, 45, 22, 33], which arises as a description of flame front flutter of gas burning in a cylindrically symmetric burner. It has become a popular example of a PDE that exhibits chaotic behavior, in particular spatiotemporal chaos [11, 21]. We will consider Kuramoto-Sivashinsky system in the following form,

[TABLE]

in periodic domain, $u(t,x)=u(t,x+L)$, and we restrict our solution to the subspace of odd solutions $u(t,-x)=-u(t,x)$. The viscosity parameter $\nu$ controls the suppression of solutions with fast spatial variations, and is set to $\nu=0.029910$ under which the system exhibit chaotic behavior [11].

Since a PDE corresponds to an infinite-dimensional dynamical system, in practice we focus on an approximate finite-dimensional representation of the system, for example, by Galerkin-projection onto basis functions as infinitely many ODE’s in the corresponding Banach space.

To develop the Galerkin projection, we follow the procedure as presented in [13], to expand a periodic solution $u(x,t)$ using a discrete spatial Fourier series,

[TABLE]

Notice that we have written this Fourier series of basis elements $e^{ikqx}$ in terms of time varying combinations of basis elements. For simplicity, consider $L=2\pi$, then $q=1$ for the following analysis. This is typical [40] with the representation of a PDE as infinitely many ODE’s in the Banach space, where orbits of these coefficients therefore become time varying patterns by Eq. (23). Substituting Eq. (23) into Eq. (22), we produce the infinitely many evolution equations for the Fourier coefficients,

[TABLE]

In general, the coefficients $b_{k}$ are complex functions of time $t$. However, by symmetry, we can reduce to a subspace by considering the special symmetry case that $b_{k}$ is pure imaginary, $b_{k}=ia_{k}$ and $a_{k}\in{\mathbb{R}}$. Then,

[TABLE]

where $k=1,..,N_{m}$. However, the assumption that there is a slow manifold (slow modes as an inertial manifold [40, 24, 38, 25]) suggests the practical matter that a finite truncation of the series Eq. (23), and correspondingly the a reduction to finitely many ODEs will suffice. Therefore we choose a sufficiently large number of modes $N_{m}$. Then we solve the resulting $N_{m}$-dimensional ODE (25) to produce the estimated solution of $u(x,t)$ by (23), and use such data for the purpose of SID, have meaning to estimate the structure and parameters of the ODE model (25).

Fig. 5 shows the first three dimensions plot under different number of modes. We see that using just a few number of modes ($N_{m}=8,...,11$) is insufficient to capture the true dynamical behavior of the system whereas too large a number of modes ($N_{m}=20,24$) may be unnecessary. In this example, an adequate but not excessive number of modes seems to be around $N_{m}=16$, as no significant information is gained by increasing $N_{m}$.

Fig. 6 shows the sparse structure of the recovered solution by different methods. Here we mention that the true non-zero parameters of KSE using $N_{m}=16$ are 200 parameters that vary in the magnitude from 0.9701 to 1705. With the second order expansion, our basis matrix will have 153 candidate functions, and it will be nearly singular with condition number $4\times 10^{7}$. Likely due to such high condition number, neither TW nor SINDy gives reasonable reconstruction. In particular, we note that the solution of SINDy is already optimized by selecting the threshold value $\lambda$ that is slightly above $\lambda_{*}$ where here $\lambda_{*}\approx 0.1731$ is the smallest magnitude of the true nonzero parameter of the full least squares solution. A larger value of $\lambda$ only worsens the reconstruction, as we found numerically.

The OLS method overcomes the disadvantage of LS by iteratively finding the most relevant “feature” variables, where relevance is measured in terms of (squared) model error; but it comes at a price: similar to LS, the OLS is sensitive to outliers in the data and such sensitivity seems to be even more amplified due to the smaller number of terms typically included in OLS as compared to LS, which cause the high false negative rate in the OLS solution. Although ER solution has few false negatives, but was completely able to recover the overall dynamic of the system as shown in Fig.( 7), while all other solutions diverges and failed to recover $u(x,t)$.

Example 3. Double Well Potential. Finally, in order to gain further insights into why standard methods fail under the presence of outliers, we consider a relatively simple double-well system, with

[TABLE]

Suppose that we measure $x$ and $f$, can we identify the function $f(x)$? We sample 61 equally spaced measurements for $x\in[-1.2,1.2]$, and we construct $\Phi$ using the $10^{th}$ order polynomial expansion with $K=11$ is the number of candidate functions. Then, we consider a single fixed value corrupted measurement to be $f(0.6)=0.2$.

Fig. 8 shows the results the double-well SID under a single outlier in the observation. We see the robustness of ER solution to the outliers while CS failed in detecting the system sparse structure. For the sake of clearness, Fig. 8 shows the results for CS and ER. The results for each solver and details are provided in Sec. (SI.4.1) in addition to more numerical examples.

Discussion

The main theme of the paper is on nonlinear system identification (SID) under noisy observations, which is to learn the functional form and parameters of a nonlinear system based on observations of its states under the presence of noise and outliers. We recast the problem into the form of an inverse problem using a basis expansion of the nonlinear functions. Such basis expansion, however, renders the resulting problem inherently high dimensional even for low-dimensional systems. In practice, the need for finite-order truncation as well as the presence of noise causes additional challenges. For instance, even under iid Gaussian observational noise for the state variables, the effective noise in the inverse problem is not necessarily so. As we demonstrate using several example systems, including the chaotic Lorenz system and the Kuramoto-Sivashinsky equations, existing SID methods are prone to noise, and can be quite sensitive to the presence of outliers. We identify the root cause of such non-robustness being the metric nature of the existing methods, as they quantify error based on metric distance, and thus a handful of data points that are “corrupted” by large noise can dominate the model fit. Each of the existing methods we considered has this property, which includes the least squares, compressive sensing, and Lasso. From a mathematical point of view, each method can be interpreted as a functional that maps input data to a model, through some optimization process. In a noisy setting, the output model should ideally change smoothly with respect to the input data, not just continuously. Our results suggest that these popular methods in fact do suffer from a sensitive dependence on outliers, as a few corrupted data can already produce very poor model estimates. Alarmingly, the now-popular CS method, which is based on sparse regression, can force to select a completely wrong sparse model under noisy input data, and this occurs even when there is just a single outlier. This is by no means contradicting previous findings of the success of CS in SID, as in such work noise is typically very small, and here we are considering a perhaps more realistic scenario with larger noise.

To fill the vacancy of SID methods that can overcome outliers, we develop an information-theoretic regression technique, called entropic regression (ER), that combines entropy measures with an iterative optimization for nonlinear SID. We show that ER is robust to noise and outliers, in the otherwise very challenging circumstances of finding a model that explains data from dynamical stochastic processes. The key to ER’s success is its ability to recover the correct and true sparsity structure of a nonlinear system under basis expansions, despite either relatively large noise, or alternatively even relatively many even larger outliers. In this sense ER is superior to any other method that we know of for such settings. Note that in the ER algorithm, least squares is used to estimate the parameters of those basis functions that are deemed relevant where relevance is detected using an information-theoretic measure that is insensitive to noise and outliers. The choice of least squares in the regression step in ER is not necessarily an optimal choice and can be potentially replaced by more advanced methods (e.g., those developed in robust regression). In the current implementation of ER we adopted least squares mainly due to its computational advantage over alternative methods. On a more fundamental level, ER’s robustness against outliers may likely be attributed to an important principle in information theory called the asymptotic equipartition property (AEP) [12]. The outcome of this principle is that sampled data can be partitioned into “typical” samples and “atypical” samples, with the rare atypical samples end up influencing the estimated entropy relatively weakly. Since ER measures relevance by entropy instead of metric distance, a few outliers, no matter how far away they are from the rest of the data points, tend to have minimal impact on the model identification process. So the general interpretion we make here is that outliers observations are likely atypical, but not part of the core of data that carry the major estimation of the entropy. This foundational concept of information theory is likely the major source of robustness of our ER method to system identification.

Methods

Existing metric-based methods for system identification

Recall (from the main text) that we recast the nonlinear system identification problem here. Given a truncated basis representation of each component of the vector field $\bm{F}$, expressed as

[TABLE]

we consider sampled data $\hat{\bm{z}}$ and the estimated vector field $\hat{\bm{F}}$, from which the coefficients (parameters) $\{a_{ik}\}$ are to be determined. In general, we use subscript “$t$” to index the sampled data, and thus the $t$-th sample satisfies the equation

[TABLE]

Here $\xi_{i}(t)$ is the effective noise that represents the accumulative impact of truncation error, state observational noise as well as approximation error in the estimation of derivatives. Consequently, an iid Gaussian noise additive to the states $\bm{z}_{i}(t)$ can translate into correlated non-Gaussian effective noise for $\xi_{i}(t)$.

Having transformed a system identification problem into an parameter estimation problem (or inverse problem) in the form of

[TABLE]

where $\bm{f}^{(i)}=[\hat{F}_{i}(\hat{\bm{z}}(1)),\dots,\hat{F}_{i}(\hat{\bm{z}}(T))]^{\top}\in\mathbb{R}^{T\times 1}$ represents the estimated function $F_{i}$ ($i$-th component of the vector field $\bm{F}$), $\Phi=[\bm{\phi}^{(1)},\dots,\bm{\phi}^{(K)}]\in\mathbb{R}^{T\times K}$ (with $\bm{\phi}^{(k)}=[\phi_{k}(\hat{\bm{z}}(1)),\dots,\phi_{k}(\hat{\bm{z}}(T))]\in\mathbb{R}^{T\times 1}$) represent sampled data for the basis functions, $\bm{\xi}^{(i)}=[\xi_{i}(1),\dots,\xi_{i}(T)]^{\top}\in\mathbb{R}^{T\times 1}$ represents effective noise, and $\bm{a}^{(i)}=[a_{i1},\dots,a_{iK}]^{\top}\in\mathbb{R}^{K\times 1}$ is the vector of parameters which is to be determined. Since the form of the equation (29) is the same for each $i$, we omit the index when discussing the general methodology, and consider the following linear inverse problem

[TABLE]

where $\bm{f}\in\mathbb{R}^{T\times 1}$ and $\Phi\in\mathbb{R}^{T\times K}$ are given, with the goal is to estimate $\bm{a}\in\mathbb{R}^{K\times 1}$ when the effective noise is not necessarily from independent multivariate Gaussian distribution.

Least Squares (LS)

The most commonly used approach to estimate $\bm{a}$ in Eq. (30) is to use the least squares criterion, which finds $\bm{a}$ by solving the following least squares minimization problem:

[TABLE]

The solution can be explicitly computed, giving

[TABLE]

where $\Phi^{\dagger}$ denotes the pseudoinverse of the matrix $\Phi$ [16]. Note that in the special case where the minimum is zero (which is unlikely under the presence of noise), the minimizer is not unique and the “least-squares” solution typically refers to a vector $\bm{a}$ that has the minimal $2$-norm and solves the equation $\Phi\bm{a}=\bm{f}$. The LS method has several advantages: it is analytically traceable and easy to solve computationally using standard linear algebra routines (e.g., SVD). However, a main disadvantage of the LS approach in system identification, as we discuss in the main text, is that it generally produces a “dense” solution, where most (if not all) components of $\bm{a}$ are nonzero, which is a severe overfitting of the actual model. This (undesired) feature also makes the method sensitive to noise, especially in the under-sampling regime.

Orthogonal Least Squares (OLS)

In orthogonal least squares (OLS) [9, 50, 29], the idea is to iteratively select the columns of $\Phi$ that minimize the (2-norm) model error, which corresponds to iterative assigning nonzero values to the components of $\bm{a}$. In particular, the first step is to select basis $\bm{\phi}_{k_{1}}$ and compute the corresponding parameter $a_{k_{1}}$ and residual $\bm{r}_{1}$ according to

[TABLE]

Then, one iteratively selects additional basis functions (until stopping cretia is met) and compute the corresponding parameter value and residual, as

[TABLE]

As for stopping criteria, there are several choices including AIC and BIC. In this work, in the absence of knowledge of the error distribution, we adopt a commonly used criterion where the iterations terminate when the norm of the residual is below a prescribed threshold. To determine the threshold, we consider 50 log-spaced candidate values in the interval $[10^{-6},100]$ and select the best using 5-fold cross validation.

Lasso

A principled way to impose sparsity on the model structure is to explicitly penalize solution vectors that are non-sparse, by formulating a regularized optimization problem:

[TABLE]

where the parameter $\lambda\geq 0$ controls the extent to which sparsity is desired: as $\lambda\rightarrow\infty$ the second term dominates and the only solution is a vector of all zeros, whereas at the other extreme $\lambda=0$ and the problem becomes identical to a least squares problem which generally yields a full (non-sparse) solution. Values of $\lambda$ in between then balances the “model fit” quantified by the 2-norm and the sparsity of the solution characterized by the 1-norm. For a given problem, the parameter $\lambda$ needs to be tuned in order to specify a particular solution. A common way to select $\lambda$ is via cross validation [20]. In our numerical experiments, we choose $\lambda$ span according to  [20], with 5-Folds cross validation and 10 values $\lambda$ span. We adopt the CVX solver [18], and from all the solutions found for each $\lambda$ we select the solution with minimum residual.

Compressed sensing (CS)

Originally developed in the signal processing literature [8, 7, 14], the idea of compressed sensing (CS) has been adopted in several recent work in nonlinear system identification [52, 51] Under the CS framework, one solves the following constrained optimization problem,

[TABLE]

where the parameter $\epsilon\geq 0$ is used to relax the otherwise strict constraint $\Phi\bm{a}=\bm{f}$, to allow for the presence of noise in data. In our numerical experiments, we choose 10 log-spaced values for $\epsilon\in[10^{-6},100]$, and 5-Folds cross validation. We adopt the CVX solver [18], and from all the solutions found for each $\epsilon$ we select the solution with minimum residual.

SINDy

In their recent contribution, Brunton, Proctor and Kutz introduced SINDy (Sparse Identification of Nonlinear Dynamics) as a way to perform nonlinear system identification [5]. Their main idea is, after formulating the inverse problem (30), to seek a sparse solution. In particular, given that Lasso can be computationally costly, they proposed to use sequential least squares with (hard) thresholding as an alternative. For a (prechosen) threshold $\lambda$, the method starts from a least squares solution and abandons all basis functions whose corresponding parameter in the solution has absolute value smaller than $\lambda$; then the same is repeated for the data matrix associated with the remaining basis functions, and so on and so forth, until no more basis function (and the corresponding parameter) is removed. For fairness of comparison, we present results of SINDy according to the best threshold parameter $\lambda$ manually chosen so that no active basis function is removed in the very first step (see KSE example); for the Lorenz system example, we choose $\lambda=0.02$ as used in a similar example as in Ref. [5].

Tran-Ward (TW)

In their recent paper [49] Tran and Ward considered the SID problem where certain fraction of data points are corrupted, and proposed a method to simultaneously identify these corrupted data and reconstruct the system assuming that the corrupted data occurs in sparse and isolated time intervals. In addition to an initial guess of the solution and corresponding residual, which can be assigned using standard least squares, the TW approach requires a pre-determiniation of three additional parameters: a tolerance value to set the stopping criterion, threshold value $\lambda$ used in each iteration to set those parameters whose absolute values are below $\lambda$ to be zero, and another parameter $\mu$ to control the extent to which data points that do not (approximately) satisfy the prescribed model are to be considered as “corrupted data” and removed. For the Lorenz system example, we used the same parameters as in Ref. [49] whereas for the KSE example, we fix $\mu=0.0125$ (the same used in Ref. [49] and select $\lambda$ similarly as for the implementation of SINDy.

Implementation Details of Entropic Regression (ER)

As described in the main text, and as shown in details in Algorithm (1), a key quantity to compute in ER is the conditional mutual information $I(X;Y|Z)$ among three (possibly multivariate) random variables $X$, $Y$ and $Z$ via samples from these variables, denoted by $(x_{t},y_{t},z_{t})_{t=1,\dots,T}$. Since the distribution of the variables and their dependences are generally unknown, we adopt a nonparametric estimator for $I(X;Y|Z)$ which is based on statistics of $k$ nearest neighbors [30]. We fix $k=2$ in all of the reported numerical experiments; we have found that the results change quite minimally when $k$ is varied from this fixed value, suggesting relative robustness of the method.

Another important issue in practice is the determination of threshold under which the conditional mutual information $I(X;Y|Z)$ should be regarded zero. In theory $I(X;Y|Z)$ is always nonnegative and equals zero if and only if $X$ and $Y$ are statistically independent given $Z$, but such absolute criterion needs to be softened in practice because the estimated value of $I(X;Y|Z)$ is generally nonzero even when $X$ and $Y$ are indeed independent given $Z$. A common way to determine whether $I(X;Y|Z)=0$ or $I(X;Y|Z)>0$ is to compare the estimated value of $I(X;Y|Z)$ against some threshold. See Sec. (SI.3) for details of robust estimation of the threshold in the context of SID.

Supplementary Material

See supplementary material for more details in information theory measures, and additional numerical results for the double-well potential, Lorenz system, and a coupled network of the logistic map.

Acknowledgements

This work was funded in part by the Simons Foundation Grant No. 318812, the Army Research Office Grant No. W911NF-16-1-0081, the Office of Naval Research Grant No. N00014-15-1-2093, and also DARPA.

SI.1 Governing Dynamics, Over-sparsity, and Sensitivity for Expansion Order

In the main text we briefly discussed the effects of polynomial expansion order choses to construct the basis matrix $\Phi$, and in this section we provide results from numerical experiments to supplement these claims.

Recall from our main text Eq.(26), in Compressive Sensing (CS) framework we solves the constrained optimization problem:

[TABLE]

where the parameter $\epsilon\geq 0$ is used to relax the otherwise strict constraint $\Phi\bm{a}=\bm{f}$, to allow for the presence of noise in data.

Fig. SI.1 shows the reconstructed model by CS for the first equation of Lorenz system regarding $\dot{x}$. We observe sensitivity on expansion order and how CS yields over-sparse solution at the $7^{th}$ order expansion. This result was obtained using 300 measurements. To extend the investigation we repeat the same numerical experiment with doubled number of measurements, and Fig. SI.2 shows that CS still over-sparse the solution. This shows the relative sensitivity of CS with respect to the expansion order of basis functions.

In order to construct a second example that clearly shows the oversparse mechanism in CS, consider the three-dimensional linear system:

[TABLE]

It is easy to find that the solution for the above system is \bm{a}=\left(\begin{array}[]{ccc}0&2&0\end{array}\right)^{T}. Now, suppose that the third “measurement” is missing, and we have the under-determined system

[TABLE]

then we have infinitely many solutions lies on the line of intersection of the two planes:

[TABLE]

Figure SI.3 shows this simple example, where the solution for $\bm{a}$ lies on the intersection of the two planes shown, and we see the true solution, the LS solution and CS solution on the solution line. We see how the least square solution is far from the true solution with a high margin of error, but we also see that it only invest in $x$ and $y$ direction where the line of intersections of the two plans lies, then, LS ignore the $z$ direction and try to invest in all feasible directions to reach the best residual.

CS have different mechanism, since within all feasible solutions, it tends to select the one with minimum $\|\bm{a}\|_{1}$, even if there is another solution with the same number of sparse that have a residual $\|A\bm{a}-b\|_{2}=0$, and it is the case in our example where \|A\left(\begin{array}[]{ccc}0&2&0\end{array}\right)^{T}-b\|_{2}=0, while the CS solution has the residual $\|A\bm{a}_{CS}-b\|_{2}=2.5\times 10^{-5}$. In other words, for the system $A\bm{a}=b$, if there exist two solutions such that $\|\bm{a_{1}}\|_{1}<\|\bm{a_{2}}\|_{1}$ and $\|A\bm{a_{2}}-b\|_{2}<\|A\bm{a_{1}}-b\|_{2}\leq\epsilon$, where $\epsilon$ is the tolerance for CS optimization, then CS will select $a_{1}$ as a solution, even it has higher residual, and regardless of the structure of the sparse or the information flow between the basis functions and the observations. Numerically, assume the system in Eq. SI.13 to be:

[TABLE]

and consider a reasonable tolerance for CS solver to be $\epsilon=1e^{-9}$, then CS will always pick [1 0 0] as a solution even it have higher residual.

For many applications…, it is accepted to have such solution since it lies on the solution line and such residual difference will have negligible effect on the final result, But in discovering the governing equations of dynamical systems, such solution can often lead to a completely wrong structure of the system.

SI.2 Information Theory

In this section we review some basic concepts in information theory that underlie the development of entropic regression. These include classical concepts such as entropy and mutual information as well as recent developments of transfer entropy and causation entropy.

SI.2.1 Entropy

Entropy is firstly known as an extensive property of a thermodynamic system. The entropy of a thermodynamic system is a function of the microscopic states consistent with the macroscopic quantities that characterize the system. Assuming equal probability of the microscopic states, the entropy is given by:

[TABLE]

where $W$ is the number of microscopic states and $k_{B}$ is Boltzmann constant named after Ludwig Eduard Boltzmann where the Eq. SI.19 curved on his gravestone. Boltzmann saw entropy as a measure of statistical disorder in the system.

An analog to thermodynamic entropy is information entropy introduced by Claude Shannon in 1948 as “measures of information, choice, and uncertainty”. To describe Shannon’s entropy, consider a discrete random variable $X$ whose probability mass function is denoted by $p(x)=Prob(X=x)$. One can calculate its entropy as [12, 44],

[TABLE]

where $K$ is positive constant, and $H(X)$ is a measure of the uncertainty or unpredictability of $X$. Note that if we assume uniform probability distribution for the states of $X$, then we have $p(x)=\frac{1}{N}$, where $N$ is the number of states, and then Eq. SI.20 can be written as $H(X)=K\log(N)$ similar to Boltzmann’s entropy under the same assumption of equal probability of the states. The constant $K$, as Shannon sates, is merely amounts to a choice of a unit of measurement, and we consider $K=1$ for the rest of this document for simplicity. Fig. SI.4 shows the entropy function for a random event with different probability.

Shannon’s work provided extended and generalized view and understanding for the entropy, and one of the extended perspectives of Shannon’s entropy is dealing with the continuous random variables, and it takes the form:

[TABLE]

where $f_{X}(x)$ is the probability density function. The entropy shown in Eq. (SI.21) is referred to the differential entropy.

SI.2.2 Mutual Information

The entropy defined in Eq. (SI.20) naturally extends to the case of multiple random variables. For example, the joint entropy $H(X,Y)$, and conditional entropy $H(X|Y)$ of two random variables $X$ and $Y$ is given, respectively, by [12, 44],

[TABLE]

[TABLE]

where $p(x,y)$ is the joint probability distribution, and $H(X|Y)$ (read as entropy of $X$ given $Y$) is the measure of the uncertainty in $X$ if $Y$ is known. Some of the main properties of the entropy, joint entropy, and conditional entropy can be summarize as follows:

• •

The entropy of a discrete variable $X$ is positive ($H(X)\geq 0$), while the differential entropy does not satisfy this property.

• •

For two independent random variables $X$ and $Y$, $H(X,Y)=H(X)+H(Y)$.

• •

The chain rule: $H(X,Y)=H(X)+H(Y|X)$.

• •

One important property is that for a random variable $X$, the conditional entropy of $X$ given any other variable $Y$ will reduce the entropy of $X$, meaning that $H(X)\geq H(X|Y)$. The equality holds when $X$ and $Y$ are independent with $H(X,Y)=0$. This property tells that the information comes from $Y$ reduces the uncertainty about $X$, and when $Y=X$, that means we have given all the information about $X$, and then we become completely certain about $X$, and that gives $H(X|X)=0$.

The joint and conditional entropies can lead to a measures that detect the statistical dependence or independence between random variables. Such measure is called the mutual information between $X$ and $Y$, and it is given by [12, 44],

[TABLE]

where the mutual information $I(X;Y)$ (reads as mutual information between $X$ and $Y$) is a measure of the mutual dependence between the two variables. In terms of joint probability distribution, mutual information can be written as,

[TABLE]

and in its continuous form,

[TABLE]

where $f_{X,Y}(x,y)$ is the joint probability density function for the two continuous random variables $X$ and $Y$.

In case of independence of the two random variables, we have

[TABLE]

and then we have

[TABLE]

The same principle holds for the continuous variables in Eq. (SI.26), while $I(X;Y)$ satisfy the inequality $I(X;Y)\leqslant min[H(X),H(Y)]$ only in the discrete variables case.

SI.2.3 Transfer Entropy and Causation Entropy

For two stochastic processes $X_{t}$ and $Y_{t}$, the reduction of uncertainty about $X_{t+1}$ due to the information of the past $\tau_{Y}$ states of $Y$ , represented by

[TABLE]

in addition to the information of the past $\tau_{X}$ states of $X$ , represented by

[TABLE]

this reduction of uncertainty about $X_{t+1}$ is measured by “Transfer Entropy” which given by [12, 46],

[TABLE]

The traditional approach of inferring causality between two stochastic processes is to perform the Granger causality test [17]. The main limitation of this test is that it can only provide information about linear dependence between two processes, and therefore fails to capture intrinsic nonlinearities that are common in real-world systems. To overcome this difficulty, Schreiber developed the concept of transfer entropy between two processes [43]. Transfer entropy measures the uncertainty reduction in inferring the future state of a process by learning the (current and past) states of another process.

In our work [46, 48], we showed by several examples that causal relationship inferred by transfer entropy are often misleading when the underlying system contains indirect connections, a dominance of neighboring dynamics, or anticipatory couplings. For example, referring to the main text and the equation $\bm{f}=\Phi\bm{a}$, we see that the approaches that consider the transfer entropy in order to find the weak terms in $\Phi$ that has no influence on $\dot{X}$ to construct the sparse matrix $\bm{a}$, these approaches neglect a simple and clear idea that the terms of $\Phi$ has an indirect influence on $\bm{f}$ through the other terms of $\Phi$. To account for these effects, we developed a measure called Causation Entropy (CSE) [46, 48], and show that its appropriate application reveals true coupling structures of the underlying dynamics.

Consider a stochastic network of $N$ processes (nodes) denoted by:

[TABLE]

where $X_{t}^{(i)}\in\mathbb{R}^{d}$ is a random variable representing the state of process (or node) $i$ at time $t$, and $i\in\mathcal{V}=\{1,2,\dots,N\}$, and let $I,J$, and $K$ be a subsets of $\mathcal{V}$, then we can define the causation entropy as the following:

Definition 1 [48]: The causation entropy from the set of processes $J$ to the set of processes $I$ conditioning on the set of processes $K$ is defined as

[TABLE]

The Causation entropy is a natural generalization of transfer entropy from measuring pairwise causal relationships to network relationships of many variables. In particular, we can list the main properties for the causation entropy, noting that if $J=\{j\}$ and $I=\{i\}$, we simplify the notation as $C_{j\rightarrow i|K}$:

• •

If $j\in K$, then the causation entropy $C_{j\rightarrow i|K}=0$, as $j$ does not carry extra information (compared to that of $K$).

• •

If $K=\{i\}$, then the causation entropy recovers the transfer entropy $C_{j\rightarrow i|i}=T_{j\rightarrow i}$ which is given by $T_{j\rightarrow i}=H(X_{t+1}^{(i)}|X_{t}^{(i)})-H(X_{t+1}^{(i)}|X_{t}^{(i)},X_{t}^{(j)})..$

In [48], we introduced the principle of optimal Causation Entropy (oCSE) in a network of $N$ processes to find the minimum subset that maximizes the causation entropy. This minimal subset can be seen as the dominant subset of a network of $N$ processes, and they rule the underlying dynamic of the network. and in the same principle, we are looking for the dominant terms of the basis function $\Phi$ on the system dynamic $\bm{f}$. See (main text Fig.1) for visualization of the reformulation of Lorenz system in a network of processes.

SI.3 Entropic Regression

In our main text we discussed the Entropic Regression method and provided its Algorithm (main text Algorithm 1). In this section we discuss the tolerance estimation and its effect on the performance of ER.

In our previous work [48] we introduced a standard shuffle test, with a “confidence” parameter $\alpha\in[0,1]$ for tolerance estimation. The shuffle test requires randomly shuffling of one of the time series $n_{s}$ times, to build a test statistic. In particular, for the $i$-th random shuffle, a random permutation $\pi^{(i)}:[T]\rightarrow[T]$ is generated to shuffle one of the time series, say, $(y_{t})$, which produces a new time series $(\tilde{y}^{(i)}_{t})$ where $\tilde{y}^{(i)}_{t}=y_{\pi^{(i)}(t)}$; $(x_{t})$ is kept the same. Then, we estimate the mutual information $I(X;\tilde{Y}^{(i)})$ using the (partially) permuted time series $(x_{t},\tilde{y}^{(i)}_{t})$, for each $i=1,\dots,n_{s}$. For given $\alpha$, we then compute a threshold value $I_{\alpha}(X;Y)$ as the $\alpha$-percentile from the values of $I(X;\tilde{Y}^{(i)})$. If $I(X;Y)>I_{\alpha}(X;Y)$, we determine $X$ and $Y$ as dependent; otherwise independent.

We showed in [48], the robustness of shuffling test for optimal causation entropy calculations specially in complex dynamics, although it is computationally expensive. For more efficient computations complexity, we considered two different approach for tolerance estimation with the confidence parameter $\alpha$; the Dynamic, and the Static approaches.

In the Dynamic tolerance estimation, we start the forward step procedure (See main text Algorithm (1)) with initial tolerance $tol=0$, and we update the tolerance value at the end of the forward step procedure by the shuffle test shown in Algorithm (2) below. Given the confidence parameter $\alpha$, we update the tolerance to be the outcome of the shuffle test on the conditional mutual information of the forward step $I(\Phi R(\Phi,\bm{f},\{S_{f},j\});\bm{f}|\Phi R(\Phi,\bm{f},\{S_{f}\}))$, with $j$ indicates the index of maximum mutual information found by the current forward step.

The Static approach is more computationally efficient, where we apply the shuffle test to the mutual information $I(\bm{f};\bm{f})$ with the confidence parameter $\alpha$, and we assign the outcome of the shuffle test to the tolerance at the beginning of the forward step with no updates follows.

SI.4 Additional Numerical Results

To demonstrate the utility of ER for nonlinear system identification under noisy observations, we compare its performance against existing methods including the standard least squares (LS), orthogonal least squares (OLS), Lasso, and compressed sensing (CS). The details of the existing approaches are described in the Methods Section. The examples we consider represent different types of systems and scenarios, including both ODEs and PDEs, differential and difference equations, and network-coupled dynamics. In addition, we consider different noise models and especially the presence of outliers in order to evaluate the robustness of the respective methods.

For each example system, we sample the state of each variable at a uniform rate of $\Delta{t}$ to obtain a multivariate time series $\{z_{k}(t_{i})\}_{k=1,\dots,N;i=1,\dots,\ell}$; then we add noise to each data point and obtain the observed time series denoted by $\{\hat{z}_{k}(t_{i})\}$, where

[TABLE]

with $\eta_{ki}$ represents noise.

SI.4.1 Double Well Potential

In analogy to the example in our main text, we consider the equation

[TABLE]

and we sample 61 equally spaced measurements for $x\in[-1.2,1.2]$, and we construct $\Phi$ using the $10^{t}h$ order polynomial expansion with $K=11$ is the number of candidate functions. Then, we consider a single fixed value corrupted measurement to be $f(0.52)=0.5$.

In this example, we see that the true solution will have a residual $\delta$ equal to outliers deviation from its true position,

[TABLE]

Fig. SI.5 shows the result for LS. The LS with its BLUE property (Best Linear Unbiased Estimator), succeed to minimize the residual to have better fitting residual than the true solution, but it is clear that the residual value does not reflect reliable solution. Practically, when the true solution gives a fitting residual $\delta$, then any other solution deviates in its residual from $\delta$ will have a reduction in the solution accuracy, no matter the direction of deviation from $\delta$. In Fig. SI.6, we see the result of OLS. We see that the results with the best residual of OLS is almost identical to LS result. Here it worth to say a detailed review for the 1000 OLS solutions under different threshold showed us a small interval that gives solutions closer in structure to the true solution more than the minimum residual solution is shown, which is if treated with suitable trade-off strategy can give a better solution.

Fig. SI.7 shows the result for CS, where it failed to find any feasible solution for all values of $\epsilon<\delta$. Such outliers makes it hard to find a parameter vector $\bm{a}$ that can fit the data including the outliers point, and even with considering high resolution for $\epsilon$ span, so, CS as discussed before tends to select the solution with minimum $\|\bm{a}\|_{1}$ within the best feasible residuals. CS solution simulation for different outliers values is provided on our YouTube channel here. Fig. SI.8 shows the result for LASSO, and it shows the sparse solution with wrong structure of LASSO. We considered the bounds of $\lambda$ to be $\lambda\in\left[\|\Phi\Phi^{\dagger}\bm{f}-\bm{f}\|,\|\bm{f}\|\right]$, where $\lambda=\|\Phi\Phi^{\dagger}\bm{f}-\bm{f}\|$ is the penalty on the solution with all entries are non-sparse and $\lambda=\|\bm{f}\|$ is the penalty on the solution with all entries are sparse. For this example, different from others (see Methods section), and because of its small dimensions, we considered very large span (1000 values) of the tununing parameter value for OLS, LASSO and CS to investigate the best expected outcomes of the methods.

Fig. SI.9 shows the accurate structure found by ER. Even with a slight difference in the magnitude of the parameters, we see how ER recovers the true basis functions. The residual of ER was $0.865$, which is higher more than most other methods, but the ER focuses on the information flow between the basis and dynamic and not the residual of solution magnitudes.

SI.4.2 Lorenz system.

SI.4.3 Coupled Network of Logistic map

Our third example is a network of coupled logistic maps which is a generalization of coupled map lattices [27] and also cellular automata [28]. This scenario of high dimensional and complex systems that have become the thrust of recent analysis including in the synchronization literature [36, 2, 23]. In this example, we assume that not only the governing dynamics are unknown, but so is the structure of the network that moderates the coupling between individual chaotic elements; both of these must be (simultaneously) identified from observed dynamic data alone. In Fig. SI.16, we compare the results of several system identification methods, including the proposed ER approach. We now offer here a rough description of why this dramatic difference in performance, in the setting particular here of noisy data subject to outliers; a more detailed mathematical analysis will be the subject of our future work.

Consider that each of these other methods we reviewed involves minimizing a functional $J({\mathbf{a}})$ of the data $a$, and that when ${\mathbf{a}}$ is subject to noise, that the functionals are each continuous with respect to their argument. We assume that the underlying system is,

[TABLE]

describing the individual elements as Logistic maps, but, the coupled network of $N$ such oscillators is of the form,

[TABLE]

where $i,j=1,...,N$, $A$ is the adjacency matrix of the coupled network, $k$ is the global coupling strength, and $f(x_{i})$ is the image of the point $x_{i}$ under the logistic map given in Eq.SI.35.

To present a specific example, let $N=50$, we construct the adjacency matrix $A$ to have simple coupling such that:

[TABLE]

where $D$ is the degree matrix of $A$, and the coupling adjacency matrix $A$ constructed randomly such that the above inequality holds. Fig. (SI.17) show the graph of the coupled network. Then if we consider only the second order expansion (where the basis matrix $\Phi$ is the second order expansion of the 50 time-series of all nodes) we will have 1326 terms in our expansion matrix. We focus on this example on solving the underdetermined system by considering 2000 measurements as maximum available measurements. So, exclude OLS which only solve overdetermined systems and cannot be investigated at a number of measurements less than 1326, and we exclude LASSO and CS for their high computation complexity. Fig. SI.16 shows the error in recovered parameters for this example. For simplicity and the computation complexity, we performed the experiment to find the parameters for one single dimension, and results are averaged over 50 runs.

This example shows the robustness of ER in recovering the coupling structure in complex coupled networks. The computations complexity in such problem can be highly reduced by considering basic and trivial assumptions. For example, we can consider each node $N_{i}$ as a default influence source for itself, and then instead of starting the forward step in ER from the empty set, we may initialize the index set with the terms that purely includes $N_{i}$.

Fig.(SI.18) shows the sparse representation of the Logistic map discussed with number of nodes $N=20$.