# Model selection of chaotic systems from data with hidden variables using
sparse data assimilation

H. Ribera, S. Shirman, A. V. Nguyen, N. M. Mangan

(August 28, 2024)

###### Abstract

Many natural systems exhibit chaotic behaviour such as the weather, hydrology,
neuroscience and population dynamics. Although many chaotic systems can be
described by relatively simple dynamical equations, characterizing these
systems can be challenging, due to sensitivity to initial conditions and
difficulties in differentiating chaotic behavior from noise. Ideally, one
wishes to find a parsimonious set of equations that describe a dynamical
system. However, model selection is more challenging when only a subset of the
variables are experimentally accessible. Manifold learning methods using time-
delay embeddings can successfully reconstruct the underlying structure of the
system from data with hidden variables, but not the equations. Recent work in
sparse-optimization based model selection has enabled model discovery given a
library of possible terms, but regression-based methods require measurements
of all state variables. We present a method combining variational annealing –
a technique previously used for parameter estimation in chaotic systems with
hidden variables – with sparse optimization methods to perform model
identification for chaotic systems with unmeasured variables. We applied the
method to experimental data from an electrical circuit with Lorenz-system like
behavior to successfully recover the circuit equations with two measured and
one hidden variable. We discuss the robustness of our method to varying noise
and manifold sampling using ground-truth time-series simulated from the
classic Lorenz system.

Significance statement

Chaos represents a challenge for studying the dynamic behavior of many
physical and biological systems. Since the 80s we have known that time-series
measurements from one variable of a chaotic system contain information about
the underlying structure of the full multi-dimensional system. However,
recovery of the full system from data with hidden variables has remained
elusive. This work develops a novel data-assimilation technique to identify
governing equations of chaotic systems from data with hidden variables. This
method identifies fairly simple, low-dimensional, and deterministic models
from seemingly incomplete data. Discovery of such equations can enable rich
mathematical study and physical insight for problems across nearly every
discipline including climate science, hydrology, neuroscience, ecology,
medicine and engineering.

##  1 Introduction

Hypothesis generation through data-driven model identification has the
potential to revolutionise science. Uncovering the interactions, structure,
and mechanisms that determine the behaviour of chaotic systems in particular
could improve scientific understanding in almost every discipline with
dynamical systems [30] including climate science [60], hydrology [59],
population dynamics [32], and neuroscience [53]. Many chaotic systems can be
informatively described by relatively simple dynamical equations. However,
characterization and control of these systems can be challenging [11], due to
sensitivity to initial conditions and difficulties in differentiating chaotic
behavior from noise [64]. Characterization through statistical, geometric, or
model-based means becomes more challenging when only a subset of the variables
are experimentally accessible. Our goal is to identify a parsimonious set of
equations to describe a chaotic system from measurements with hidden
variables.

Much data-analysis for chaotic systems has focused on learning the attracting
manifold structure from time-series. In the early 80s, Takens’s theorem [65]
describes the conditions under which one can use the time-delay embedding from
a single variable to construct a manifold that preserves the topological
properties of the full system. Takens’s result formalized the idea that the
information of the manifold structure, and therefore chaotic dynamics, could
be recovered from the time-history of a single state variable. Manifold
reconstruction methods [50, 28, 41] based on partial information provide
insight into the system structure, dimensionality, and statistics of chaotic
systems. By constructing manifolds from time-delays, Sugihara _et al._
developed methods discriminating chaos from noise [64] and detecting causality
between measured variables [63]. Methods including reservoir computing [67,
27], other deep learning frameworks [72], data assimilation combined with
neural networks [15], support vector machine [48], and nearest neighbours [5]
can accurately predict the dynamics of chaotic systems using a data-trained
model with no specific physical knowledge of the system. For a review of
predictive methods see [6]. Assuming a reasonable model structure is known,
data-assimilation methods [8, 9] including variational annealing [70] can
estimate model parameters for chaotic systems from incomplete, indirect, and
noisy measurements. Although these methods are designed to assimilate
information from data-streams with hidden variables and learn about chaotic
systems, they are not designed for the purpose of hypothesizing parsimonious
models or identifying model structure.

Data-driven discovery of parsimonious dynamical systems models to describe
chaotic systems is by no means new. Early on, least-squares fitting of
combinatorial sets of polynomial basis functions to time-series data followed
by information-theory based selection produced models that reproduced manifold
structure and statistics of the system [21]. Symbolic regression demonstrated
successful recovery of the widely accepted equations for the chaotic double-
pendulum system [58]. More recently sparse regression [68, 17, 34], motivated
model selection techniques such as SINDy [16], which recover the ground-truth
equations for chaotic systems from a relatively large library of functions,
without needing a computationally intensive combinatorial search. Other
sparsity-promoting frameworks have improved upon robustness for chaotic
systems equation recovery through integral formulations [57, 51, 47], data
assimilation methods [12], Bayesian frameworks [13], and entropic regression
[4]. However, all these methods require measurements of all state-variables
that significantly impact the desired dynamic. Notably, Champion _et al._
recently used an autoencoder framework for automatic discovery system
coordinates and equations, but required input time-series of a higher
dimension than intrinsic dimension of the system [18].

Model selection with hidden variables require different methodology. By
‘hidden variables’ we mean that the number of measured variables is smaller
than the intrinsic dimension of the system. Measured variables are not
considered hidden if they are corrupted by noise or indirectly sampled through
a measurement function. A few methods address the problem of model selection
with hidden variables, but they have not been demonstrated for chaotic
systems. For example, Daniels _et al._ [22, 23] combinatorially fit each model
in a predefined model space using data assimilation and subsequently use
Bayesian inference to select the best model. Successful recovery of mass-
action kinetic systems for chemical reactions was demonstrated with hidden
variables using a neural network approach [37]. A recent method uses LASSO to
select predictive models for chaotic systems from a library with higher order
derivatives given a single state variable [61]. This method effectively finds
a higher-order ODE representation of the Lorenz and Rössler systems, but it is
unclear how the recovered structures relate to the ground truth models.

In this paper we present a new method to perform model selection in dynamical
systems with hidden variables. This method combines the data assimilation
technique variational annealing, which has been used to estimate parameters
when the structure of the system is known, with sparse model selection via
hard thresholding. We call this method Data Assimilation for Hidden, Sparse
Inference (DAHSI). To demonstrate that our method could identify interpretable
models for chaotic systems, we followed the philosophy of earlier works [58,
16] and demonstrated recovery of accepted parsimonious models from
experimental data and simulated time-series where the ground truth is known.
In the Results section DAHSI successfully selected a set of models for a
circuit that has Lorenz-like behaviour from experimental data of two state
variables (one hidden). One of the identified models has the same structure as
the Lorenz system. The other identified models with high AIC/BIC support
exhibit nearly indistinguishable dynamics and suggest novel terms which may
better represent the experimental circuit system. Moreover, we used ground
truth simulations of the canonical Lorenz system to study how our method
performs with varying data size and noise. In the Materials and Methods
section we describe the DAHSI algorithm for model selection with hidden
variables.

##  2 Results: Model selection for chaotic systems

###  2.1 Identification of models for the Lorenz circuit from experimental
data

The Lorenz system [42] was originally developed to forecast the weather and
has become a canonical example when developing new methods to characterize
chaotic systems. To demonstrate model selection on experimental data with
hidden variables, we considered high-quality data from the electrical circuit
in Blakely _et al._ [10] (Fig. 1(a)). This system exhibits similar structure
and behavior to the highly studied Lorenz system and is well described by
relatively simple circuit equations

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=\hat{\sigma}(y-x),$ (1)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=\hat{\rho}x-\hat{\gamma}y-\hat{\varepsilon}xz,$ (2)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=-\hat{\beta}z+\hat{\eta}xy.$ (3)

The structure of this system is similar to the Lorenz system, but in the
standard Lorenz formulation $\hat{\varepsilon}=\hat{\eta}=\hat{\gamma}$. Here,
$\mathbf{X}=(x,y,z)$ denote the voltages across the capacitors $C_{1}$,
$C_{2}$ and $C_{3}$ in the circuit (Fig. 1(a)). The measured variables are $x$
and $z$, and $y$ is unmeasured or hidden. We denote the noisy measurements of
$x$ and $z$ by $x_{e}$ and $z_{e}$, respectively, and the measurement function
$\mathbf{h}((x,z))=(x_{e},z_{e})=\mathbf{Y}$. The experimental sampling rate
is $\Delta t_{e}=80$ ns resulting in $55,000$ time points. A low-pass filter
was applied to remove high-frequency measurement error [10]. We re-scaled the
experimental time by $\Delta t=\frac{\Delta t_{e}}{1.6\times
10^{-5}\textrm{ns}}=0.005$ so that the right hand side terms of (1)-(3) are
around $\mathcal{O}(1)$. We trained our method with $N=501$ time points (Fig.
1(a)), at a sampling rate $2\Delta t=0.01$ (re-scaled). The attractor is
reasonably well sampled with 501 points (SI Appendix, Fig. 2), and we retain
the remaining data for validation.

We demonstrated model identification with hidden variables of the Lorenz-like
system ((1)-(3)) using DAHSI (Fig. 1). First, we constructed a model library
based on domain-knowledge. In this case we used monomials up to degree two in
three variables, representing $10^{9}$ possible models composed of subsets of
possible terms. From this library we generated a generic governing equation
for each variable via the linear combination of all the candidate functions
(Fig. 1(b1)). Our goal was to find a small subset of candidate functions
within the library which describe the dynamics in the data. We did not assume
that we knew the “correct” model complexity a priori, and searched for the set
of models which balance error and simplicity.

To perform model selection, we minimised a cost function composed of the
measurement error, $A_{E}$, model error, $A_{M}$, and sparse penalty,
$\lambda\|\mathbf{p}\|_{1}$ as a function of the parameters, $\mathbf{p}$ and
library of functions $\boldsymbol{\mathbf{\Theta}}$ (Fig. 1 (b), and Materials
and Methods section). The model error contains the coupling between variables,
taking advantage of the information about hidden variables in the time-history
of the measured variables. The measurement error only depends on the
measurements and measured variables estimated from the model. Model selection
is enabled through the sparse penalty which determines the number of
parameters, $p_{k,j}$, that will be active in the model or zero.

To minimize the cost function, we used variational annealing (VA) [70], a
data-assimilation technique for non-convex parameter estimation in nonlinear,
chaotic systems. The problem is highly non-convex, with many local minima, due
to the incoherence between the data and the model [52, 1]. Decreasing the
information or measurements [38] and increasing the number of terms in the
library will both increase the number of local minima (SI Appendix, Fig. 1).
VA works by varying $R_{f}$ which sets the balance between model error and
measurement error (Fig. 1(b2)). When $R_{f}=0$, only measurement error
contributes leading to a convex cost function with an easy to find global
minima. As the model is enforced by gradually increasing $R_{f}$, the
landscape increases in complexity and many local minima appear. By
initialising the search near the minima for the previous $R_{f}$ the solution
remains near the first minima found. Varying $\lambda$ leads to different
model structures or candidate models. As the penalty strength, $\lambda$,
increases, the global minima moves to 0 in a larger number of parameters (Fig.
1(b2)). Because there are many local minima, we need to choose $N_{I}=500$
random initial guesses to fully explore the landscape.

