How entropic regression beats the outliers problem in nonlinear system identification

Following the standard routine in nonlinear SID,¹⁰ 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.^3,7 A common choice is the standard polynomial basis

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 is referred to as the structure of the model.

Consider time series data ${z(t)=[z1(t),…,zm(t)]⊤}t=t0,…,tℓ$ and corresponding ${F(z(t)}t=t0,…,tℓ$ generated from a nonlinear, high-dimensional dynamical system (1), possibly subject to observational noise. From $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 $Fi(z(tk))=zi(tk+1)−zi(tk)τk+O((tk+1−tk))$ or central difference which has improved accuracy: $Fi(z(tk))=zi(tk+1)−zi(tk−1)tk+1−tk−1+O((tk+1−tk−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 $F$⁠.

Under the basis representation (2), the identification of $F$ becomes equivalent to estimating all the parameters ${aik}$⁠. In practice, the empirically observed state variable is subject to noise: $z^(t)=z(t)+η(t)$ with $η(t)$ representing the (multivariate) noise and $F^i$ denoting the approximated value of $Fi$⁠. For noisy observations $z^(t)$⁠, the difference between $F^i(z^(t))$ and $Fi(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 $ξ(t)$ term and express the forward model as

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

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,

where $f(i)=[F^i(z^(t1)),…,F^i(z^(tℓ))]⊤∈Rℓ×1$ represents the $i$th component of the estimated vector field from the observational data, $Φ=[ϕ(1),…,ϕ(K)]∈Rℓ×K$ [with $ϕ(k)=[ϕk(z^(t1)),…,ϕk(z^(tℓ))]∈Rℓ×1$] represent sampled data for the basis functions, $ξ(i)=[ξi(t1),…,ξi(tℓ)]⊤∈Rℓ×1$ represents noise, and $a(i)=[ai1,…,aiK]⊤∈RK×1$ is the vector of parameters, which is to be determined. Note that the form of Eq. (6) is the same for each $i$⁠, and solving each $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:

where $f∈Rℓ×1$ and $Φ∈Rℓ×K$ are given, with the goal to estimate $a∈RK×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 (OLS), lasso, and 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.

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 Akaike information criterion (AIC),¹¹ ours combines model construction with selection. In addition, even though it is not uncommon for entropy measures to be adopted in system identification,^8,12 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^13,14 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 (7). For an index set $S⊂N∪{0}$⁠, the estimated parameters can be thought of as a mapping from the joint space of $Φ$⁠, $f$⁠, and $S$ to a vector denoted as $a^=R(Φ,f,S)$⁠. For instance, under a least squares criterion, the mapping is given by $R(Φ,f,S)S=ΦS†f$ (⁠$ΦS$ denotes the columns of matrix $Φ$ indexed by $S$⁠) and $R(Φ,f,S)i=0$ for all $i∉S$⁠. Using the estimated parameters, the recovered signal can be computed as $ΦR(Φ,f,S)$⁠. In the ER algorithm, we start by selecting a basis function $ϕk1$ that maximizes its mutual information with $f$⁠, compute the corresponding parameter $ak1$ using the least squares method, and obtain the corresponding regression model output $z1$ according to

Here, $I(x;y)$ denotes mutual information between $x$ and $y$⁠, which is a model-free measure of the statistical dependence between two distributions (that is, $x$ and $y$ are independent if and only if their mutual information equals zero; however, in practice, due to finite samples and statistical estimation, we wish to distinguish that Mutual Information (MI) is statistically insignificantly indistinguishable from zero and noting that it is never negative as well).¹⁵ Next, in each iteration of the forward stage, we perform the following computations and updates for $i=2,3,…$⁠:

The process terminates when $maxkI(ΦR(Φ,f,k);f|zi−1)≈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={k1,…,km}$ together with the corresponding parameters $ak1,…,akm$ (⁠$aj=0$ for $j∉S$⁠) and model $f≈ak1ϕk1+⋯+akiϕki$⁠. Finally, we turn to the backward stage, where the terms that had previously been included are re-examined for their information-theoretic relevance and those that are redundant will be removed. In particular, we sequentially check for each $j=ki∈S$ to determine if the basis term $ϕj$ is redundant by computing

and updating $S→S/{j}$ (that is, remove $j$ from the set $S$⁠) if $I(ΦR(Φ,f,S);f|z¯j)≈0$⁠. The result of the backward ER is the reduced set of indices $S={ℓ1,…,ℓn}$ with $n≤m$⁠, together with the corresponding parameters $aℓ1,…,aℓn$ (⁠$aj=0$ for $j∉S$⁠) computed as $a=R(Φ,f,S)$⁠, and accordingly the recovered model $f≈ϕa=ϕSaS=aℓ1ϕℓ1+⋯+aℓnϕℓ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 section (also see Supplementary Materials).

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 Ordinary Differential Equations (ODEs) and Partial Differential Equations (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 $Δt$ to obtain a multivariate time series ${z(ti)}k=1,…,N;i=1,…,ℓ$⁠, where $z=[z1,…,zd]⊤∈Rd$⁠; then, we add noise to each state variable and obtain the noisy empirical time series denoted by ${z^(ti)}$⁠, where

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

First, we compare several nonlinear SID methods in reconstructing the Lorenz system when the state observational noise is drawn independently from a Gaussian distribution, $η∼N(0,ε2)$⁠. As we discussed before, this translates into effective noise that is not necessarily Gaussian or even independent. Figure 2 shows the error in the estimated parameters, where $error=‖atrue−aestimated‖2$⁠. As shown in Fig. 2, even with observational noise as low as $ε=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 Refs. 9 and 2 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 $104$⁠. See Fig. 3 for the results of the TW method with a large span of tuning parameters.

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 $η$⁠, where $η∼N(0,ε12)$ with probability $1−p$ and $η∼N(0,ε12+ε22)$ with probability $p$⁠. Here, we use $ε1=10−5$⁠, $ε2=0.2$⁠, and $p=0.2$⁠. The results of SID are shown in Figs. 4 and 5. 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 retains its capacity of accurately identifying the underlying system. As an example, in each of the side panels of Fig. 4, 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 undersampled regime, the ER method successfully recovers the system dynamic. As an ample amount of data becomes available, we note that the TW method starts to produce excellent reconstruction, which is consistent with recent findings reported in Ref. 9.

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 structures identified by the different methods? In Fig. 6, we compare the structure of the identified model using different methods across a range of parameter values for $ρ$⁠. 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 appear 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 discuss in the supplementary material, 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.

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 ODEs in the corresponding Banach space.

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

where $q=2πL$⁠.

Notice that we have written this Fourier series of basis elements $eikqx$ in terms of time varying combinations of basis elements. For simplicity, consider $L=2π$⁠, then $q=1$ for the following analysis. This is typical²⁶ with the representation of a PDE as infinitely many ODEs in the Banach space, where orbits of these coefficients, therefore, become time varying patterns by Eq. (13). Substituting Eq. (13) into Eq. (12), we produce the infinitely many evolution equations for the Fourier coefficients,

In general, the coefficients $bk$ are complex functions of time $t$⁠. However, by symmetry, we can reduce to a subspace by considering the special symmetry case that $bk$ is purely imaginary, $bk=iak$ and $ak∈R$⁠. Then,

where $k=1,…,Nm$⁠. However, the assumption that there is a slow manifold (slow modes as an inertial manifold^26–29) suggests the practical matter that a finite truncation of the series Eq. (13), and correspondingly the reduction to finitely many ODEs will suffice. Therefore, we choose a sufficiently large number of modes $Nm$⁠. Then, we solve the resulting $Nm$-dimensional ODE (15) to produce the estimated solution of $u(x,t)$ by (13), and use such data for the purpose of SID, so as to estimate the structure and parameters of the ODE model (15).

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

Figure 8 shows the sparse structure of the recovered solution by different methods. Here, we mention that the true nonzero parameters of Kuramoto-Sivashinsky equations (KSE) using $Nm=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×107$⁠. 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 $λ$ that is slightly above $λ∗$⁠, where $λ∗≈0.1731$ is the smallest magnitude of the true nonzero parameter of the full least squares solution. A larger value of $λ$ 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 the ER solution has a few false negatives, it was completely able to recover the overall dynamic of the system as shown in Fig. 9, while all other solutions diverges and failed to recover $u(x,t)$⁠.

Figure 10 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. 10 shows the results for CS and ER. The results for each solver and details are provided in the supplementary material.

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 nonrobustness to 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).¹⁵ The outcome of this principle is that sampled data can be partitioned into “typical” samples and “atypical” samples, with the rare atypical samples ending 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 interpretation 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.

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 $F$⁠, expressed as

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

Here, $ξ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 $zi(t)$ can translate into correlated non-Gaussian effective noise for $ξi(t)$⁠.

A system identification problem can be transformed into parameters estimation problem (or inverse problem) in the form of

where $f(i)=[F^i(z^(1)),…,F^i(z^(T))]⊤∈RT×1$ represents the estimated function $Fi$ (⁠$i$-th component of the vector field $F$⁠), $Φ=[ϕ(1),…,ϕ(K)]∈RT×K$ (with $ϕ(k)=[ϕk(z^(1)),…,ϕk(z^(T))]∈RT×1$⁠) represent sampled data for the basis functions, $ξ(i)=[ξi(1),…,ξi(T)]⊤∈RT×1$ represents effective noise, and $a(i)=[ai1,…,aiK]⊤∈RK×1$ is the vector of parameters, which is to be determined. Since the form of the Eq. (19) is the same for each $i$⁠, we omit the index when discussing the general methodology, and consider the following linear inverse problem:

where $f∈RT×1$ and $Φ∈RT×K$ are given, with the goal is to estimate $a∈RK×1$ when the effective noise is not necessarily from independent multivariate Gaussian distribution.

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

The solution can be explicitly computed, giving

where $Φ†$ denotes the pseudoinverse of the matrix $Φ$⁠.³⁰ 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 $a$ that has the minimal $2$-norm and solves the equation $Φa=f$⁠. The LS method has several advantages: it is analytically traceable and easy to solve computationally using standard linear algebra routines (e.g., Singular Value Decomposition [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 $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.

In orthogonal least squares (OLS),^4,31,32 the idea is to iteratively select the columns of $Φ$ that minimize the (2-norm) model error, which corresponds to iterative assigning nonzero values to the components of $a$⁠. In particular, the first step is to select basis $ϕk1$ and compute the corresponding parameter $ak1$ and residual $r1$ according to

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

As for stopping criteria, there are several choices including AIC and Bayesian information criterion (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 the 5-fold cross validation.

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

where the parameter $λ≥0$ controls the extent to which sparsity is desired: as $λ→∞$⁠, the second term dominates and the only solution is a vector of all zeros, whereas at the other extreme, $λ=0$ and the problem becomes identical to a least squares problem, which generally yields a full (nonsparse) solution. Values of $λ$ 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 $λ$ needs to be tuned in order to specify a particular solution. A common way to select $λ$ is via cross validation.³³ In our numerical experiments, we choose $λ$ span according to Ref. 33, with the 5-fold cross validation and 10 values of $λ$ span. We adopt the CVX solver,³⁴ and from all the solutions found for each $λ$⁠, we select the solution with minimum residual.

Originally developed in the signal processing literature,^35–37 the idea of compressed sensing (CS) has been adopted in several recent works in the nonlinear system identification.^6,7 Under the CS framework, one solves the following constrained optimization problem:

where the parameter $ε≥0$ is used to relax the otherwise strict constraint $Φa=f$⁠, to allow for the presence of noise in data. In our numerical experiments, we choose 10 log-spaced values for $ε∈[10−6,100]$⁠, and the 5-fold cross validation. We adopt the CVX solver,³⁴ and from all the solutions found for each $ε$⁠, we select the solution with minimum residual.

In their recent contribution, Brunton, Proctor, and Kutz introduced SINDy (Sparse Identification of Nonlinear Dynamics) as a way to perform nonlinear system identification.² Their main idea is, after formulating the inverse problem (20), 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 $λ$⁠, the method starts from a least squares solution and abandons all basis functions whose corresponding parameter in the solution has absolute value smaller than $λ$⁠; 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 $λ$ 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 $λ=0.02$ as used in a similar example as in Ref. 2.

In their recent paper,⁹ 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 occur 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 predetermination of three additional parameters: a tolerance value to set the stopping criterion, threshold value $λ$ used in each iteration to set those parameters whose absolute values are below $λ$ to be zero, and another parameter $μ$ 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. 9, whereas for the KSE example, we fix $μ=0.0125$ (the same used in Ref. 9 and select $λ$ similarly as for the implementation of SINDy.

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 $(xt,yt,zt)t=1,…,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.³⁸ 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 non-negative and equals zero if and only if $X$ and $Y$ are statistically independent given $Z$⁠, but such an 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. (??) for details of robust estimation of the threshold in the context of SID.

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

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 (ONR) (Grant No. N00014-15-1-2093), and the DARPA.