Regime switching is ubiquitous in many complex dynamical systems with multiscale features, chaotic behavior, and extreme events. In this paper, a causation entropy boosting (CEBoosting) strategy is developed to facilitate the detection of regime switching and the discovery of the dynamics associated with the new regime via online model identification. The causation entropy, which can be efficiently calculated, provides a logic value of each candidate function in a pre-determined library. The reversal of one or a few such causation entropy indicators associated with the model calibrated for the current regime implies the detection of regime switching. Despite the short length of each batch formed by the sequential data, the accumulated value of causation entropy corresponding to a sequence of data batches leads to a robust indicator. With the detected rectification of the model structure, the subsequent parameter estimation becomes a quadratic optimization problem, which is solved using closed analytic formulas. Using the Lorenz 96 model, it is shown that the causation entropy indicator can be efficiently calculated, and the method applies to moderately large dimensional systems. The CEBoosting algorithm is also adaptive to the situation with partial observations. It is shown via a stochastic parameterized model that the CEBoosting strategy can be combined with data assimilation to identify regime switching triggered by the unobserved latent processes. In addition, the CEBoosting method is applied to a nonlinear paradigm model for topographic mean flow interaction, demonstrating the online detection of regime switching in the presence of strong intermittency and extreme events.
Online nonlinear system identification with sequential data has recently become important in many applications, e.g., extreme weather events, climate change, and autonomous systems. In this work, we developed a causation entropy boosting (CEBoosting) framework for online nonlinear system identification. For each sequential data batch, this framework calculates the causation entropy that evaluates the contribution of each function in a large set of candidate functions to the system dynamics. The causation entropies based on multiple data batches are then aggregated to identify a few candidate functions that have significant impacts on the system dynamics. With the identified sparse set of functions, the framework further fits a model of the system dynamics. We tested the proposed framework for complex systems with features including chaotic behavior, high dimensionality, intermittency and extreme events, and partial observations. The results show that the CEBoosting method can capture the regime switching and then fit models of system dynamics for various types of complex dynamical based on a limited amount of sequential data.
I. INTRODUCTION
Regime switching is ubiquitous in many complex dynamical systems in geoscience, engineering, neural science, and material science.1–7 The switching is usually associated with sudden changes in internal states or appears when specific external forcing is exerted. Dynamical systems often display distinct behavior with regime switching. One example is the atmospheric jets, which meander in different directions when the atmosphere alternates between blocked and unblocked regimes.8,9 Similarly, an excitable medium is susceptible to finite perturbations, which triggers regime switching from a quiescent state to one with various wave patterns.10–12 Regime switching can also induce an increased occurrence of extreme events, leading to, for example, extreme weather and climate patterns,13,14 bursting neurons,15 or extreme ductile damages.16 Detecting regime switching and the corresponding underlying dynamics, which rely on appropriate model identification methods, have significant social and scientific impacts. Challenges in detecting regime switching are associated with the intrinsic properties of many complex dynamical systems, including high dimensionality, partial or incomplete observations, and the intermittent occurrence of rare and extreme events.17–21
Efficient model identification has received significant attention. Both physical knowledge and observational data facilitate learning the underlying model dynamics. The model structures are often established utilizing physical intuitions for traditional knowledge-based model identification. The primary process then becomes the estimation of model parameters. Linear models are natural candidates for simple problems22,23 and can potentially be skillful for short-term forecasts. Other families of models with pre-determined structures, such as the physics-constrained nonlinear regression models24,25 and conditional Gaussian nonlinear systems,26,27 are alternative nonlinear models aiming to capture specific underlying dynamical features. On the other hand, recent progress has been made in data-driven model identification. Data-driven reduced-order models have been widely used in scientific and engineering applications.28–31 Sparse model discovery methods also appear as advanced model identification tools that allow automatic learning of the model structure and parameters from data and lead to nonlinear models with parsimonious structures via sparse regression.32–39 Recently, an information-flow-based approach for system identification has gained popularity, which utilizes causation entropy to identify the sparse model structure of an unknown dynamical system.40 An iterative approach called entropic regression which is robust with noise and outliers has further been developed and applied to the area of cognitive neuroscience to improve the accuracy of recovery of the structure of structural and functional connectivity between anatomical regions.41–43 With a limited amount of indirect data, derivative-free optimization methods44–46 have been explored as model identification tools. In addition, non-parametric and machine learning models have been built to characterize complex dynamical systems.47–52
Among various model discovery approaches, online model identification is a particularly useful method in practice, which sequentially determines model structure and estimates model parameters when new observation arrives.53–58 It is the primary model identification strategy in many geophysical and engineering problems, where the limited amount of historical data is insufficient to robustly discover the underlying dynamics. It should be noted that online model identification can be further combined with data assimilation to handle noisy observations or recover the unobserved state variables in the situation with partial observations.59–61 However, unlike the online parameter estimation that can be efficiently addressed by standard filtering methods,62 the lack of knowledge about the proper model structure poses a unique challenge in online model identification. The sequentially arriving data play a vital role in progressively rectifying the model and reducing the uncertainty in the identified system. Although existing system identification methods can identify the proper model structure via promoting sparsity in the offline setting by fitting a model to abundant data, promoting sparsity relies on the model fitting may eliminate some important model structures in the context of sequential learning. To address such a challenge, we incorporate causation entropy to achieve robust online model identification. As regime switching often occurs and completes within a short transient period, developing suitable online identification methods for discovering regime switching exploiting transition data is essential with practical importance.
In this paper, a causation entropy boosting (CEBoosting) strategy is developed. It is incorporated into an online model identification method to detect regime switching and discover the nonlinear dynamics associated with the new regime. Different from many existing sparse model identification algorithms, such as those relying on LASSO (least absolute shrinkage and selection operator) regression46,63 or thresholding,32,35 the method developed here separates the estimation of model parameters from the recurrent identification of nonlinear model structure. Such a separation allows using closed analytic formulas for the entire online learning algorithm and, therefore, the overall computational cost is significantly reduced. In this new strategy, causation entropy64,65 is utilized to provide a logic value (i.e., true or false) of each candidate function in a pre-determined library throughout the online learning process. By examining the causation entropy on the newly arrived data, the reversal of one or a few such causation entropy indicators associated with the model calibrated for the current regime implies the detection of regime switching. In other words, the causation entropy indicator, which can be efficiently calculated, is employed to decide if the existing terms in the current model need to be rectified and if the system demands additional terms as a response to regime switching. Note that the sequential data in online learning are collected within a short time window to form a batch of time series, which is utilized to compute the causation entropy. As each batch contains a short amount of data, it may embody only part of the dynamical properties. Nevertheless, as time evolves, the accumulated value of causation entropy corresponding to a sequence of batches leads to a robust indicator of the model structure in response to regime switching. The concept of accumulating causation entropy calculated from sequential data relates to the statistical method of bagging.66 With the detected rectification of the model structure, the subsequent parameter estimation becomes a quadratic optimization problem, which is solved using closed analytic formulas. For multiple times of regime switching, a summation of residual models is calibrated, which relates to the statistical method of boosting.67–70
The proposed new strategy has several unique features. First, causation entropy takes into account the interdependence between all the candidate functions in the pre-determined library and, therefore, it can eliminate the superficial causal relationship. Model identification exploiting the causation entropy has been shown to reach a higher selection accuracy than LASSO regression or elastic net.65 The causation-based learning approach also indicates robust results in the presence of indirect coupling between features and stochastic noise,71 which are crucial features of complex systems. Second, causation entropy is only utilized to indicate the terms that need to be added or removed from the existing model. In other words, although computing the exact value of the causation entropy is challenging, closed analytic formulas are available for efficiently approximating this causation entropy indicator, which allows an effective detection of the model structure. Third, the parameter estimation only needs to be carried out after the model structure is entirely determined. Therefore, the overall computational cost is reduced compared with applying LASSO regression, which requires detecting the model structure and estimating model parameters simultaneously for each batch of data. Lastly, compared to the pioneering works of using causation entropy for system identification (e.g., Refs. 40, 64, and 72), this work focuses on dynamical systems with regime switching, and the proposed CEBoosting method utilizes online data to (i) identify the regime switching and (ii) estimate the sparse model structure to calibrate the unknown model parameters for the new regime. Therefore, the key contribution of this work is to demonstrate the use of causation entropy for a robust online sparse system identification with the sequential and limited amount of data, by leveraging statistical concepts of bagging and boosting. It is worth highlighting that the causation entropy indicator is easy to calculate and applicable to moderately large dimensional systems. The method developed here is also adaptive to the situation with partial observations, where utilizing data assimilation to recover the unobserved state variables can be incorporated into the learning process for identifying regime switching resulting from the latent processes. In addition, the method is not limited to the Gaussian data. It can be applied to dynamical regimes with strong intermittency and extreme events. Applications to nonlinear dynamical systems with moderately large dimensions, partial observations, and extreme events are all studied in the paper.
The remainder of the paper is organized as follows. The development of the online identification method utilizing the CEboosting strategy is presented in Sec. II. Section III includes four test cases. In addition to a standard chaotic model as a proof-of-concept, the other three test cases emphasize the method applying to systems with moderately large dimensions, partial observations, and extreme events, respectively. The paper is concluded in Sec. IV.
II. METHODOLOGY
As regime switching often results from a sudden change of a small number of the model parameters or specific components of the model structure, the residual model typically has a sparse structure that can be calibrated using relatively short data, which is precisely the case of the online identification problems. Therefore, instead of learning the entire model associated with the new regime, the focus is to estimate the residual part . Once the residual part is identified, it is then added to the existing model that provides the new system as a response to the regime switching.
The limited amount of data is assumed to arrive sequentially in the form of batches. The th batch represents data [or a subset of the vector in the partial observation case] for . It should be noted that, as is typically unknown in practice, the identification algorithm usually does not start from (namely, the left point of the first interval ). Yet, for the simplicity of presentation, is chosen to be for the numerical examples in this work. This will not affect the identification algorithm as applying the algorithm to those batches prior to will not indicate regime switching. But this setup facilitates counting for the length of the data that is needed to detect the regime switching once it occurs at .
A. Causation entropy
B. Causation entropy boosting (CEBoosting) algorithm
With this stable aggregated causation entropy matrix , we then impose sparsity into by setting when and extract a set of remaining coefficients .
After the sparse parameter matrix for residual dynamics model (2.9) is obtained, the current model is updated by adding the information from the residual dynamics . The CEBoosting algorithm repeats the above procedure when a new batch of data arrives. A schematic illustration of Lorenz 63 system with regime switching is displayed in Fig. 1, and more details of the CEBoosting algorithm can be found in Algorithm 1 presented in Appendix A. It should be noted that the notations adopted in this section assume the regime switching time for the simplicity of the illustration. In practice, regime switching can happen at with , for which Algorithm 1 presents the detailed procedures of detecting regime switching, aggregating causation entropy matrix, identifying a sparse model structure, and then fitting the model parameters.
Schematic of CEBoosting algorithm (based on the Lorenz 63 model). Panel (a): data are generated from the Lorenz 63 system with regime switching. The parameter in regime 1 is , while is changed to 38 in regime 2. The index of incoming batch data starts with . Panel (b): is the current model parameter matrix defined in (2.3) with being the basis functions and being the dynamic of state . With , the residual dynamics and the causation entropy between and can be calculated for each batch. is the aggregated causation entropy from batch to . is the binary matrix of defined in (2.10) indicating the sparse structure of the residual model. is defined in (2.11) indicating number that becomes stable, i.e., the pattern is consistent with the previous aggregated causation entropy matrix. Panel (c): with and data batches from the batch with new regime detected to the one with a stable , the residual model is calibrated by least square estimation.
Schematic of CEBoosting algorithm (based on the Lorenz 63 model). Panel (a): data are generated from the Lorenz 63 system with regime switching. The parameter in regime 1 is , while is changed to 38 in regime 2. The index of incoming batch data starts with . Panel (b): is the current model parameter matrix defined in (2.3) with being the basis functions and being the dynamic of state . With , the residual dynamics and the causation entropy between and can be calculated for each batch. is the aggregated causation entropy from batch to . is the binary matrix of defined in (2.10) indicating the sparse structure of the residual model. is defined in (2.11) indicating number that becomes stable, i.e., the pattern is consistent with the previous aggregated causation entropy matrix. Panel (c): with and data batches from the batch with new regime detected to the one with a stable , the residual model is calibrated by least square estimation.
C. Comparison with SINDy
In recent years, SINDy has been a popular framework for sparsely identifying nonlinear dynamics from data.32 SINDy exploits an iterative thresholding regularization to determine the model structure and estimate the model parameters simultaneously. The iterative thresholding regularization guarantees the parsimonious model structure. Using abundant training data (e.g., long time-series), the original SINDy method was designed for offline model identification. To better work with a limited amount of data, a recent extension of SINDy34 leveraged the statistical approach of bagging and achieved a more robust learning performance.
Compared to the SINDy methods, the CEBoosting algorithm is mainly designed with a low computational cost for robust online learning with limited sequential data. The key difference from the SINDy methods is that imposing sparsity is decoupled from the parameter estimation in the CEBoosting algorithm. Specifically, causation entropy is utilized only to determine a parsimonious model structure without dealing with the model calibration. The concept of bagging is further introduced with sequential batches of data to ensure robust estimation of causation entropy before imposing a parsimonious model structure. It results in utilizing the minimum data to determine the model structure. With a robust estimation of the parsimonious model structure, the CEBoosting algorithm exploits a simple least square estimate to solve the parameter values, a quadratic optimization problem with a closed analytic solution. This avoids repeatedly estimating the parameters in the LASSO-type or thresholding-based regression approaches when determining the model structure by examining each batch data. In online learning, the CEBoosting algorithm, by design, can also avoid accidentally removing essential basis functions due to the incorrectly calibrated model based on a limited amount of data.
III. NUMERICAL EXPERIMENTS FOR SYSTEMS WITH MODERATELY LARGE DIMENSIONS, PARTIAL OBSERVATIONS, AND EXTREME EVENTS
In this section, the performance of the CEBoosting method is demonstrated utilizing four different chaotic or turbulent systems, which are models that mimic crucial features in many science or engineering disciplines. The presentation starts with the three-dimensional Lorenz 63 system, which is a classical chaotic system. This test is utilized as a proof of concept to demonstrate the detailed steps of the method in the online identification of the non-linear dynamics with regime switching. To illustrate the efficiency of the algorithm in capturing the regime switching behavior in a relatively high-dimensional case, the 40-dimensional Lorenz 96 system is utilized as a second test, where the localization technique is incorporated to mitigate the curse of dimensionality. As strongly non-Gaussian statistics, intermittency and extreme events appear in many climate, atmosphere, and ocean science problems, a multi-mode layered topographic model that captures these crucial turbulent features is adopted as the next test model. The last test case aims to deal with a more realistic scenario where only a subset of the state variables is observed. Data assimilation is, therefore, essential in such a partial observational case. A stochastic parameterized extended Kalman filter (SPEKF) model is used to estimate and simulate the hidden system state variables. Below, assuming the starting model and its parameters are available for all the numerical examples. The goal is to learn the regime switching and the corresponding residual model that adjusts the original model to a new regime based on online sequential data.
For all four systems, the numerical simulation time step size is chosen as . The hyper-parameter in (2.11) is chosen as 4 for the numerical examples of this work.
A. The Lorenz 63 system: A classical chaotic system
Figure 2(a) compares the trajectories in the phase space between the original (regime 1) and new (regime 2) regimes. It can be seen that the two systems are on different manifolds in the phase space. Figure 2(c) shows the time series of the three state variables. With the regime switching at , it can be seen in Fig. 2(c) that the patterns of time series change accordingly, especially for the variable that demonstrates a shift of its mean value. The autocorrelation function (ACF) of each state variable for the original system is presented in Fig. 2(b), which shows a rapid decay of correlation within one time unit, except for the variable that has some oscillations in its ACF. Similar behavior of the ACF is observed in Fig. 2(d) after the regime switching. Finally, Fig. 2(e) shows the ensemble mean of the state variable as a function of time from an independent simulation with 5000 ensemble members, which reveals that the transition time of the regime switching is about 20 time units. This also indicates that the transition time depends on the property of the transient feature and is very different from the decorrelation time. Nevertheless, the decorrelation time of the original system provides a natural way to determine the batch size. Thus, the batch size is chosen as one time unit here.
Lorenz 63 system with regime switching. (a) Trajectories of the original system and the new one in phase space. (c) Time series of system state variables. (b) and (d) Autocorrelation function (ACF) of the original system and the new regime. (e) Ensemble mean of variable before and after regime switching at .
Lorenz 63 system with regime switching. (a) Trajectories of the original system and the new one in phase space. (c) Time series of system state variables. (b) and (d) Autocorrelation function (ACF) of the original system and the new regime. (e) Ensemble mean of variable before and after regime switching at .
The CEBoosting algorithm is employed to detect the regime switching and identify the sparse structure of the residual model. The candidate basis functions include all the linear and quadratic nonlinear functions: . With the identified sparse structure, the coefficients of the residual model are then determined via the least square estimation. To detect the regime switching and identify the sparse structure of the residual model, the aggregated causation entropy matrix is gradually updated until its structure gets stable, i.e., the causation entropy values of some candidate basis functions keep being significantly greater than others. The stable causation entropy matrix structure is obtained after 6 time units and is presented in Table I, where only one entry has a significant value. It is worthwhile to highlight that the total time units (6 units) to discover the regime switching and determine the model parameters in the new regime is shorter than the transition time (20 units). In other words, utilizing the information from the transient period is sufficient for the CEBoosting algorithm to determine the model response to the regime switching in this test case.
The L63 model—The causation entropy matrix (CEM) after 6 time units. The pattern of CEM becomes stable and does not change with incorporating more batch data. The entry with a significant value of the causation entropy is highlighted using the bold font. According to the pattern of this CEM, a residual model will be built.
10−4 . | x . | y . | z . | xy . | xz . | yz . | x2 . | y2 . | z2 . |
---|---|---|---|---|---|---|---|---|---|
11.0437 | 10.6292 | 5.1156 | 4.6996 | 6.7041 | 9.2123 | 8.9308 | 3.0866 | 4.0888 | |
63.1474 | 3.7597 | 5.0283 | 2.5584 | 2.6370 | 5.5017 | 4.3689 | 2.4472 | 4.4099 | |
9.0692 | 6.6605 | 9.5399 | 11.3740 | 8.5529 | 10.6165 | 10.8299 | 10.8725 | 10.5916 |
10−4 . | x . | y . | z . | xy . | xz . | yz . | x2 . | y2 . | z2 . |
---|---|---|---|---|---|---|---|---|---|
11.0437 | 10.6292 | 5.1156 | 4.6996 | 6.7041 | 9.2123 | 8.9308 | 3.0866 | 4.0888 | |
63.1474 | 3.7597 | 5.0283 | 2.5584 | 2.6370 | 5.5017 | 4.3689 | 2.4472 | 4.4099 | |
9.0692 | 6.6605 | 9.5399 | 11.3740 | 8.5529 | 10.6165 | 10.8299 | 10.8725 | 10.5916 |
Finally, based on the sparse structure identified in Table I, a residual model can be calibrated via the least square estimation based on the residuals of the original model. The residual model is then added to the original model as a correction term. The coefficients of the corrected model for the new regime are listed in Table II. It can be seen that the corrected model successfully captures the crucial parameter value that induces the new regime.
The L63 model—Updated model parameters for the new regime.
. | x . | y . | z . | xy . | xz . | yz . | x2 . | y2 . | z2 . |
---|---|---|---|---|---|---|---|---|---|
−9.9749 | 9.9671 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
38.0357 | −1.0129 | 0 | 0 | 0 | −1.0016 | 0 | 0 | 0 | |
0 | 0 | −2.6670 | 1.0017 | 0 | 0 | 0 | 0 | 0 |
. | x . | y . | z . | xy . | xz . | yz . | x2 . | y2 . | z2 . |
---|---|---|---|---|---|---|---|---|---|
−9.9749 | 9.9671 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
38.0357 | −1.0129 | 0 | 0 | 0 | −1.0016 | 0 | 0 | 0 | |
0 | 0 | −2.6670 | 1.0017 | 0 | 0 | 0 | 0 | 0 |
B. The Lorenz 96 system: Strategy for applying CEBoosting to relatively high-dimensional systems
The original system (regime 1) has a forcing coefficient , which makes the system have a strongly chaotic behavior. After the regime switching at , the new system (regime 2) takes the forcing value . Meanwhile, the coefficient of , representing the damping effect, changes from to . This leads to a fully turbulent regime.111 Some weak coherent structures can still be observed in the strongly chaotic regime, but they disappear in the fully turbulent one. See Fig. 3 for the two regimes. Figure 4 shows the detailed statistical properties of the state variable . Note that the model has homogeneous dynamics, so the statistical properties for different state variables are the same. Figure 4(a) illustrates how the time series of change with the regime switching. From the strongly chaotic to the fully turbulent regime, the variance of all the state variables becomes more extensive, and the decorrelation time becomes shorter. The batch size of 1 time unit is chosen here, which is again of the same order as the decorrelation time. The total transition time is about 1–2 units, according to Fig. 4(b).
The Lorenz 96 system with regime switching. Top: 40-dimensional system state evolving at the time interval with a regime switching at . Regime 1 is chaotic with , and regime 2 is turbulent with and the linear terms . Bottom: zoom-in view of the regime switching at the time interval .
The Lorenz 96 system with regime switching. Top: 40-dimensional system state evolving at the time interval with a regime switching at . Regime 1 is chaotic with , and regime 2 is turbulent with and the linear terms . Bottom: zoom-in view of the regime switching at the time interval .
Statistical properties of the state variable of the Lorenz 96 system with regime switching. Panel (a): time series of in regimes 1 and 2. Panels (b) and (c): the PDFs of in both regimes. Panels (d) and (e): the ACFs of in both regimes. Panels (f) and (g): the ensemble mean and the ensemble variance of before and after the regime switching at .
Statistical properties of the state variable of the Lorenz 96 system with regime switching. Panel (a): time series of in regimes 1 and 2. Panels (b) and (c): the PDFs of in both regimes. Panels (d) and (e): the ACFs of in both regimes. Panels (f) and (g): the ensemble mean and the ensemble variance of before and after the regime switching at .
A set of candidate basis functions is built by exploiting polynomials up to the second order for each state . Note that the residual dynamics of each state is assumed to have contributions from the basis functions that consist of only its adjacent states. The CEBoosting algorithm identifies the sparse structure of the residual model using about two batches of data (i.e., the causation entropy matrix pattern does not change after two time units). Table III shows the CEM based on utilizing two batches of data, which is about the same length as the estimated transition time of Lorenz 96 system [see Figs. 4(f) and 4(g)]. It is seen that the causation entropy entries associated with the actual residual terms are at least one order more significant than others. Therefore, the CEBoosting algorithm successfully identifies the correct sparse structure of the residual model based on data within the time interval of the transient period.
The L96 model—The CEM after 2 time units. The entries with significant values of the causation entropy are highlighted using the bold font. The pattern of CEM becomes stable and does not change with incorporating more data batches. According to the pattern of this CEM, a residual model will be built.
10−4 . | x1 . | . | x2x40 . | x39x40 . | x2 . | . | x3x1 . | x40x1 . | ⋅⋅ ⋅ . | x40 . | . | x1x39 . | x38x39 . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
436.0597 | 4.2015 | 2.2964 | 1.4823 | 0.0976 | 3.4716 | 1.3792 | 0.9251 | ⋅⋅ ⋅ | 0.1869 | 3.0008 | 0.7236 | 8.0802 | |
8.4920 | 9.8888 | 6.0428 | 8.0793 | 19.8015 | 8.3095 | 10.1492 | 9.6844 | ⋅⋅ ⋅ | 0 | 0 | 0 | 0 | |
⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ |
2.7772 | 3.1363 | 4.5714 | 8.4405 | 9.6327 | 5.3498 | 4.8701 | 7.6239 | ⋅⋅ ⋅ | 200.3710 | 2.5200 | 7.6225 | 2.7748 |
10−4 . | x1 . | . | x2x40 . | x39x40 . | x2 . | . | x3x1 . | x40x1 . | ⋅⋅ ⋅ . | x40 . | . | x1x39 . | x38x39 . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
436.0597 | 4.2015 | 2.2964 | 1.4823 | 0.0976 | 3.4716 | 1.3792 | 0.9251 | ⋅⋅ ⋅ | 0.1869 | 3.0008 | 0.7236 | 8.0802 | |
8.4920 | 9.8888 | 6.0428 | 8.0793 | 19.8015 | 8.3095 | 10.1492 | 9.6844 | ⋅⋅ ⋅ | 0 | 0 | 0 | 0 | |
⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ |
2.7772 | 3.1363 | 4.5714 | 8.4405 | 9.6327 | 5.3498 | 4.8701 | 7.6239 | ⋅⋅ ⋅ | 200.3710 | 2.5200 | 7.6225 | 2.7748 |
After obtaining the sparse structure according to the CEM in Table III, the residual dynamics are calibrated via a linear combination of candidate basis functions, and the coefficients for the linear combination are obtained by the least square estimation. The L96 model in the new regime is then updated by adding this residual model to the original model. The coefficients of the updated model are summarized in Table IV. The identified model shows a good agreement with the true system of the new regime with and the linear term .
The L96 model—Coefficients of the updated model for the new regime. The entries with bold font indicate the parameters that contribute to the regime switching.
. | 1 . | x1 . | . | x2x40 . | x39x40 . | x2 . | . | x3x1 . | x40x1 . | ⋅⋅ ⋅ . | x40 . | . | x1x39 . | x38x39 . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
16.0248 | −1.4938 | 0 | 1.0000 | −1.0016 | 0 | 0 | 0 | 0 | ⋅⋅ ⋅ | 0 | 0 | 0 | 0 | |
16.0490 | 0 | 0 | 0 | 0 | −1.5077 | 0 | 1.0000 | −1.0000 | ⋅⋅ ⋅ | 0 | 0 | 0 | 0 | |
⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ |
16.0027 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ⋅⋅ ⋅ | −1.5017 | 0 | 1.0000 | −1.0000 |
. | 1 . | x1 . | . | x2x40 . | x39x40 . | x2 . | . | x3x1 . | x40x1 . | ⋅⋅ ⋅ . | x40 . | . | x1x39 . | x38x39 . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
16.0248 | −1.4938 | 0 | 1.0000 | −1.0016 | 0 | 0 | 0 | 0 | ⋅⋅ ⋅ | 0 | 0 | 0 | 0 | |
16.0490 | 0 | 0 | 0 | 0 | −1.5077 | 0 | 1.0000 | −1.0000 | ⋅⋅ ⋅ | 0 | 0 | 0 | 0 | |
⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ | ⋅⋅ ⋅ |
16.0027 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ⋅⋅ ⋅ | −1.5017 | 0 | 1.0000 | −1.0000 |
C. The topographic model: Dynamical system with intermittency and extreme events
Figure 5 displays the regime switching and statistical properties of all five variables in the topographic model. The pattern and amplitude of each time series demonstrate a noticeable change after the regime switching happens at . Note that the long time series of both regimes are utilized here, ensuring that the associated statistics are computed with a sufficiently large number of sample points. The PDFs behave like Laplace distributions with heavy tails, indicating many extreme events in the model simulation. The ACF of in the original model (regime 1) decays very slowly, taking about 30 time units until the ACF approaches zero. Because of this, a batch size of 30 time units is utilized in the CEBoosting algorithm. Due to the strengthening of the topographic effect after the regime switching, the time series of and in the new model (regime 2) display multiscale features, where the ACFs have a quick decay at the beginning but then relax to zero slowly after 50 time units.
Topographic model with regime switching. Middle column: time series of each state in the time interval [0, 5000] with regime switching at . First and fifth columns: probability density function (PDF) of each state in both regimes. Second and fourth columns: autocorrelation function (ACF) of each state in both regimes.
Topographic model with regime switching. Middle column: time series of each state in the time interval [0, 5000] with regime switching at . First and fifth columns: probability density function (PDF) of each state in both regimes. Second and fourth columns: autocorrelation function (ACF) of each state in both regimes.
Similar to the previous test models, the library contains the polynomials up to the second order. In other words, the following 20 basis functions are employed to build the library: .
As the noise levels of state variables are significantly different from each other, the magnitudes of causation entropy for different residual dynamics are not the same either. Therefore, the selection of the candidate functions is based on the causation entropy values for the dynamics of each state variable separately. To determine the entire residual model, 60 time units is used. See Table V.
The topographic model—The CEM after 60 time units. The pattern of CEM does not change with incorporating more batch data. The entries with significant values of the causation entropy are highlighted using the bold font.
10−3 . | v1 . | v2 . | v3 . | v4 . | u . | v1u . | v2u . | v3u . | v4u . | ⋅⋅ ⋅ . |
---|---|---|---|---|---|---|---|---|---|---|
0.0115 | 0.0227 | 0.0165 | 0.0614 | 120.9879 | 0.0216 | 0.0365 | 0.0407 | 0.0249 | ⋅⋅ ⋅ | |
0.0124 | 0.0327 | 0.0066 | 0.0150 | 0.0313 | 0.0002 | 0.0090 | 0.0039 | 0.0284 | ⋅⋅ ⋅ | |
0.0023 | 0.0033 | 0.0142 | 0.0099 | 51.3949 | 0.0537 | 0.0291 | 0.0572 | 0.0082 | ⋅⋅ ⋅ | |
0.0359 | 0.0619 | 0.0150 | 0.0134 | 0.0253 | 0.0180 | 0.0127 | 0.0099 | 0.0210 | ⋅⋅ ⋅ | |
0.1080 | 0.0022 | 0.0471 | 0.0146 | 0.0065 | 0.0068 | 0.0097 | 0.0049 | 0.0063 | ⋅⋅ ⋅ |
10−3 . | v1 . | v2 . | v3 . | v4 . | u . | v1u . | v2u . | v3u . | v4u . | ⋅⋅ ⋅ . |
---|---|---|---|---|---|---|---|---|---|---|
0.0115 | 0.0227 | 0.0165 | 0.0614 | 120.9879 | 0.0216 | 0.0365 | 0.0407 | 0.0249 | ⋅⋅ ⋅ | |
0.0124 | 0.0327 | 0.0066 | 0.0150 | 0.0313 | 0.0002 | 0.0090 | 0.0039 | 0.0284 | ⋅⋅ ⋅ | |
0.0023 | 0.0033 | 0.0142 | 0.0099 | 51.3949 | 0.0537 | 0.0291 | 0.0572 | 0.0082 | ⋅⋅ ⋅ | |
0.0359 | 0.0619 | 0.0150 | 0.0134 | 0.0253 | 0.0180 | 0.0127 | 0.0099 | 0.0210 | ⋅⋅ ⋅ | |
0.1080 | 0.0022 | 0.0471 | 0.0146 | 0.0065 | 0.0068 | 0.0097 | 0.0049 | 0.0063 | ⋅⋅ ⋅ |
Based on the sparse model structure identified in Table V, the residual model is calibrated via the least square estimation. The coefficients of the calibrated model for the new regime are summarized in Table VI. It can be seen that the corrected model successfully updates the changed parameters and in the new regime. The resulting model can reproduce the strong non-Gaussian features with intermittency and extreme events.
The topographic model—Updated model for the new regime.
. | v1 . | v2 . | v3 . | v4 . | u . | v1u . | v2u . | v3u . | v4u . | ⋅⋅ ⋅ . |
---|---|---|---|---|---|---|---|---|---|---|
−0.0496 | −1.0008 | 0 | 0 | −4.2463 | 0 | 0.9998 | 0 | 0 | ⋅⋅ ⋅ | |
0.9988 | −0.0501 | 0 | 0 | 0 | −1.0009 | 0 | 0 | 0 | ⋅⋅ ⋅ | |
0 | 0 | −0.0554 | −0.4968 | −2.1200 | 0 | 0 | 0 | 2.0015 | ⋅⋅ ⋅ | |
0 | 0 | 0.4998 | −0.0523 | 0 | 0 | 0 | −1.9996 | 0 | ⋅⋅ ⋅ | |
2.1123 | 0 | 4.5713 | 0 | −0.0599 | 0 | 0 | 0 | 0 | ⋅⋅ ⋅ |
. | v1 . | v2 . | v3 . | v4 . | u . | v1u . | v2u . | v3u . | v4u . | ⋅⋅ ⋅ . |
---|---|---|---|---|---|---|---|---|---|---|
−0.0496 | −1.0008 | 0 | 0 | −4.2463 | 0 | 0.9998 | 0 | 0 | ⋅⋅ ⋅ | |
0.9988 | −0.0501 | 0 | 0 | 0 | −1.0009 | 0 | 0 | 0 | ⋅⋅ ⋅ | |
0 | 0 | −0.0554 | −0.4968 | −2.1200 | 0 | 0 | 0 | 2.0015 | ⋅⋅ ⋅ | |
0 | 0 | 0.4998 | −0.0523 | 0 | 0 | 0 | −1.9996 | 0 | ⋅⋅ ⋅ | |
2.1123 | 0 | 4.5713 | 0 | −0.0599 | 0 | 0 | 0 | 0 | ⋅⋅ ⋅ |
It is worthwhile to remark that, as the noise levels in the dynamics of are lower than that in , a shorter time series (and shorter batch length) with in total of only 15 time units can be utilized to reach a stable CEM for the components. See Table VII. Figure 6 shows the time evolution of the ensemble mean and the ensemble variance of the topographic model, including the time instant of regime switching. The results here can be utilized to infer the transition time, which is about 20 time units. Therefore, in this topographic model with a large noise in the dynamics, the stable pattern of all variables can be identified by the CEBoosting algorithm within the transition period. Yet, 60 time units of data are needed for identifying the model structure associated with the variable , mainly due to the larger noises of .
Ensemble mean and variance of each state variable in the topographic model. Regime switching happens at .
Ensemble mean and variance of each state variable in the topographic model. Regime switching happens at .
The topographic model—The CEM after 15 time units. The pattern of all v variables is stable. The entries with significant values of the causation entropy are highlighted using the bold font.
10−3 . | v1 . | v2 . | v3 . | v4 . | u . | v1u . | v2u . | v3u . | v4u . | ⋅⋅ ⋅ . |
---|---|---|---|---|---|---|---|---|---|---|
0.0182 | 0.0233 | 0.0322 | 0.0679 | 22.665 | 0.0244 | 0022 | 0.1726 | 0.0285 | ⋅⋅ ⋅ | |
0.0768 | 0.0071 | 0.0002 | 0.0004 | 0.0879 | 0.0582 | 0.0007 | 0.0171 | 0.1059 | ⋅⋅ ⋅ | |
0.0187 | 0.0006 | 0.1123 | 0.0411 | 10.3880 | 0.0133 | 0.0024 | 0.0843 | 0.0166 | ⋅⋅ ⋅ | |
0.0087 | 0.1908 | 0.0171 | 0.0023 | 0.000 | 0.0001 | 0.0520 | 0.0015 | 0.0016 | ⋅⋅ ⋅ | |
0.0008 | 0.0007 | 0.0078 | 0.0194 | 0.0167 | 0.0250 | 0.0054 | 0.0223 | 0.0107 | ⋅⋅ ⋅ |
10−3 . | v1 . | v2 . | v3 . | v4 . | u . | v1u . | v2u . | v3u . | v4u . | ⋅⋅ ⋅ . |
---|---|---|---|---|---|---|---|---|---|---|
0.0182 | 0.0233 | 0.0322 | 0.0679 | 22.665 | 0.0244 | 0022 | 0.1726 | 0.0285 | ⋅⋅ ⋅ | |
0.0768 | 0.0071 | 0.0002 | 0.0004 | 0.0879 | 0.0582 | 0.0007 | 0.0171 | 0.1059 | ⋅⋅ ⋅ | |
0.0187 | 0.0006 | 0.1123 | 0.0411 | 10.3880 | 0.0133 | 0.0024 | 0.0843 | 0.0166 | ⋅⋅ ⋅ | |
0.0087 | 0.1908 | 0.0171 | 0.0023 | 0.000 | 0.0001 | 0.0520 | 0.0015 | 0.0016 | ⋅⋅ ⋅ | |
0.0008 | 0.0007 | 0.0078 | 0.0194 | 0.0167 | 0.0250 | 0.0054 | 0.0223 | 0.0107 | ⋅⋅ ⋅ |
D. A stochastic parameterized extended Kalman filter (SPEKF) model: Incorporating data assimilation into the CEBoosting algorithm
The SPEKF model (3.4) has been widely used as an approximate model to describe a spectral mode of a complex turbulent system, especially in the context of data assimilation and ensemble prediction.111,114,116 Physically, the variable represents one of the resolved modes (i.e., observable) in the turbulent signal, while the three hidden variables , and are surrogates for the nonlinear interaction between and other unobserved modes in the original governing equation after applying the spectral decomposition. The idea of the SKEPF model is that the small or unresolved scale variables are stochastically parameterized by inexpensive linear and Gaussian processes, representing stochastic damping , stochastic phase , and stochastic forcing . Despite the model error in using such Gaussian approximations for the original unresolved nonlinear dynamics, these Gaussian processes succeed in providing accurate statistical feedback from the unresolved scales to the resolved ones. Thus, the intermittency and non-Gaussian features observed in the resolved variables can be accurately recovered. The statistics in the SPEKF model can also be solved with exact and analytic formulas, which allow an accurate and efficient estimation of the model states. The SPEKF type of model has been used for filtering multiscale turbulent dynamical systems,111 stochastic superresolution,117 and filtering Navier–Stokes equations with model error.118 It has been shown that the SPEKF model has much higher skill than classical Kalman filters using the so-called mean stochastic model (MSM) to capture the irregularity and intermittency in nature.
In the following, the system will experience regime switching twice. In the starting regime (regime 1), all three hidden variables have significant contributions to the dynamics of the observed process . In regime 2, will be set to zero; therefore, the stochastic phase does not influence the dynamics of the resolved variable . Finally, in regime 3, the stochastic damping will be removed from the dynamics of , but the contribution from will be added back. Figure 7 shows the observed trajectories and the associated statistics of the real part of the observed variable. Due to the stochastic damping, the time series in regimes 1 and 2 are intermittent with multiple extreme events when the overall damping becomes positive. As a result, the PDFs are non-Gaussian fat-tailed. In contrast, very few extreme events are found in regime 3, and the associated PDF is Gaussian. Similarly, the ACF in regime 2 shows a clear oscillatory pattern when it decays. This is due to a dominant frequency of the time series coming from the constant phase . Such a regular oscillation in the ACF becomes less significant when the stochasticity is added to the phase term via .
The SPEKF model with three different regimes (twice the regime switching) as time evolves. In regime 1, all three hidden variables , , and have significant contributions to the dynamics of . In regime 2, is set to zero. In regime 3, the stochastic damping is removed from the dynamics of u and the contribution from is added back. First row: the time series of the real part of the observed variable . Second row: the PDF of in each regime. Third row: the ACF of in each regime.
The SPEKF model with three different regimes (twice the regime switching) as time evolves. In regime 1, all three hidden variables , , and have significant contributions to the dynamics of . In regime 2, is set to zero. In regime 3, the stochastic damping is removed from the dynamics of u and the contribution from is added back. First row: the time series of the real part of the observed variable . Second row: the PDF of in each regime. Third row: the ACF of in each regime.
The goal here is to detect the regime switching and reveal the dynamics in each regime. It is worth highlighting that the dimensions of the system in the three regimes are different. In the absence of or , the dimension reduces from to from regime 1 to regimes 2 and 3. However, such a change is unknown in practice and relies on the learning algorithm to detect it. In particular, there is no observed time series of , , and . Therefore, recovering these variables becomes an essential step in the online learning algorithm. To this end, the following procedure is adopted that incorporates data assimilation into the CEBoosting algorithm. Assume the model in regime 1 is known. Each time when the new batch of time series of is obtained, such a model is utilized to sample a trajectory of the three unobserved variables , , and conditioned on the observed trajectory of . The sampled trajectory can be thought of as the analog of one ensemble member of the ensemble Kalman smoother solution of the system conditioned on the observed signal of ,119 although the sampled trajectory using the SPEKF model can be written down using closed analytic formulas.27 See Ref. 39 for the implementation details of using such closed analytic formulas for data assimilation. This augments the unobserved components of each batch of data. Then, the CEBooting algorithm is utilized to compute the causal relationships in light of the observed time series of and the sampled time series of , , and . If the causation entropy from any function involving to the dynamics of is zero, then the equation is eliminated from the final model structure. A similar logic applies to and . Note that as and are complex-valued variables, the actual computation regards the real and imaginary parts as two processes in computing the causation entropy.
The main focus here is on identifying the dynamics of . In particular, we aim to investigate if all three unobserved variables contribute to the observed process of . To this end, the linear Gaussian models of , , and are assumed to be fixed. A function library is designed with { , , , , , , , , , , } to learn the structure of .
Figure 8 shows the real part of with regime switching and sampled trajectory of , , across the three regimes. When the unobserved variables , , or contribute to the dynamics of , the sampled processes match the truth quite well. On the other hand, when and disappear in regimes 2 and 3, respectively, the sampling result provides random trajectories. The causation entropy from the terms involving these trajectories has no contribution to the dynamics of . The first row of Table VIII, showing the causation entropies of regime 1, is based on a time series with 200 time units. The second and the third rows show the causation entropies for the detected regimes 2 and 3, respectively, using one batch of data with a length of 20 units. It is seen that the term has a nearly zero causation entropy to in regime 2, where is sampled. Similarly, the sampled leads to a nearly zero causation entropy from to in regime 3.
Trajectory sampling for the SPEKF model with twice the regime switching. The black curve represents the conditional sampling of each hidden process given the observed variable in all three regimes.
Trajectory sampling for the SPEKF model with twice the regime switching. The black curve represents the conditional sampling of each hidden process given the observed variable in all three regimes.
The SPEKF model—The CEM with Data Assimilation. First row: causation entropy with existing regime 1 data (200 time units). Second row: first batch (20 time units) causation entropy of regime 2. Third row: first batch (20 time units) causation entropy of regime 3.
. | uγ . | uω . | b . | u2 . | γω . | γ2 . | ω2 . | b2 . | γb . | ωb . | ub . |
---|---|---|---|---|---|---|---|---|---|---|---|
0.2527 | 0.5272 | 0.2289 | 0.0077 | 0.0006 | 0.0028 | 0.0003 | 0.0011 | 0.0090 | 0.0084 | 0.0137 | |
0.4525 | 0.0208 | 0.1163 | 0.0491 | 0.0070 | 0.0080 | 0.0019 | 0.0015 | 0.0113 | 0.0035 | 0.0066 | |
0.0844 | 0.7032 | 0.3486 | 0.0048 | 0.0027 | 0.0019 | 0.0017 | 0.0001 | 0.0068 | 0.0106 | 0.0063 |
. | uγ . | uω . | b . | u2 . | γω . | γ2 . | ω2 . | b2 . | γb . | ωb . | ub . |
---|---|---|---|---|---|---|---|---|---|---|---|
0.2527 | 0.5272 | 0.2289 | 0.0077 | 0.0006 | 0.0028 | 0.0003 | 0.0011 | 0.0090 | 0.0084 | 0.0137 | |
0.4525 | 0.0208 | 0.1163 | 0.0491 | 0.0070 | 0.0080 | 0.0019 | 0.0015 | 0.0113 | 0.0035 | 0.0066 | |
0.0844 | 0.7032 | 0.3486 | 0.0048 | 0.0027 | 0.0019 | 0.0017 | 0.0001 | 0.0068 | 0.0106 | 0.0063 |
Figure 9 shows the conditional mean and uncertainty (two standard deviations) of and from the smoother solution. It is seen that when and contribute to the dynamics of in regime 1, the conditional mean follows the truth quite well. However, the conditional mean of and behaves like a random trajectory in regimes 2 and 3, respectively, with relatively large uncertainty. Thus, the sampled trajectories are formed from randomness and uncertainty, with no causal inference on the observed variable .
The smoother mean and the smoother uncertainty (two standard deviations) of and from the SPEKF model with twice the regime switching. The black curve shows the posterior mean of each hidden process conditioned on the observed . The shaded gray area corresponds to the 95% confidence interval (two standard deviations) of the corresponding mean estimate.
The smoother mean and the smoother uncertainty (two standard deviations) of and from the SPEKF model with twice the regime switching. The black curve shows the posterior mean of each hidden process conditioned on the observed . The shaded gray area corresponds to the 95% confidence interval (two standard deviations) of the corresponding mean estimate.
IV. CONCLUSION
Online nonlinear system identification with sequential data has recently become important in many applications, e.g., extreme weather events, climate change, and autonomous systems. In this work, we developed a causation entropy boosting (CEBoosting) framework for online nonlinear system identification. The CEBoosting algorithm aims to (i) discover a sparse residual model structure based on the aggregated causation entropy calculated from sequential data and (ii) calibrate the residual model with the identified sparse structure via least square estimation. If the true system experiences multiple regime switching, the proposed framework gradually identifies a summation of residual models, which has a close analogy to the statistical technique of boosting. We tested the proposed framework for complex systems with features including chaotic behavior, high dimensionality, intermittency and extreme events, and partial observations. The results show that the CEBoosting method can capture the regime switching and then calibrate residual models for various types of complex dynamical based on a limited amount of sequential data. It is worth noting that constraints can be naturally added to the learning algorithm. One important constraint is the so-called physics constraint,24 which requires the total energy in the quadratic nonlinear terms to be conserved. It guarantees the long-term stability of the identified system. Such a constraint has yet to be explicitly incorporated into the current framework, although the resulting parameters in various non-Gaussian test cases shown in this work already roughly satisfy this constraint. Adding constraints can be easily achieved by imposing simple relationships between model parameters in the parameter estimation step, which still allows using closed analytic formulas for finding the parameters. See, for example, Ref. 39 for details. Other future work includes the uncertainty quantification of the aggregated causation entropy and the further study of other causality metrics.
ACKNOWLEDGMENTS
N.C. is partially funded by Office of Naval Research N00014-19-1-2421 and Army Research Office W911NF-23-1-0118. J.W. and C.C. are supported by the University of Wisconsin-Madison, Office of the Vice Chancellor for Research and Graduate Education with funding from the Wisconsin Alumni Research Foundation.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Chuanqi Chen: Conceptualization (supporting); Investigation (lead); Methodology (supporting); Software (lead); Writing – original draft (supporting); Writing – review & editing (supporting). Nan Chen: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Jin-Long Wu: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. The codes and examples that support the findings of this study are available in the link: https://github.com/ChuanqiChenCC/CEBoosting.120
APPENDIX A: CEBoosting ALGORITHM
Algorithm 1 presents the detailed procedures of the CEBoosting method, including (i) detecting regime switching, (ii) aggregating causation entropy matrix (CEM) until a consistent pattern of CEM is obtained, (iii) identifying a sparse model structure according to the aggregated CEM, and (iv) fitting the model parameters. The consistent pattern of CEM is determined by (2.11), i.e., the aggregated CEM does not change for data batches. In this work, we choose for all the numerical examples.
Causation Entropy Boosting.
1: Current Model |
2: for do |
3: |
4: |
5: Residual dynamics |
6: |
7: if then Same regime; Output current model |
8: |
9: else New regime detected |
10 , , |
11: while do Aggregate CEM until it is consistent for iterations |
12: |
13: |
14: Threshold aggregated CEM to get with values |
15: if then |
16: |
17: else |
18: |
19: end if |
20: end while |
21: Smallest that satisfies the criterion in (2.11) |
22: Select a sparse model structure and estimate parameters |
23: |
24: |
25: |
26: end if |
27: return |
28: end for |
1: Current Model |
2: for do |
3: |
4: |
5: Residual dynamics |
6: |
7: if then Same regime; Output current model |
8: |
9: else New regime detected |
10 , , |
11: while do Aggregate CEM until it is consistent for iterations |
12: |
13: |
14: Threshold aggregated CEM to get with values |
15: if then |
16: |
17: else |
18: |
19: end if |
20: end while |
21: Smallest that satisfies the criterion in (2.11) |
22: Select a sparse model structure and estimate parameters |
23: |
24: |
25: |
26: end if |
27: return |
28: end for |
APPENDIX B: DATA ASSIMILATION
The data assimilation technique27 used in Sec. III D (SPEKF model) of numerical experiments has been summarized here.
(The notation is the complex conjugate transpose.)
[The notation means solving (B4) backward from end point of interval [0,T] with and .]
The SPEKF model (3.4) is following the structure of conditional Gaussian non-linear system (B2): corresponds to and corresponds to with , , , , , , , .
With the observed variable , three hidden variables will be sampled according to the sampling formula (B5) and presented in Fig. 8. The smoother estimation in Fig. 9 is from the smoother formula (B4). The sampled trajectory is the one that satisfies the posterior distribution at each time instant and considers the temporal correlation at different time instants. The formula is used in generating the trajectories of the unobserved variables for computing the causation entropy in Sec. III D.