Figure 1: DAHSI model selection for the Lorenz-like system. (a) Electrical
circuit from [10], training data of measured variables $x$ and $y$, and time-
delay embedding of test data ($\tau=0.02$). (b1) Model library and generic
governing equation for each variable. (b2) Cost function as model error weight
and the sparsity constraint vary. (c) Local minima with high cost (light
grey), low cost (dark grey) and Lorenz structure (blue) as function of
$\lambda$. (d) 25 low cost models are down-selected. (e) Model structure
identified near the Pareto front. (f) Time series, error, and relative AIC for
identified models (coloured), and higher error models (grey).

The sparse-variational annealing process generates 169 candidate models, which
must be further down-selected and validated to complete the model-selection
process. We down-selected to the 25 models (SI Appendix) with a cost function
value less than $10^{-3}$ (Fig. 1(c)). In our system there is a clear gap in
cost-function value at this value, but the criteria and gap size will be
system dependent. To ensure we have the best parameter fit for each down-
selected model we performed parameter estimation via VA without sparsity
constraint.

To validate the models, we needed to estimate an initial condition for the
hidden variable $y$, for which there is no experimental data. We used an 8th
order finite difference approximation of the time derivative of $x$ for each
model structure and solve the resulting algebraic equation for $y_{0}$ (SI
Appendix). We used the dynamic equation for $x$ since all down selected models
contain $y$ but not any higher order $y$ terms. Estimation of the initial
condition for hidden variables is only possible after the candidate models are
found and must be done for the initial condition of each segment of validation
data. This procedure takes advantage of Takens’s theorem that the information
in $y$ is available in the time-delay of $x$.

Validation within the Lyapunov time ensures that the time-series do not
diverge due to the inherent sensitivity to differences in initial conditions
introduced by measurement and numerical error. All down-selected models have a
similar Lyapunov time around 0.9 time units. We considered $S=1083$ segments
of the experimental data (excluding the training set), each of length 1/4 of a
Lyapunov time to calculate the sum of the average error for each model (Fig.
1(d)). We discarded the first four points of each time segment as these points
were used to predict the initial condition for $y$. The average error for the
$s$-th time segment of the $m$-th model is defined as
$E^{s}_{av,m}=\frac{1}{2M}\sum_{i=1}^{M}(x_{i,e}^{s}-x_{i}^{s})^{2}+(z_{i,e}^{s}-z_{i}^{s})^{2}$,
where $x_{i,e}$ and $z_{i,e}$ are the $x$ and $z$ components of the
experimental data, respectively, and $i$ is the time index. The sum of all
average errors over the time segments $S$ of the $m$-th model is
$E_{av,m}=\sum_{s=1}^{S}E^{s}_{av,m}$.

The candidate models on the Pareto front (Fig. 1(d), and SI Appendix, Table
S1) best balance model complexity and error (Fig. 1(e)). We successfully
recovered the Lorenz-like structure derived by Blakely _et al._ [10], which
has the lowest average error of recovered models with 7 active terms. For
$\lambda=3.9$ the system presented in [10] is selected for 10.6% of the
$N_{I}=500$ randomly chosen initialisation. However, we have no guarantee that
this model is the "true model" for the circuit system. All models have a
similar manifold structure (Fig. 1(e)) and low error within a Lyapunov time
(Fig. 1(f)). We believe the main limitation of our prediction window is the
uncertainty introduced by the hidden variable into the parameter estimation
during VA. This uncertainty then propagates into the $y_{0}$ estimate required
for each validation data set and magnifies noise (SI Appendix, Figs. S4 and
S5). Given the difficulty in selecting between proposed chaotic models that
exhibit such similar behaviour [2, 3], the primary goal of DAHSI as a model
identification method is to generate possible models. However, we were able to
consistently identify a unique model (salmon with 11 terms, Fig. 1(f)) with
the most support using Akaike information criteria as done in [45] and
Bayesian information criteria (SI Appendix, Fig. 6), as well as identifying a
weakly supported model (gold with 10 terms, Fig. 1(f)). By generating multiple
models that lie near the Pareto front DAHSI has effectively generated
hypothesis for additional terms, which could be tested with further
experimentation.

While DAHSI identified the same equation terms as Lorenz and the circuit
formulation from [10], the parameters fit through the final step of VA are not
the same. We compare the ability to predict the experimental data with the
classical Lorenz system, the circuit formulation from [10] and the DAHSI-
recovered models, each of which have a different number of free parameters
(Table 1). We perform parameter estimation via VA for each model and use the
validation data-set described above to calculate $E_{av}$, $\Delta$AIC, and
$\Delta$BIC. Although the average error is similar for the VA-estimated
circuit formulation and all DAHSI models, the DAHSI recovered models with 10
and 11 terms have substantially more $\Delta$AIC and $\Delta$BIC support. The
classical Lorenz parameter structure, which only has 4 free parameters, is
unable to capture the dynamics of the system. The parameters estimated via VA
for the circuit model with 6 free parameters perform much better than those
estimated from first principles [10]. Notably, the parameters estimated for
the 7-term DAHSI model are very close to the parameters estimated for the
original circuit model. Further experimentation is needed to determine if the
coefficients in the $\dot{x}$ equation should be equal, $p_{1,2}=p_{1,3}$, and
if the coefficient on the $y$ term in the $\dot{y}$ equation, $p_{2,3}$ should
be positive, negative, or zero (Fig. 1(e)). The additional terms suggested by
the 10 and 11 term DAHSI recovered models are strongly supported by the
AIC/BIC calculations, but would require further experimentation to
conclusively validate. They may represent parasitic resistances or other
physical effects which have a small but real impact on the circuit dynamics
and were neglected during the original derivation by Blakely et al. [10].
Recovery of the Lorenz-like model and identification of other models with
AIC/BIC support demonstrates that DAHSI can successfully identify parsimonious
models for chaotic systems.

Table 1: Parameter estimation for the classical Lorenz formulation (4 free
parameters, $p_{1,2}=p_{1,3}$, $p_{2,3}=p_{2,7}=-p_{3,6}$); the circuit
formulation in [10] (6 free parameters, $p_{1,2}=p_{1,3}$); and the DAHSI-
recovered models.

|  |  |  | circuit formulation | DAHSI-recovered
---|---|---|---|---|---
| Term | Parameter | classical | as in [10] | estimated | 7-terms | 10-terms | 11-terms
eq. $\dot{x}$ | $1$ | $p_{1,1}$ | – | – | – | – | – | $-0.2514$
$x$ | $p_{1,2}$ | $-29.7560$ | $-12.9032$ | $-16.5369$ | $-16.9554$ | $-17.0172$ | $-17.0582$
$y$ | $p_{1,3}$ | $29.7560$ | $12.9032$ | $16.5369$ | $18.7853$ | $19.9884$ | $19.9840$
$z$ | $p_{1,4}$ | – | – | – | – | $0.1596$ | $0.1833$
eq. $\dot{y}$ | $x$ | $p_{2,2}$ | $68.5427$ | $54.2903$ | $28.0876$ | $24.3535$ | $22.6028$ | $22.6017$
$y$ | $p_{2,3}$ | $-12.8815$ | $-1.2903$ | $-0.0763$ | $0.2580$ | $0.3346$ | $0.3567$
$xy$ | $p_{2,6}$ | – | – | – | – | $-0.0906$ | $-0.0843$
$xz$ | $p_{2,7}$ | $-12.8815$ | $-14.2857$ | $-7.6252$ | $-6.7054$ | $-6.2507$ | $-6.2561$
eq. $\dot{z}$ | $z$ | $p_{3,4}$ | $-3.4168$ | $-3.8259$ | $-3.6547$ | $-3.6835$ | $-3.6954$ | $-3.6966$
$xy$ | $p_{3,6}$ | $12.8815$ | $3.4843$ | $4.3315$ | $4.8273$ | $5.1412$ | $5.1292$
| $xz$ | $p_{3,7}$ | – | – | – | – | $0.0903$ | $0.0791$
$E_{av}$ | – | – | 2165 | 319 | 10.37 | 9.7441 | 9.0995 | 9.0345
$\Delta$AIC | – | – | 5920 | 3852 | 139.315 | 73.887 | 5.758 | 0
$\Delta$BIC | – | – | 5885 | 3827 | 114.377 | 53.937 | 0.77 | 0

###  2.2 Robustness study on the simulated Lorenz system

To study the robustness of our method to varying noise and manifold sampling
we used ground-truth time series simulated from the classic Lorenz system,

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=\sigma(y-x),$ (4)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=x(\rho-z)-y,$ (5)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=-\beta z+xy,$ (6)

where $\sigma=10$, $\rho=28$, and $\beta=8/3$. We numerically simulated the
system using Runge-Kutta 4th order and a time step of $\Delta t=0.01$ and
$N=501$, producing time-series similar to the experimental data set. As in the
experimental data set, we considered $y$ to be the hidden variable. We studied
the recovery rates of DAHSI as a function of the VA tuning parameter,
$\alpha$, and found trends similar to previous work [55], (SI Appendix, Table
S3).

First, we studied the robustness of our method to measurement error modeled as
additive Gaussian noise of mean zero and varying standard deviation, $\omega$.
Therefore, the measurement function is
$\mathbf{h}(\mathbf{X})=\mathbf{X}+\mathcal{N}(0,\omega)$. We expect that
different noise instances, controlled by the random number generator seed,
will change our recovery rate due to random corruption of essential parts of
the data or overall poor manifold sampling.

We calculated recovery for 3 different standard deviations of noise with 20
noise seeds each and calculate the cumulative distribution function of the
recovery rate (Fig. 2(a)). The random noise seeds produced wide variation in
recovery rate between 10-90% for the lowest noise, indicating that the minimal
data set used here is not very robust. As the noise strength increased, the
cumulative distributions shifted left as more seeds have lower recover.
Setting $\omega=0.01$ produced a binomial distribution, with either a high
recovery rate (> 80%) (the majority of simulations), or a low recovery rate (<
15%). For $\omega=0.05$ there were some seeds with intermediate recovery
rates, more low recovery rates, and a few seeds that with a very high recovery
rate. The noise level dramatically affected the recovery rate for
$\omega=0.1$. The vast majority of simulations led to less than 10% recovery.
More than half had 0% recovery, and only one had higher than 85% recovery.

Next, we investigated how manifold sampling affected the recovery rate of our
system. We chose 3 different noise seeds, and varied the number number of time
points $N$ by increasing the length of the time-series (Fig. 2(b)). Varying
the length of the time-series changed the sampling of the manifold,
demonstrating that sampling lobe transitions is crucial for accurate model
recovery. For one seed (light blue line) the recovery was high for $N=501$
through $N=401$. There were sharp drops in recovery of $\approx 60$% and
$\approx 15$% when the data-set lost a lobe crossing in the attractor, as
happens at $N=351$ and $N=301$, respectively. Sharp drops in another seed
(dark blue line) also occurred when the sampling of the crossing between lobes
is reduced at $N=460$ and $301$. Decreased sampling of each lobe did not
appear to have as dramatic an effect ($N=401$ to $460$). The increase of
recovery rate for the dark-blue noise instance at $N=351$ suggests that
optimal sampling requires some nontrivial balance of different dynamic
regions. The specific corruption of noise instance had a big impact on how
many crossings are needed to get a high recovery as the recovery was
consistently high for one seed (cyan). These results suggest that optimal
manifold sampling to counter noise corruption would vastly improve DAHSI
performance on data sets with high noise.

Figure 2:  Robustness to noise and manifold sampling. (a) Cumulative
distribution function of the recovery rate for three noise levels and 20 noise
seeds. $\omega=0.01$ (light grey); $\omega=0.05$ (dark grey); $\omega=0.1$
(black). (b) Recovery rate on eight different manifold sizes, for three
different noise seeds (colors).

##  3 Discussion

In this paper we have presented DAHSI, a method to identify non-linear
dynamical systems from data with hidden variables. DAHSI combines variational
annealing, a data assimilation technique, with sparse thresholding. We applied
DAHSI to an experimental data set from a circuit that exhibits Lorenz-like
dynamics [10]. The outcome is a set of candidate models, including a model
with the same Lorenz-like structure derived by Blakely _et al._ from circuit
equations [10]. Two additional parsimonious models with strong support based
on AIC/BIC-based validation were also identified. The unanticipated terms
suggested by these models may represent real physical processes in the
circuit, such parasitic resistances or other factors not included in the
idealized model derivation. Through this example, we demonstrated that DAHSI
works as an effective tool for generating models and functional hypothesis
from data.

To analyze recovery and the effects of noise and manifold sampling in a system
where we know the ground truth, we studied the performance of DAHSI applied to
simulated time-series from the classical Lorenz system. Notably, we
successfully selected the ground truth model as most likely from those
generated by DAHSI using information-criteria based validation techniques (SI
Appendix. Fig. 3). Our noise studies showed recovery rates of 80% for
$\sim\mathcal{N}(0,0.01)$ and 10% or lower for $\sim\mathcal{N}(0,0.1)$.
Therefore we anticipate that the current formulation of DAHSI will have
reasonable recovery rates for noise levels $<10\%$ of the signal value.
Further robustness to noise could be achieved through integral formulations
similar to those used for sparse regression, rather than the discretized
mapping between time-points used here [57, 51, 47]. Manifold sampling impacts
recovery and we conclude that recovery is especially sensitive to sampling at
the saddle point transition between the lobes. Moreover, the noise seed used
to generate the synthetic data impacts the recovery and we suspect this is due
to random corruption of measurements from different regions of the manifold.
For chaotic systems, increasing the time of experiment will eventually ensure
robust sampling of the manifold. However, the computational time of DAHSI
scales with the length of the input time-series [31]. Therefore, we anticipate
that short bursts of time-series designed to optimally sample the manifold
would provide optimal sampling and computational efficiency. Further metrics
for analyzing the information content of our data and minimal data-
requirements for recovering models [35] would lead to optimal manifold
sampling.

One of the main benefits of a sparse model selection framework is that we
identify likely model structures while avoiding combinatorial testing and
validation of all potential models. For example, the number possible models
described three variables with monomials up to degree two is approximately
$10^{9}$. Doing parameter estimation on each of these models and validating
would be computationally intractable, taking at least $10^{5}$ processor-days
with our setup. For comparison, our entire model selection and validation
process took just over a day of computational time. Running one initialisation
of the problem and sweeping through $\lambda=2.5:0.1:5.5$ with $N=501$ (as
done in Example 2.1) took 4 hours. We parallelized simulations using
Northwestern’s High Performance Computing Cluster Quest, running about 100 of
simulations at a time, leading to a total computational time of roughly 20
hours. Performing parameter estimation without thresholding on a single model
takes between 15 second and 15 minutes, depending on model structure.
Parameter estimation on 25 down-selected models took 5 hours with our set up.
Times estimates are for a Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz processor.
In order to understand the impact of library size on a call to IPOPT, the
optimiser used in DAHSI, we tested model libraries with 7, 10, 13, 16, 19, and
30 terms (SI Appendix, Fig. 7). The computational time does not scale
monotonically with library size. Instead, we find that a library with 10 terms
can take 100 times longer to run than the library of 30 monomials. We suspect
that the variation in optimization time depends on correlations between
library functions [43], model symmetries, and other structural features.

In addition to the chaotic systems presented in the results, we have applied
DAHSI on two non-chaotic systems: on time-series data from a Lotka-Volterra-
like system with no hidden variables and on simulated time-series for a mass
action kinetics systems with hidden variables. Although DAHSI recovered
reasonable models for both systems, there are several caviats. Recover of
Lotka-Volterra required an iterative formulation (SI Appendix, Fig. 15). We
also compared DAHSI to SINDy [16] for the Lotka-Volterra system and found that
SINDy was far superior in speed when all variables are observed. Recall that a
comparison between DAHSI and SINDy is not possible for Lorenz-like circuit
system, as SINDy requires access to the unmeasured $y$ variable. The mass
action kinetic system modeled a semiconductor with two trap levels differing
by one electronic unit of charge (SI Appendix). The recovery rate for the
ground truth model was low, around 3%. Unlike chaotic systems, which are
highly non-convex, the mass-action kinetic system has a very flat cost
function due to structural parameter identifiability issues (SI Appendix,
Figs. S11-S14), [7, 29, 46, 25]. Stochastic gradient decent algorithms such as
IPOPT are known to perform poorly for flat cost functions so switching to an
optimiser designed for such systems [40] may improve recovery. Other data-
assimilation methods for parameter estimation with hidden variables such as
3D-Var, 4D-Var, Kalman filtering, and hybrid methods [8] may be more cost-
effective if VA is unnecessary to navigate to the global minimum of a highly
non-convex function.

The formulation of cost function and sparsity constraint also likely impacts
recovery. Different methods for sparse model-selection include stepwise and
all-subsets regression, ridge regression [36], LASSO [68], least angle
regression [24], and SR3 [73]. SR3 accelerates convergence and has been shown
to outperform other methods and improves performance but has an extra tuning
parameter. The parameter path for the first four methods is shown to be
different in [34] and therefore, we expect that different regularisation
methods will lead to different model identification. Comparison between
different sparsity-enforcement mechanisms within DAHSI framework could improve
recovery but may be somewhat system dependent.

We anticipate many future applications and extensions of DAHSI. The framework
for DAHSI does not have any intrinsic restrictions about the functional form
of the equations, in particular the function library need not be linear in the
unknown parameters. Variational annealing is designed to handle stochasticity
through the model error. In addition, data assimilation is commonly used for
PDE systems, including PDE discovery [19]. Therefore, we anticipate we can
apply or extend our framework to broader applications, without reformulation
as was needed in sparse-regression based frameworks for rational functions
[44], stochastic systems [14], and PDEs [56, 39]. Modifications to the
optimization methodology and further investigation of optimal data-sampling
strategies could improve the computational efficiency of DAHSI, opening up
higher dimensional problems to model selection with hidden variables.

##  4 Methods: Mathematical formulation of cost function and algorithm

The dynamics of many physical systems can be described by models with only a
few terms. Our goal is to retrieve the sparse system representation of these
type of systems given the measurements of some, but not all, of the state
variables. We consider a dynamical system with unknown governing equations

$\frac{\mathop{}\\!\mathrm{d}\mathbf{X}}{\mathop{}\\!\mathrm{d}t}=\mathbf{F}(\mathbf{X}(t),\mathbf{p}),$
(7)

where $\mathbf{X}=(x_{1},x_{2},\dots,x_{D})\in\mathbb{R}^{D}$ are the state
variables, $\mathbf{F}=(F_{1},\,F_{2},\dots,F_{D})$ are the unknown functions
that govern the dynamics of the system and $\mathbf{p}$ is a set of unknown
parameters.

For a system with hidden variables, the measurements
$\mathbf{Y}=(y_{1},y_{2},\dots,y_{L})\in\mathbb{R}^{L}$ are lower dimensional
$L\leq D$ than the underlying variables. The measurement function
$\mathbf{h}(\mathbf{X})=\mathbf{Y}$ is a known transformation of a subset of
the state variables in (86). In principle, the measurement function could map
some combination of state variables to a lower-dimension, as in
$\mathbf{h}(\mathbf{X})=x_{1}+x_{2}$. In this work we assume $\mathbf{h}$
captures Gaussian experimental noise such that,
$\mathbf{Y}=\mathbf{X}+\mathcal{N}(0,\omega)$. The measurements are taken at
$N$ equally spaced point in time between $[t_{1},\,t_{N}]$.

The function capturing the nonlinear dynamics of each state variable, $F_{k}$,
is assumed to be sparse in function space as has been done previously [33,
16]. Given a library of possible functions
$\boldsymbol{\mathbf{\Theta}}=(\theta_{1},\theta_{2},\dots,\theta_{q})$, we
can write a candidate function $\hat{F}_{k}$ as

$\hat{F}_{k}\coloneqq\hat{F}_{k}(\mathbf{X},\mathbf{p})=p_{k,1}\theta_{1}(\mathbf{X})+p_{k,2}\theta_{2}(\mathbf{X})+\cdots+p_{k,q}\theta_{q}(\mathbf{X}),$
(8)

for $k=1,2,\dots,D$. There is no inherent restriction that the functions be
linearly additive. The set of $p_{k,j}$ defines the vector
$\mathbf{p}\in\mathbb{R}^{P}$, where $P=Dq$ is the total number of unknown
parameters.

We want to estimate the unknown parameters $p_{k,j}$ and all state variables
$\mathbf{X}$ using only the measurements $\mathbf{Y}$ with the constraint that
$\mathbf{p}$ is sparse. This is equivalent to minimising the negative log
likelihood

$\begin{split}&A(\mathbf{X},\mathbf{p})=\frac{1}{N}\sum_{i=1}^{N}\|\mathbf{X}(t_{i})-\mathbf{Y}(t_{i})\|^{2}\\\
&+\frac{1}{N}\sum_{i=1}^{N-1}R_{f}\left\\{\|\mathbf{X}(t_{i+1})-\mathbf{f}(\mathbf{X}(t_{i}),\mathbf{p},\mathbf{\hat{F}})\|^{2}\right\\}+\lambda\|\mathbf{p}\|_{1}.\end{split}$
(9)

Here,
$\mathbf{f}(\mathbf{X}(t_{i}),\mathbf{p},\mathbf{\hat{F}})=\mathbf{X}(t_{i+1})$
defines the discrete time model dynamics and is obtained by discretising (86)
using a Hermite-Simpson collocation. We note that if $\lambda=0$ in 9 we
obtain the cost function used in VA. Following the statistical derivation in
[66, 26, 1], the experimental error,
$A_{E}(\mathbf{X},\mathbf{Y})=\frac{1}{N}\sum_{i=1}^{N}\|\mathbf{X}(t_{i})-\mathbf{Y}(t_{i})\|^{2}$
assumes Gaussian noise and the model error,
$A_{M}(\mathbf{X},\mathbf{p},\mathbf{\hat{F}})=\frac{1}{N}\sum_{i=1}^{N-1}\left\\{\|\mathbf{X}(t_{i+1})-\mathbf{f}(\mathbf{X}(t_{i}),\mathbf{p},\mathbf{\hat{F}})\|^{2}\right\\}$
assumes a relaxed delta function. We assume that the state at the $t_{i+1}$
depends only on the state at $t_{i}$. We assume that each element in
$\mathbf{p}$ follows a Laplace distribution with _mean_ $0$ (SI Appendix). The
details and necessary background to minimise (9) are presented in the
following sections.

###  4.1 DAHSI: Data Assimilation for Hidden Sparse Inference

Our algorithm, Data Assimilation for Hidden Sparse Inference (DAHSI), performs
model identification for chaotic systems from data with hidden variables. It
combines the data assimilation technique VA with sparse thresholding (Fig.
3(a)). The code base for DAHSI can be found at [54].

As the desired model complexity is unknown ahead of time, DAHSI sweeps through
different hard-threshold values, $\lambda$. For each $\lambda$, the cost
function (9), is minimized by iterating between VA [70, 71] and hard-
thresholding of the parameters. We chose the iterative framework over direct
incorporation of the $\ell_{1}$ penalty into the minimized cost function,
based on the results that show that least square with thresholding converges
locally, often outperforming convex variants [73, 20], and recent
demonstrations that LASSO makes mistakes early in the recovery pathway [62].

At each VA step, we minimize $A_{E}+R_{f}A_{M}$, which is 4DVar in its "weak"
formulation [66, 26], over $\mathbf{X}$ and $\mathbf{p}$ given $R_{f}$ using
IPOPT, an optimisation package that uses a gradient descent method [69]. The
state variables $\mathbf{X}^{\text{ini}}$ are initialized as $\mathbf{Y}$ for
the measured states and random values from a uniform distribution within
specified bounds for the unmeasured states. Since we expect the parameter
vector $\mathbf{p}$ to be sparse, it is initialized as
$\mathbf{p^{\text{ini}}}=0$.

Initially $R_{f}$ takes some small value $R_{f,0}=\epsilon$, as $R_{f}=0$
would lead to an unconstrained solution on the unmeasured states and
$\mathbf{p}$. At each step $\beta=0,1,2,\dots,\beta_{\max}$ of VA, $R_{f}$ is
updated to $R_{f}=R_{f,0}\alpha^{\beta}$, for $\alpha>1$. After each step
$\beta$ of VA, we enforce sparsity by applying a hard threshold, $\lambda$, to
$\mathbf{p}^{(\beta)}$. The solution,
$\\{\mathbf{X}^{(\beta)},\mathbf{p}^{(\beta)}\\}$, at each step of the VA
process is used as the initialization for the next step. We choose
$\beta_{\max}$ so that the cost function plateaus, Fig. 3(b), and our final
solution is $\\{\mathbf{X}^{\text{fin}},\mathbf{p}^{\text{fin}}\\}$. Because
there are many local minima, we run $N_{I}$ different initial guesses to fully
explore the landscape of $A_{E}+R_{f}A_{M}$. It is important to note that the
same $\lambda$ yields multiple models due to the $N_{I}$ different
initializations of the unmeasured states. For example, if we consider
$N_{I}=500$ with a fixed $\lambda=3.9$ in our Example 2.1, we find a total of
20 models (Fig. 3(b)).

To produce candidate models with varying sparsity, the entire $\beta$ sweep
with VA and thresholding is repeated for each $\lambda$. As with other model
identification methods, different $\lambda$ will yield different models (for
the same initialisation of unmeasured states). For one particular
initialisation in Example 2.1, with $\lambda=3.8$ the term $z$ is selected in
the first equation of the system. With larger $\lambda=3.9$, the term $z$ is
no longer selected (Fig. 3(c)). Although the same $\lambda$ yields multiple
models due to the difference of the initial choice of unmeasured states, as we
would expect, higher values of $\lambda$ produce models with fewer active
terms (Fig. 3(d)).

Figure 3: DAHSI Algorithm. (a) Schematic of Algorithm 1. (b) Action paths as
function of $\beta$ for $N_{I}=500$ and$\lambda=3.9$ (left). Final action
values (right) for high (light grey) and low (dark grey) action values; and
the Lorenz-like structure (blue). (c) Parameter $p_{1,4}$ in the last steps of
VA for $\lambda=3.8$ and $3.9$ (d) Model complexity as function of $\lambda$.
Algorithm 1 DAHSI Algorithm.

1:procedure DAHSI

2: Input: measurements $\mathbf{Y}$, generic model library $\mathbf{\Theta}$,
$\lambda_{max}$, $\beta_{max}$, $\alpha$

3: Calculate discrete function $\mathbf{\hat{F}}$ from $\mathbf{\Theta}$

4: for $l=1:L$ do

5: $x_{l}=y_{l}$ $\triangleright$ Fit measurements to data

6: Randomly initialise unobserved variables $\\{x_{l+1},\dots,x_{D}\\}$

7: $\mathbf{X}^{\text{ini}}=\\{x_{1},x_{2},\dots,x_{l},x_{l+1},\dots,x_{D}\\}$

8: Initialise $\mathbf{p}^{\text{ini}}=0$ $\triangleright$ Force sparsity

9: Assemble pair $\\{\mathbf{X}^{\text{ini}},\mathbf{p}^{\text{ini}}\\}$

10: $R_{f,0}=\epsilon$

11: while $\lambda<\lambda_{\max}$ do

12: for $\beta=0:\beta_{\max}$ do $\triangleright$ Variational Annealing

13: $R_{f}=R_{f,0}\alpha^{\beta}$

14: $\\{\mathbf{X}^{(\beta)},\mathbf{p}^{(\beta)}\\}$ =
$\min_{\mathbf{X},\mathbf{p}}A_{E}(\mathbf{X},\mathbf{Y})+R_{f}A_{M}(\mathbf{X},\mathbf{p},\mathbf{\hat{F}})$
$\triangleright$ Minimize via IPOPT

15: if $p^{(\beta)}_{k,j}<\lambda$ then $\triangleright$ Hard-threshold
$\mathbf{p}$

16: $p^{(\beta)}_{k,j}=0$

17: model${}^{(\lambda)}\leftarrow\mathbf{p}^{(\beta)}$ $\triangleright$ Store
models

18: $\lambda=2\lambda$ $\triangleright$ Increase $\lambda$

## References

  * [1] H. Abarbanel, Predicting the future: completing models of observed complex systems, Springer, 2013.
  * [2] L. A. Aguirre and S. Billings, Validating identified nonlinear models with chaotic dynamics, International Journal of Bifurcation and Chaos, 4 (1994), pp. 109–125.
  * [3] L. A. Aguirre and C. Letellier, Modeling nonlinear dynamics and chaos: a review, Mathematical Problems in Engineering, 2009 (2009).
  * [4] A. A. R. AlMomani, J. Sun, and E. Bollt, How entropic regression beats the outliers problem in nonlinear system identification, Chaos: An Interdisciplinary Journal of Nonlinear Science, 30 (2020), p. 013107.
  * [5] N. S. Altman, An introduction to kernel and nearest-neighbor nonparametric regression, The American Statistician, 46 (1992), pp. 175–185.
  * [6] P. Amil, M. C. Soriano, and C. Masoller, Machine learning algorithms for predicting the amplitude of chaotic laser pulses, Chaos: An Interdisciplinary Journal of Nonlinear Science, 29 (2019), p. 113111.
  * [7] J. F. Apgar, D. K. Witmer, F. M. White, and B. Tidor, Sloppy models, parameter uncertainty, and the role of experimental design, Molecular BioSystems, 6 (2010), pp. 1890–1900.
  * [8] R. N. Bannister, A review of operational methods of variational and ensemble-variational data assimilation, Quarterly Journal of the Royal Meteorological Society, 143 (2017), pp. 607–633.
  * [9] B. P. Bezruchko, D. A. Smirnov, and I. V. Sysoev, Identification of chaotic systems with hidden variables (modified bock’s algorithm), Chaos, Solitons & Fractals, 29 (2006), pp. 82–90.
  * [10] J. N. Blakely, M. B. Eskridge, and N. J. Corron, A simple lorenz circuit and its radio frequency implementation, Chaos: An Interdisciplinary Journal of Nonlinear Science, 17 (2007), p. 023112.
  * [11] S. Boccaletti, The control of chaos: theory and applications, Physics Reports, 329 (2000), pp. 103–197.
  * [12] M. Bocquet, J. Brajard, A. Carrassi, and L. Bertino, Data assimilation as a learning tool to infer ordinary differential equation representations of dynamical models, Nonlinear Processes in Geophysics, 26 (2019), pp. 143–162.
  * [13] M. Bocquet, J. Brajard, A. Carrassi, and L. Bertino, Bayesian inference of chaotic dynamics by merging data assimilation, machine learning and expectation-maximization, Foundations of Data Science, 2 (2020), pp. 55–80.
  * [14] L. Boninsegna, F. Nüske, and C. Clementi, Sparse learning of stochastic dynamical equations, The Journal of chemical physics, 148 (2018), p. 241723.
  * [15] J. Brajard, A. Carassi, M. Bocquet, and L. Bertino, Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observations: a case study with the lorenz 96 model, arXiv preprint arXiv:2001.01520, (2020).
  * [16] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences, 113 (2016), pp. 3932–3937.
  * [17] E. J. Candès and M. B. Wakin, An introduction to compressive sampling, IEEE signal processing magazine, 25 (2008), pp. 21–30.
  * [18] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton, Data-driven discovery of coordinates and governing equations, Proceedings of the National Academy of Sciences, 116 (2019), pp. 22445–22451.
  * [19] H. Chang and D. Zhang, Identification of physical processes via combined data-driven and data-assimilation methods, Journal of Computational Physics, 393 (2019), pp. 337–350.
  * [20] R. Chartrand and V. Staneva, Restricted isometry properties and nonconvex compressive sensing, Inverse Problems, 24 (2008), p. 035020.
  * [21] J. P. Crutchfield and B. McNamara, Equations of motion from a data series, Complex systems, 1 (1987), p. 121.
  * [22] B. C. Daniels and I. Nemenman, Automated adaptive inference of phenomenological dynamical models, Nature communications, 6 (2015), p. 8133.
  * [23] B. C. Daniels, W. S. Ryu, and I. Nemenman, Automated, predictive, and interpretable inference of caenorhabditis elegans escape dynamics, Proceedings of the National Academy of Sciences, 116 (2019), pp. 7226–7231.
  * [24] B. Efron, T. Hastie, I. Johnstone, R. Tibshirani, et al., Least angle regression, The Annals of statistics, 32 (2004), pp. 407–499.
  * [25] M. C. Eisenberg and M. A. Hayashi, Determining identifiable parameter combinations using subset profiling, Mathematical biosciences, 256 (2014), pp. 116–126.
  * [26] G. Evensen, Data assimilation: the ensemble Kalman filter, Springer Science & Business Media, 2009.
  * [27] H. Fan, J. Jiang, C. Zhang, X. Wang, and Y.-C. Lai, Long-term prediction of chaotic systems with machine learning, Physical Review Research, 2 (2020), p. 012080.
  * [28] A. M. Fraser and H. L. Swinney, Independent coordinates for strange attractors from mutual information, Physical Review A, 33 (1986), pp. 1134–1140.
  * [29] A. Gábor, A. F. Villaverde, and J. R. Banga, Parameter identifiability analysis and visualization in large-scale kinetic models of biosystems, BMC systems biology, 11 (2017), pp. 1–16.
  * [30] L. Gardini, C. Grebogi, and S. Lenci, Chaos theory and applications: a retrospective on lessons learned and missed or new opportunities, Nonlinear Dynamics, 102 (2020), pp. 643–644.
  * [31] J. Gondzio, Interior point methods 25 years later, European Journal of Operational Research, 218 (2012), pp. 587–601.
  * [32] M. P. Hassell, H. N. Comins, and R. M. Mayt, Spatial structure and chaos in insect population dynamics, Nature, 353 (1991), pp. 255–258.
  * [33] T. Hastie, R. Tibshirani, and J. Friedman, The elements of statistical learning: data mining, inference, and prediction, Springer Science & Business Media, 2009.
  * [34] T. Hesterberg, N. H. Choi, L. Meier, C. Fraley, et al., Least angle and l1 penalized regression: A review, Statistics Surveys, 2 (2008), pp. 61–93.
  * [35] L. S. T. Ho, H. Schaeffer, G. Tran, and R. Ward, Recovery guarantees for polynomial coefficients from weakly dependent data with outliers, Journal of Approximation Theory, 259 (2020), p. 105472.
  * [36] A. E. Hoerl and R. W. Kennard, Ridge regression: applications to nonorthogonal problems, Technometrics, 12 (1970), pp. 69–82.
  * [37] W. Ji and S. Deng, Autonomous discovery of unknown reaction pathways from data by chemical reaction neural network, arXiv preprint arXiv:2002.09062, (2020).
  * [38] N. Kadakia, The Dynamics of Nonlinear Inference, PhD thesis, UC San Diego, 2017.
  * [39] S. H. Kang, W. Liao, and Y. Liu, Ident: Identifying differential equations with numerical time evolution, arXiv preprint arXiv:1904.03538, (2019).
  * [40] V. Kantabutra and E. Zheleva, Gradient descent with fast gliding over flat regions: a first report, in IEEE 2002 28th Annual Conference of the Industrial Electronics Society. IECON 02, IEEE.
  * [41] M. B. Kennel, R. Brown, and H. D. I. Abarbanel, Determining embedding dimension for phase-space reconstruction using a geometrical construction, Physical Review A, 45 (1992), pp. 3403–3411.
  * [42] E. N. Lorenz, Deterministic nonperiodic flow, Journal of the atmospheric sciences, 20 (1963), pp. 130–141.
  * [43] N. M. Mangan, T. Askham, S. L. Brunton, J. N. Kutz, and J. L. Proctor, Model selection for hybrid dynamical systems via sparse regression, Proceedings of the Royal Society A, 475 (2019), p. 20180534.
  * [44] N. M. Mangan, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Inferring biological networks by sparse identification of nonlinear dynamics, IEEE Transactions on Molecular, Biological and Multi-Scale Communications, 2 (2016), pp. 52–63.
  * [45] N. M. Mangan, J. N. Kutz, S. L. Brunton, and J. L. Proctor, Model selection for dynamical systems via sparse regression and information criteria, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473 (2017), p. 20170009.
  * [46] N. Meshkat, M. Eisenberg, and J. J. DiStefano III, An algorithm for finding globally identifiable parameter combinations of nonlinear ode models using gröbner bases, Mathematical biosciences, 222 (2009), pp. 61–72.
  * [47] D. A. Messenger and D. M. Bortz, Weak sindy for partial differential equations, arXiv preprint arXiv:2007.02848, (2020).
  * [48] S. Mukherjee, E. Osuna, and F. Girosi, Nonlinear prediction of chaotic time series using support vector machines, in Neural Networks for Signal Processing VII. Proceedings of the 1997 IEEE Signal Processing Society Workshop, IEEE, 1997, pp. 511–520.
  * [49] E. P. Odum and G. W. Barrett, Fundamentals of ecology, vol. 3, Saunders Philadelphia, 1971.
  * [50] N. H. Packard, J. P. Crutchfield, J. D. Farmer, and R. S. Shaw, Geometry from a time series, Physical Review Letters, 45 (1980), pp. 712–716.
  * [51] Y. Pantazis and I. Tsamardinos, A unified approach for sparse dynamical system inference from temporal measurements, Bioinformatics, 35 (2018), pp. 3387–3396.
  * [52] L. M. Pecora and T. L. Carroll, Synchronization in chaotic systems, Physical Review Letters, 64 (1990), pp. 821–824.
  * [53] M. Rabinovich and H. Abarbanel, The role of chaos in neural systems, Neuroscience, 87 (1998), pp. 5–14.
  * [54] H. Ribera, DAHSI code base.  https://github.com/hribera/DAHSI, 2021.
  * [55] P. J. Rozdeba, Nonlinear Inference in Partially Observed Physical Systems and Deep Neural Networks, PhD thesis, UC San Diego, 2018.
  * [56] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Data-driven discovery of partial differential equations, Science Advances, 3 (2017), p. e1602614.
  * [57] H. Schaeffer and S. G. McCalla, Sparse model selection via integral terms, Physical Review E, 96 (2017).
  * [58] M. Schmidt and H. Lipson, Distilling free-form natural laws from experimental data, science, 324 (2009), pp. 81–85.
  * [59] B. Sivakumar, Chaos theory in hydrology: important issues and interpretations, Journal of hydrology, 227 (2000), pp. 1–20.
  * [60] J. Slingo and T. Palmer, Uncertainty in weather and climate prediction, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 369 (2011), pp. 4751–4767.
  * [61] A. Somacal, L. Boechi, M. Jonckheere, V. Lefieux, D. Picard, and E. Smucler, Uncovering differential equations from data with hidden variables, arXiv preprint arXiv:2002.02250, (2020).
  * [62] W. Su, M. Bogdan, and E. Candès, False discoveries occur early on the lasso path, The Annals of Statistics, 45 (2017).
  * [63] G. Sugihara, R. May, H. Ye, C.-h. Hsieh, E. Deyle, M. Fogarty, and S. Munch, Detecting causality in complex ecosystems, science, 338 (2012), pp. 496–500.
  * [64] G. Sugihara and R. M. May, Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series, Nature, 344 (1990), pp. 734–741.
  * [65] F. Takens, Detecting strange attractors in turbulence, in Dynamical systems and turbulence, Warwick 1980, Springer, 1981, pp. 366–381.
  * [66] O. Talagrand and P. Courtier, Variational assimilation of meteorological observations with the adjoint vorticity equation. i: Theory, Quarterly Journal of the Royal Meteorological Society, 113 (1987), pp. 1311–1328.
  * [67] Y. Tang, J. Kurths, W. Lin, E. Ott, and L. Kocarev, Introduction to focus issue: When machine learning meets complex systems: Networks, chaos, and nonlinear dynamics, Chaos: An Interdisciplinary Journal of Nonlinear Science, 30 (2020), p. 063151.
  * [68] R. Tibshirani, Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society: Series B (Methodological), 58 (1996), pp. 267–288.
  * [69] A. Wächter, An interior point algorithm for large-scale nonlinear optimization with applications in process engineering, PhD thesis, PhD thesis, Carnegie Mellon University, 2002.
  * [70] J. Ye, N. Kadakia, P. Rozdeba, H. Abarbanel, and J. Quinn, Improved variational methods in statistical data assimilation., Nonlinear Processes in Geophysics, 22 (2015).
  * [71] J. Ye, D. Rey, N. Kadakia, M. Eldridge, U. I. Morone, P. Rozdeba, H. D. Abarbanel, and J. C. Quinn, Systematic variational method for statistical nonlinear state and parameter estimation, Physical Review E, 92 (2015), p. 052901.
  * [72] K. Yeo, Model-free prediction of noisy chaotic time series by deep learning, arXiv preprint arXiv:1710.01693, (2017).
  * [73] P. Zheng, T. Askham, S. L. Brunton, J. N. Kutz, and A. Y. Aravkin, A unified framework for sparse relaxed regularized regression: Sr3, IEEE Access, 7 (2018), pp. 1404–1423.

Supplementary Information for

Model selection of chaotic systems from data with hidden variables using
sparse data assimilation

H. Ribera, S. Shirman, A. V. Nguyen and N. M. Mangan

###### Contents

  1. 1 Introduction
  2. 2 Results: Model selection for chaotic systems
    1. 2.1 Identification of models for the Lorenz circuit from experimental data
    2. 2.2 Robustness study on the simulated Lorenz system
  3. 3 Discussion
  4. 4 Methods: Mathematical formulation of cost function and algorithm
    1. 4.1 DAHSI: Data Assimilation for Hidden Sparse Inference
  5. S0 Cost function analysis
  6. S0 Time-delay embedding of used training data
  7. S0 AIC calculation for synthetic data
    1. S0.1 Initial condition choice for unmeasured y
    2. S0.2 Prediction window
  8. S0 Down-selected models
    1. S0.1 Models identified in the Pareto front edge
    2. S0.2 AIC and BIC on the 25 down-selected models
  9. S0 Action derivation
  10. S0 Computational time
  11. S0 Semiconductor
    1. S0.1 1 hidden variable
    2. S0.2 Parameter identifiability
  12. S0 Predator-Prey
  13. S0 $\alpha$ parameter in VA algorithm

##  S0 Cost function analysis

Our aim is now to explore the landscape of the cost function as to understand
the problem that we are solving and why it is very challenging. For
illustrative purposes, in the following discussion we are only considering two
dimensions of the cost function $\hat{A}=A_{E}+R_{f}A_{M}$. We use the
classical Lorenz system and take all parameters in the structure fixed and we
add two extra parameters (highlighted in red),

$\displaystyle\dot{x}$
$\displaystyle=\sigma(y-x)+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}p_{1,1}},$
(1) $\displaystyle\dot{y}$
$\displaystyle=x(\rho-z)-y+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}p_{2,1}},$
(2) $\displaystyle\dot{z}$ $\displaystyle=xy-\beta z.$ (3)

We then vary these two parameters and plot what the cost function looks like,
for three different values of $R_{f}$. The cost function that we want to
minimise, is the one that has a large $R_{f}$ value (Fig. 1, right).

Figure 1: Varying $p_{1,1}$ and $p_{2,1}$ for three different $R_{f}$.

The cost function $\hat{A}$ is highly non-convex and the task of finding its
global minima a priori is a difficult task.

##  S0 Time-delay embedding of used training data

Figure 2: Time-delay embedding of training data ($\tau=0.02$).

##  S0 AIC calculation for synthetic data

We want to find which model is the one that best represents the data synthetic
data generated (in which we added some noise $\sim\mathcal{N}(0,0.01)$). Since
we are working with chaotic systems, we only expect prediction up to the
Lyapunov time of the system. We consider 1/4 of the shortest Lyapunov time out
of all the down-selected models for the synthetic data, $t_{M}\approx 0.3$. We
use $S=300$ time series of length $t_{M}$ as our validation set, but discard
the first four points as they will be used to predict the initial condition
for $y_{0}$ (as shown in the following section). To calculate the AIC score,
we define the residual sum of squares of the $m$-th model as

$\text{RSS}_{m}=\sum_{s=1}^{S}E^{s}_{av,m}(\mathbf{Y}_{s},\mathbf{F}_{m},\mathbf{p}_{m}),$
(4)

where $\mathbf{Y}_{s}=[x_{e},z_{e}]_{s}$ is the synthetic data of the time-
series $s$, $\mathbf{F}_{m}$ the governing equations of the $m$-th model, and
$\mathbf{p}_{m}$ denotes the parameters found via parameter estimation for the
$m$-th model. $E^{s}_{av,m}$ is the average absolute error over the time-
series $s$ and is defined as

$E^{s}_{av,m}(\mathbf{Y}_{s},\mathbf{F}_{m},\mathbf{p}_{m})=\frac{1}{2M}\sum_{i=1}^{M}(x_{i,e}^{s}-x_{i}^{s})^{2}+(z_{i,e}^{s}-z_{i}^{s})^{2},$
(5)

where $x^{s}$ and $z^{s}$ denote the $x$ and $z$ component, respectively, of
the solution of the $m$-th model in the $s$ time series, found via RK4 with
$\Delta t=0.01$. $M$ denotes 1/4 of a Lyapunov time, excluding the first four
points as we have mentioned before.

Finally, we can define the AIC of the $m$-th model as

$\text{AIC}_{m}=S\log\left(\frac{\sum_{s=1}^{S}E^{s}_{av,m}(\mathbf{Y}_{s},\mathbf{F}_{m},\mathbf{p}_{m})}{S}\right)+2N_{p,m},$
(6)

where $N_{p,m}$ is the number of free parameters in the $m$-th model.

We finally re-scale by the minimum AIC value, denoted by AIC${}_{\text{min}}$,
and so $\Delta\text{AIC}_{m}=\text{AIC}_{m}-\text{AIC}_{\text{min}}$.

###  S0.1 Initial condition choice for unmeasured y

We need an initial condition for each time series to be able to simulate each
model. We have an initial condition for both $x$ and $z$ given by the
experimental data, but we do not have any information for the $y$ component.
We cannot use the VA to estimate $y_{0}$ and parameters simultaneously (which
would lead to better prediction windows see next section) because our
validation data will then have been used for training. Let us consider the 8th
order finite difference approximation of the time derivative of $x$

$\frac{\mathop{}\\!\mathrm{d}x(t)}{\mathop{}\\!\mathrm{d}t}\approx\frac{\begin{multlined}3x(t+4\Delta
t)-32x(t+3\Delta t)+168x(t+2\Delta t)-672x(t+\Delta t)\\\ +672x(t+\Delta
t)-168x(t+2\Delta t)+32x(t+3\Delta t)-3x(t+4\Delta
t)\end{multlined}3x(t+4\Delta t)-32x(t+3\Delta t)+168x(t+2\Delta
t)-672x(t+\Delta t)\\\ +672x(t+\Delta t)-168x(t+2\Delta t)+32x(t+3\Delta
t)-3x(t+4\Delta t)}{840\Delta t}$ (7)

For each model, we have that

$\frac{\mathop{}\\!\mathrm{d}x(t)}{\mathop{}\\!\mathrm{d}t}=F_{1,m}(x(t),y(t),z(t),\mathbf{p}_{m}).$
(8)

Putting (7) and (8) together we have

$\frac{-x(t+2\Delta t)+8x(t+\Delta t)-8x(t-\Delta t)+x(t-2\Delta t)}{12\Delta
t}\approx F_{1,m}(x(t),y(t),z(t),\mathbf{p}_{m}).$ (9)

We need to solve for $y(0)$. We note that for the down-selected models in
Example A in our manuscript the only terms with $y$ in the first equation in
all the models is just the first order term, so for this case this is a
particularly simple equation to solve.

The results in the synthetic data indicate that there are only four candidate
models that best represent the data. Even though the $\Delta$AIC from
incorrect models (Fig. 3, red, green and yellow lines) does not increase as we
add more time series $S$ in the calculation of AIC, we consistently pick the
correct model structure (blue line) as the one with lowest $\Delta AIC$.

Figure 3: $\Delta$AIC from the different models DAHSI found using the
synthetic data.

###  S0.2 Prediction window

We compare how the prediction window changes from having two observed
variables to having three observed variables. For noise $\omega=0.01$, having
one hidden variable (Fig. 4, top row) and using the real value of $y_{0}$,
leads to no prediction at all. However, using the estimated $y_{0}$ calculated
as in the previous section leads to a prediction window of about 3.5 Lyapunov
times. This shows that the parameter estimates and $y_{0}$ estimate are
compensating for each other. For the case of all variables observed (Fig. 4,
bottom row), we see that using the real value of $y_{0}$ leads to a prediction
window of about 6 Lyapunov times. If we estimate $y_{0}$ the prediction window
reduces to about 3.5 Lyapunov times. For a higher noise $\omega=0.1$, having
one hidden variable (Fig. 5, top row) and using the real value of $y_{0}$,
again leads to no prediction at all. Moreover, using the estimated $y_{0}$
calculated as in the previous section leads to a shorter prediction window
than for lower noise, about 1 Lyapunov time. This shows that with increased
noise amplified by hidden variables, the $y_{0}$ estimate cannot compensate
for the parameter estimate that well. For the case of all variables observed
with $\omega=0.1$ (Fig. 5, bottom row), we see that using the real value of
$y_{0}$ leads to a prediction window of about 3.5 Lyapunov times. If we
estimate $y_{0}$ the prediction window reduces to about 1 Lyapunov time.

Figure 4: Top: Prediction of the model with parameters estimated using 2
observed variables ($x$ and $z$), when using the real $y_{0}$ and the
estimated $y_{0}$ (as shown in §S0.1). Bottom: Prediction of the model with
parameters estimated using 3 observed variables, when using the real $y_{0}$
and the estimated $y_{0}$. The noise added to the synthetic data is
$\mathcal{N}(0,\omega)$, $\omega=0.01$. Figure 5: Top: Prediction of the model
with parameters estimated using 2 observed variables ($x$ and $z$), when using
the real $y_{0}$ and the estimated $y_{0}$ (as shown in §S0.1). Bottom:
Prediction of the model with parameters estimated using 3 observed variables,
when using the real $y_{0}$ and the estimated $y_{0}$. The noise added to the
synthetic data is $\mathcal{N}(0,\omega)$, $\omega=0.1$.

##  S0 Down-selected models

We present the structure of the 25 down-selected models, but we do not provide
the parameter estimation (the parameter values can be found in the code
package [54]).

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y,$ (10)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (11)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy.$ (12)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y,$ (13)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (14)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy.$ (15)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (16)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (17)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy.$ (18)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y,$ (19)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (20)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy.$ (21)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (22)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (23)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy.$ (24)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (25)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,1}+p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (26)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy.$ (27)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y,$ (28)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (29)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=0.$ (30)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y,$ (31)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (32)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,6}xy.$ (33)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y,$ (34)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (35)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=0.$ (36)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y,$ (37)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz+,$ (38)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=0.$ (39)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y,$ (40)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (41)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,6}xy.$ (42)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (43)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (44)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=0.$ (45)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y,$ (46)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (47)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,6}xy.$ (48)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (49)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (50)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy.$ (51)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (52)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (53)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,6}xy.$ (54)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (55)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (56)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,6}xy.$ (57)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (58)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (59)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy+p_{3,7}xz.$ (60)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (61)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (62)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,6}xy.$ (63)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (64)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,6}xy+p_{2,7}xz,$ (65)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy+p_{3,7}xz.$ (66)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (67)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,7}xz,$ (68)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy+p_{3,7}xz.$ (69)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (70)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,1}+p_{2,2}x+p_{2,3}y+p_{2,7}xz,$ (71)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,6}xy.$ (72)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (73)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,6}xy+p_{2,7}xz,$ (74)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy+p_{3,7}xz.$ (75)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (76)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,6}xy+p_{2,7}xz,$ (77)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy+p_{3,7}xz.$ (78)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (79)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,2}x+p_{2,3}y+p_{2,6}xy+p_{2,7}xz,$ (80)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy+p_{3,7}xz.$ (81)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z,$ (82)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,1}+p_{2,2}x+p_{2,3}y+p_{2,6}xy+p_{2,7}xz,$ (83)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,4}z+p_{3,6}xy+p_{3,7}xz.$ (84)

###  S0.1 Models identified in the Pareto front edge

Table 1: Models identified in the Pareto front in Figure 1(d) in the main
text.

|  | number of active terms
---|---|---
| Term | 6 | 7 | 8 | 9 | 10 | 11 | 12
eq. $\dot{x}$ | $1$ | 0 | 0 | -0.8112 | 0 | 0 | -0.2514 | -2.4053
$x$ | -16.5556 | -16.9554 | -16.4666 | -16.5603 | -17.0172 | -17.0582 | -17.0627
$y$ | 19.8000 | 18.7853 | 19.8120 | 16.7514 | 19.9884 | 19.9840 | 19.9862
$z$ | 0 | 0 | 0.0276 | 0.1486 | 0.1596 | 0.1833 | 1.4595
$x^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$xy$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$xz$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$y^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$yz$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$z^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
eq. $\dot{y}$ | $1$ | 0 | 0 | 0 | 0 | 0 | 0 | 0.7892
$x$ | 23.2613 | 24.3535 | 23.0763 | 27.3789 | 22.6028 | 22.6017 | 22.6061
$y$ | 0 | 0.2580 | 0 | 0 | 0.3346 | 0.3567 | 0.3298
$z$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$x^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$xy$ | 0 | 0 | 0 | -0.0922 | -0.0906 | -0.0843 | -0.3647
$xz$ | -6.3345 | -6.7054 | -6.2868 | -7.4621 | -6.2507 | -6.2561 | -6.2691
$y^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$yz$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$z^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
eq. $\dot{z}$ | $1$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$x$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$y$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$z$ | -3.6646 | -3.6835 | -3.6736 | 4.3951 | -3.6954 | -3.6966 | -3.6941
$x^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$xy$ | 5.1948 | 4.8273 | 5.2315 | -3.6660 | 5.1412 | 5.1292 | 5.1326
$xz$ | 0 | 0 | 0 | 0.0883 | 0.0903 | 0.0791 | 0.2900
$y^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$yz$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$z^{2}$ | 0 | 0 | 0 | 0 | 0 | 0 | 0
$E_{av}$ |  | 10.1693 | 9.7441 | 9.7174 | 9.6778 | 9.0995 | 9.0345 | 9.5765

###  S0.2 AIC and BIC on the 25 down-selected models

Bayesian information criteria (BIC) is defined as

$\text{BIC}_{m}=S\log\left(\frac{\sum_{s=1}^{S}E^{s}_{av,m}(\mathbf{Y}_{s},\mathbf{F}_{m},\mathbf{p}_{m})}{S}\right)+S\log(N_{p,m}).$
(85)

In the same way when we defined $\Delta$AIC in a previous section, we re-scale
by the minimum BIC value, denoted by BIC${}_{\text{min}}$, and so
$\Delta\text{BIC}_{m}=\text{BIC}_{m}-\text{BIC}_{\text{min}}$.

We will now calculate how AIC ((6)) and BIC ((85)) change as we add more time
series into the calculation. For each $S$ that we use to calculate both AIC
and BIC ($S\leq 1083$, which is the total number of time segments we have
available that are of length 1/4 of a Lyapunov time), we will pick $S$ random
time-segments to ensure that the $S$ time-segments used in the calculation are
independent samples.

For both $\Delta\text{AIC}_{m}$ and $\Delta\text{BIC}_{m}$ we are able to
consistently identify a unique model (Fig. 6). If we just look at the Pareto
front (Fig. 1(d) in the main text), one might ask if the decrease between 9
and 10 terms is meaningful. Both AIC and BIC say that it is.

Figure 6: $\Delta\text{AIC}_{m}$ and $\Delta\text{BIC}_{m}$ from the different
models DAHSI found using the experimental data in [10].

##  S0 Action derivation

We consider a dynamical system with unknown governing equations

$\frac{\mathop{}\\!\mathrm{d}\mathbf{X}}{\mathop{}\\!\mathrm{d}t}=\mathbf{F}(\mathbf{X}(t),\mathbf{p}),$
(86)

where $\mathbf{X}=(x_{1},x_{2},\dots,x_{D})\in\mathbb{R}^{D}$ are the state
variables, $\mathbf{F}=(F_{1},\,F_{2},\dots,F_{D})$ are the unknown functions
that govern the dynamics of the system and $\mathbf{p}$ is a set of unknown
parameters. Te measurements
$\mathbf{Y}=(y_{1},y_{2},\dots,y_{L})\in\mathbb{R}^{L}$ are lower dimensional
$L\leq D$ than the underlying variables.

Our goal is to find $\mathbf{X}$ and $\mathbf{p}$ that maximise the
probability $P(\mathbf{X},\mathbf{p}\;|\;\mathbf{Y},\mathbf{\hat{F}})$. We
have that [1]

$P(\mathbf{X},\mathbf{p}\;|\;\mathbf{Y},\mathbf{\hat{F}})=\int\exp\left[-A_{0}(\mathbf{X},\mathbf{Y})\right]\;\mathop{}\\!\mathrm{d}\mathbf{X}.$
(87)

Furthermore,

$A_{0}(\mathbf{X},\mathbf{Y})=-\sum_{i=1}^{N}\text{CMI}\left[\mathbf{X}(t_{i}),\mathbf{Y}(t_{i})\;|\;\mathbf{Y}(t_{0}),\dots,\mathbf{Y}(t_{i-1})\right]-\sum_{i=1}^{N-1}\log\left[P(\mathbf{X}(t_{i+1}),\mathbf{p}\;|\;\mathbf{X}(t_{i}),\mathbf{\hat{F}})\right],$
(88)

We make the following assumptions:

  1. 1.

The measurements $\mathbf{Y}$ have uncorrelated Gaussian error and that there
is no correlation between errors in measuring different quantities or at
varying time points [1];

  2. 2.

The state at the next time point depends only on the state at the current time
point, and that our model can have some error by widening the $\delta$
function it would follow otherwise using a Gaussian approximation of it [1];

  3. 3.

Each element in $\mathbf{p}$ follows a Laplace distribution with _mean_ $0$
and diversity $b$.

With assumption 1 it can be shown that

$\text{CMI}\left[\mathbf{X}(t_{i}),\mathbf{Y}(t_{i})\;|\;\mathbf{Y}(t_{0}),\dots,\mathbf{Y}(t_{i-1})\right]=\frac{1}{2\sigma_{m}^{2}}\sum_{l=1}^{L}\left(x_{l}(t_{i})-y_{l}(t_{i})\right)^{2}.$
(89)

For the second term in the sum, we need to find an expression for
$P(\mathbf{X}(t_{i+1}),\mathbf{p}\;|\;\mathbf{X}(t_{i}),\mathbf{\hat{F}})$.
Let us now focus on the $k$-th component of $\mathbf{X}(t_{i+1})$, and so our
goal is to find an expression for
$P(x_{k}(t_{i+1}),\mathbf{p}_{k}\;|\;\mathbf{X}(t_{i}),F_{k})$.

We consider the library of $q$ possible functions and the generic expression
for each equation of our model:

$\hat{F}_{k}\coloneqq\hat{F}_{k}(\mathbf{X},\mathbf{p})=p_{k,1}\theta_{1}(\mathbf{X})+p_{k,2}\theta_{2}(\mathbf{X})+\cdots+p_{k,q}\theta_{q}(\mathbf{X}),$
(90)

for $k=1,2,\dots,D$.

We can rewrite the probability we are seeking as

$P(x_{k}(t_{i+1}),\mathbf{p}_{k}\;|\;\mathbf{X}(t_{i}),F_{k})=P(\mathbf{p}_{k}\;|\;x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k})P(x_{k}(t_{i+1})\;|\;\mathbf{X}(t_{i}),F_{k}).$
(91)

Now each term in the right hand side can also be rewritten as

$\displaystyle P(\mathbf{p}_{k}\;|\;x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k})$
$\displaystyle=\frac{P(x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k}\;|\;\mathbf{p}_{k})P(\mathbf{p}_{k})}{P(x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k})},$
(92) $\displaystyle P(x_{k}(t_{i+1})\;|\;\mathbf{X}(t_{i}),F_{k})$
$\displaystyle=\frac{P(x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k})}{P(\mathbf{X}(t_{i}),F_{k})}.$
(93)

Thus, (91) becomes

$P(x_{k}(t_{i+1}),\mathbf{p}_{k}\;|\;\mathbf{X}(t_{i}),F_{k})=\frac{P(x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k}\;|\;\mathbf{p}_{k})P(\mathbf{p}_{k})}{P(\mathbf{X}(t_{i}),F_{k})}.$
(94)

We can rewrite the first therm on the right hand side in (94) as a likelihood,

$P(x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k}\;|\;\mathbf{p}_{k})=\mathcal{L}(\mathbf{p}_{k}\;|\;x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k}).$
(95)

Assuming that our next state follows a normal distribution with mean $f_{k}$
and standard deviation $\sigma^{2}$,

$\mathcal{L}(\mathbf{p}_{k}\;|\;x_{k}(t_{i+1}),\mathbf{X}(t_{i}),F_{k})=\frac{1}{\sigma\sqrt{2\pi}}\exp{\left(-\frac{\left[x_{k}(t_{i+1})-f_{k}(\mathbf{X},\mathbf{p},F_{k})\right]^{2}}{2\sigma^{2}}\right)}.$
(96)

With assumption 3, we know that each $p_{k,j}$ follows a Laplace distribution,

$p_{k,j}\sim\text{Laplace}(0,b)=\frac{1}{2b}\exp{\left(-\frac{|p_{k,j}|}{b}\right)},$
(97)

and so

$P(\mathbf{p}_{k})=\prod_{j=1}^{q}\frac{1}{2b}\exp{\left(-\frac{|p_{k,j}|}{b}\right)}.$
(98)

With this we can write (94) as

$\begin{split}P(&x_{k}(t_{i+1}),\mathbf{p}_{k}\;|\;\mathbf{X}(t_{i}),F_{k})\propto\\\
&\propto\frac{1}{\sigma\sqrt{2\pi}}\exp{\left(-\frac{\left[x_{k}(t_{i+1})-f_{k}(\mathbf{X},\mathbf{p},F_{k})\right]^{2}}{2\sigma^{2}}\right)}\prod_{j=1}^{q}\frac{1}{2b}\exp{\left(-\frac{|p_{k,j}|}{b}\right)}.\end{split}$
(99)

Note that since we are going to be minimising the action $A_{0}$ ((88)) we
forget about the constant term $P(\mathbf{X}(t_{i}),F_{k})$ in the denominator
and we just have a proportionality instead of an equality.

Note that because the $k$-th current state only depends upon the previous one,

$P\left(\mathbf{X}(t_{i+1}),\mathbf{p}\;|\;\mathbf{X}(t_{i}),\mathbf{\hat{F}}\right)=\prod_{k=1}^{D}P(x_{k}(t_{i+1}),\mathbf{p}_{k}\;|\;\mathbf{X}(t_{i}),F_{k}),$
(100)

and so, finally, we can write

$\begin{split}P(&\mathbf{X}(t_{i+1}),\mathbf{p}\;|\;\mathbf{X}(t_{i}),\mathbf{\hat{F}})\propto\\\
&\propto\prod_{k=1}^{D}\left\\{\frac{1}{\sigma\sqrt{2\pi}}\exp{\left(-\frac{\left[x_{k}(t_{i+1})-f_{k}(\mathbf{X},\mathbf{p},F_{k})\right]^{2}}{2\sigma^{2}}\right)}\prod_{j=1}^{q}\frac{1}{2b}\exp{\left(-\frac{|p_{k,j}|}{b}\right)}\right\\}.\end{split}$
(101)

Upon taking the logarithm to this expression above,

$\log(P(\mathbf{X}(t_{i+1}),\mathbf{p}\;|\;\mathbf{X}(t_{i}),\mathbf{\hat{F}}))\propto\sum_{k=1}^{D}\left\\{-\frac{\left[x_{k}(t_{i+1})-f_{k}(\mathbf{X},\mathbf{p},F_{k})\right]^{2}}{2\sigma^{2}}-\lambda\|\mathbf{p}_{k}\|_{1}\right\\}+\frac{D}{\sigma\sqrt{2\pi}}+\frac{\lambda
D}{2},$ (102)

where $\lambda=q/b$.

We have seen that (88) becomes

$A(\mathbf{X},\mathbf{p})=\frac{1}{N}\sum_{i=1}^{N}\|\mathbf{X}(t_{i})-\mathbf{Y}(t_{i})\|^{2}+\frac{1}{N}\sum_{i=1}^{N-1}R_{f}\left\\{\|\mathbf{X}(t_{i+1})-\mathbf{f}(\mathbf{X}(t_{n}),\mathbf{p},\mathbf{\hat{F}})\|^{2}\right\\}+\lambda\|\mathbf{p}\|_{1},$
(103)

which is what we wanted to show.

##  S0 Computational time

We use the Lorenz system, with all variables observed, $N=1001$ time points,
$\Delta t=0.01$, and no noise. The more terms our library
$\boldsymbol{\mathbf{\Theta}}$, the more time it takes to evaluate the cost
function associated, its Jacobian and its Hessian (Fig. 7(left)). However, due
to model symmetries and other structural features, the time to run our
algorithm does not monotonically increase with increasing number of terms in
our library. A library with 10 terms can take 100 times more to run than the
full library of 30 monomials (Fig. 7(right)).

Figure 7: 7 terms: parameter estimation. 10 terms: in blue $x^{2}$ in each
equation; in red $1$ in each equation. 13 terms: $x^{2}$ and $y^{2}$ in each
equation. 16 terms: $x^{2}$, $y^{2}$ and $z^{2}$ in each equation. 19 terms:
$1$, $x^{2}$, $y^{2}$ and $z^{2}$ in each equation. 30 terms: model selection.

##  S0 Semiconductor

We consider this semiconductor model ($T$ trap levels with two possible states
differing by one electronic unit of charge),

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=e_{n,01}y-R_{n,10}xz,$ (104)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=-e_{n,01}y+R_{n,10}xz,$ (105)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=e_{n,01}y-R_{n,10}xz.$ (106)

$x$ denotes the number of electrons in the conduction band, $y$ denotes the
number of traps with 2 electrons, and $z$ denotes the number of traps with 1
electron. We chose $e_{n,01}=0.5$ and $R_{n,10}=0.25$.

Figure 8: Dynamics from the original system (104)-(106).

Instead of using the library of all monomials in three variables up to degree
two, we know that there are only a few terms make sense physically. Our
generic model for this example is

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x+p_{1,3}y+p_{1,4}z+p_{1,5}x^{2}+p_{1,7}xz,$
(107) $\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,1}+p_{2,2}x+p_{2,3}y+p_{2,4}z+p_{2,5}x^{2}+p_{2,7}xz,$
(108) $\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,1}+p_{3,2}x+p_{3,3}y+p_{3,4}z+p_{3,7}xz.$ (109)

We first consider three observed variables, $D=L=3$. We consider a time series
of $N=101$ equally spaced time points, with $\Delta t=0.01$. The $\lambda$
sweep results in a different amount of active terms for each value. See Fig. 9
(left). Since we know the model from which our data comes from, we just want
to see if the model that has the right number of terms (highlighted in red)
corresponds to our original one, which it does.

###  S0.1 1 hidden variable

We consider two observed variables, $L=2$. We pick $x$ and $y$. We run
$N_{I}=1,000$ different initialisations. There is a question in this
particular case on how the initial guess should be picked (see Algorithm 2. We
do a $\lambda$ sweep from $\lambda=0.1$ through $\lambda=0.3$. Out of all the
1,000 different initialisations, we recover the right sparsity pattern 68
times. The optimal $\lambda=0.19$, for which we recover the right sparsity
pattern 33 times (see Fig. 10).

observed | hidden | N |  $\Delta$ t | $\beta_{\max}$ | $\lambda$ | recovery
---|---|---|---|---|---|---
2 | 1 ($z$) | 101 | 0.01 | 30 | 0.19 | 3.3%
Table 2: Recovery of the semiconductor system with one hidden variable.

Figure 9: Left: all observed variables. Right: one hidden variable.
Highlighted in red are the $\lambda$ that lead to model recovery. Figure 10:
Percentage of recovery rate for different $\lambda$ values for 1,000 different
initialisations. Initial guess for unmeasured variables is obtained through
the derivatives of the measured variables. Algorithm 2 Algorithm for picking
an initial guess for unobserved variables in the semiconductor case

1:for $d=1:(D-L)$ do $\triangleright$ loop in unmeasured variables

2: Pick at random one of the observed variables.

3: $dX_{d}\leftarrow$ Calculate gradient vector from its time series.

4: while $Z_{d}$ out of bounds do $\triangleright$ Make sure unmeasured
variable is within bounds

5: $Z_{d}(t_{1})\leftarrow$ Random initial condition for unobserved variable
within bounds.

6: for $i=1:N-1$ do

7: $Z_{d}(t_{i+1})=\Delta t\times dX_{d}(t_{i})+Z_{d}(t_{i})$

###  S0.2 Parameter identifiability

There are two main reasons of why a parameter might not be identifiable: said
parameter does not influence the model output; there is a interdependence
among different parameters, that is, one can compensate the change of one
parameter (that would influence the model output) by changing other
parameter(s) and have the output be the same. In this section, we focus on the
latter.

One way to detect pairwise interplay is by plotting contours of the cost
function versus pairs of parameters. Largely eccentric contours or _valleys_
show that the cost function is almost unchanged in one direction, and the two
parameters are highly correlated. The main drawback for our case in particular
is that we will be limited to find relationships only between pairs of
parameters instead of higher dimensional interactions.

Consider the generic model (except that the right terms are fixed –
highlighted in red; $e_{n,01}=0.5$, $R_{n,10}=0.25$)

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{1,1}+p_{1,2}x{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\,+\,e_{n,01}}y+p_{1,4}z+p_{1,5}x^{2}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\,-\,R_{n,10}}xz,$
(110) $\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{2,1}+p_{2,2}x{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\,-\,e_{n,01}}y+p_{2,4}z+p_{2,5}x^{2}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\,+\,R_{n,10}}xz,$
(111) $\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=p_{3,1}+p_{3,2}x{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\,+\,e_{n,01}}y+p_{3,4}z{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\,-\,R_{n,10}}xz.$
(112)

We will now add only two extra parameters (two of the black terms) at a time.
Each Fig. 11-14 is obtained by picking one term (parameter 1, which is the
first extra term in the system), and then study the cost function by adding
another term (parameter 2, which is the second extra term in the system). We
study this for all the possibles “parameter 2”.

Take Fig. 11. Parameter 1 here is the term $1$ in the equation
$\mathop{}\\!\mathrm{d}x/\mathop{}\\!\mathrm{d}t$, that is, $p_{1,1}$. This
extra term is fixed for all subplots. Then parameter 2 (the second extra term)
corresponds to (in order of subplots)
$p_{1,2},\,p_{1,4},\,p_{1,5},\,p_{2,1},\,p_{2,2},\,p_{2,4},\,p_{2,5},\,p_{3,1},\,p_{3,2}\,p_{3,4}$.
Figures 12-14 follow the same logic. These four figures already show
identifiability problems.

Figure 11: parameter 1 is $1$ on the
$\mathop{}\\!\mathrm{d}x/\mathop{}\\!\mathrm{d}t$ equation.

Figure 12: parameter 1 is $x$ on the
$\mathop{}\\!\mathrm{d}x/\mathop{}\\!\mathrm{d}t$ equation.

Figure 13: parameter 1 is $1$ on the
$\mathop{}\\!\mathrm{d}y/\mathop{}\\!\mathrm{d}t$ equation.

Figure 14: parameter 1 is $z$ on the
$\mathop{}\\!\mathrm{d}z/\mathop{}\\!\mathrm{d}t$ equation.

##  S0 Predator-Prey

Although it is not very common to find _pure_ predator-prey interactions in
nature, there is a classical set of data by the Hudson Bay company which
corresponds the number of snowshoe hares and Canadian lynxes trapped in
Canada, which in turn shows the relative population of both [49]. The data is
recorded yearly, so $\Delta t=1$. We use data between 1900 and 1920, thus
$N=21$. In this particular case we really do not know the dynamics behind the
system although we know that the snowshoe hare is the primary food of the
lynx. Therefore, we can assume that we have a predator-prey system, and there
is the classical Lotka-Volterra model to describe these type of dynamics. We
consider $L=D=2$ (Fig. 15(a)). We build the library of functions with all the
monomials up to degree two in two variables, and with it we construct our
generic model (Fig. 15(b)). We run our algorithm and varying $\lambda$ we
obtain a list of possible models. By looking at the corresponding AIC values
for each one, we find that the model with 7 active terms is the best one (Fig.
15(d)). We now consider that our generic model is the resulting model with 7
active terms. Again, we run the algorithm to find that the best model is one
containing only 5 terms (Fig. 15(f-h)). We iterate this process, and run the
algorithm considering the model with 5 active terms as the generic one. We
find that the best model is the one containing 4 terms (Fig. 15(i-k)). This
identified model corresponds to the Lotka-Volterra one. Once we do only
parameter estimation on it, we obtain the dynamical system shown in Fig.
15(k). We compare the original data (dashed) with the resulting model (solid),
which show an excellent match.

Figure 15: The recovery of the Lotka-Volterra system required an iterative
formulation which consisted of down-selecting relevant monomials to describe
the dynamics via AIC at the end of the variational annealing and start the
algorithm again with less terms in the generic model description.

##  S0 $\alpha$ parameter in VA algorithm

We study how the parameter $\alpha$ used to increase the value of
$R_{f}=R_{f,0}\alpha^{\beta}$ during the VA algorithm affects the recovery. We
use the class Lorenz system,

$\displaystyle\frac{\mathop{}\\!\mathrm{d}x}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=\sigma(y-x),$ (113)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}y}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=x(\rho-z)-y,$ (114)
$\displaystyle\frac{\mathop{}\\!\mathrm{d}z}{\mathop{}\\!\mathrm{d}t}$
$\displaystyle=-\beta z+xy,$ (115)

where $\sigma=10$, $\rho=28$, and $\beta=8/3$. We numerically simulate the
system using Runge-Kutta 4th order and a time step of $\Delta t=0.01$,
producing time-series similar to the experimental data set. We add some error
modeled as additive Gaussian noise of mean zero and standard deviation
$\omega=0.01$. Therefore, the measurement function is
$\mathbf{h}(\mathbf{X})=\mathbf{X}+\mathcal{N}(0,\omega)$. We consider
$N=501$, and $y$ to be the hidden variable.

As we increase $\alpha$ the recovery rate decreases, and for $\alpha\geq 1.3$
the recovery is 0% (Table 3).

Table 3: Recovery rates for varying $\alpha$. $\alpha$ | recovery rate (%)
---|---
1.1 | 93
1.2 | 87
1.25 | 20
1.3 | 0
1.4 | 0
1.5 | 0