response

# Kernel-based parameter estimation of dynamical systems with unknown
observation functions

Ofir Lindenbaum  Program in Applied Mathematics, Yale University, 51 Prospect
Street, New Haven, CT 06511, USA<EMAIL_ADDRESS>,  Amir Sagiv
Department of Applied Physics and Applied Mathematics, Columbia University,
500 West 120th Street, New York, NY 10027, USA<EMAIL_ADDRESS>,  Gal
Mishne  Halicioglu Data Science Institute, UC San Diego 9500 Gilman Drive MS
0555 SDSC 215E La Jolla, CA 92093-0555<EMAIL_ADDRESS>and  Ronen Talmon
Faculty of Electrical Engineering, Technion - Israel Institute of Technology,
Haifa 32000, Israel<EMAIL_ADDRESS>
###### Abstract.

A low-dimensional dynamical system is observed in an experiment as a high-
dimensional signal; for example, a video of a chaotic pendulums system.
Assuming that we know the dynamical model up to some unknown parameters, can
we estimate the underlying system’s parameters by measuring its time-evolution
only once? The key information for performing this estimation lies in the
temporal inter-dependencies between the signal and the model. We propose a
kernel-based score to compare these dependencies. Our score generalizes a
maximum likelihood estimator for a linear model to a general nonlinear setting
in an unknown feature space. We estimate the system’s underlying parameters by
maximizing the proposed score. We demonstrate the accuracy and efficiency of
the method using two chaotic dynamical systems - the double pendulum and the
Lorenz ’63 model.

*OL and AS contributed equally to this work.

## Lead Paragraph

The purpose of many experimental designs is to measure a quantity of interest
by observing a dynamical system and comparing the observations to a known
model. This procedure can be difficult when the system is chaotic and can be
even more challenging when the correspondence between the measurements and the
model, the observation function, is unknown. For a single experiment, i.e.,
when the observed data is a single time-series, learning the unknown
observation function is not an option. Instead, we construct a kernel-based
score in a way that is agnostic to the unknown observation function. We can
derive a maximum likelihood estimator for identifying the observed dynamical
system’s parameters by maximizing that score. The intuition behind our
approach is that even though the map between the coordinates/model and the
observations is unknown, the dynamics of the data convey enough information on
the dynamics of the model with the true parameter, thus facilitating an
informed parameter estimation procedure. We propose two optimization schemes
for maximizing our score. Finally, we demonstrate that our method can
accurately estimate the governing parameters for two chaotic dynamical systems
from complex and high-dimensional data.

##  1\. Introduction

Consider a common situation in experimental sciences - an experiment is
designed to measure a quantity of interest by observing a dynamical system and
comparing the observations to a known model. But can this measurement be
performed when the model is complex and chaotic? Furthermore, is this
procedure possible when the correspondence between the measurement and the
model, the observation function, is unknown?

Figure 1. The schematic settings of our problem. We are given the observations
$y(t)$ and a mechanism to generate $x(t;\omega)$ for every $\omega\in\Omega$
(an ODE). What is the true underlying parameter $\omega^{*}$ driving $y(t)$?

For example, it is straightforward to estimate the gravitational free
acceleration $g$ by observing a pendulum; the angle $x(t)$ of a pendulum of
length $\ell$ varies periodically according to the harmonic oscillator
ordinary differential equation (ODE) $\ddot{x}(t)=-(g/\ell)x$. By solving the
ODE, $g$ may be estimated using $g=\nu^{2}\ell$, where the frequency $\nu$ can
be directly observed from the pendulum’s oscillations. But can such a
measurement scheme be applied to the chaotic double pendulum, where no easily
observable parameter like the frequency $\omega$ exists?

The double pendulum example illustrates a more general class of problems; see
Fig. 1. In an experiment, an observed signal $y(t)$ is related to its
governing model $x(t;\omega^{*})$ by specific yet unknown parameters (or
parameter-vector) $\omega^{*}$ and an unknown observation function $G$, i.e.,

$y(t)=G(x(t;\omega^{*});\zeta)\,,$

where $\zeta$ is a noise source. The purpose of this study is to estimate the
system’s parameters $\omega^{*}$ among all possible parameters $\omega$ in a
parameter space $\Omega$, using the observation $y(t)$ and the general model
$x(t;\omega)$. Critically, we note that even though the map $\omega\mapsto
x(t;\omega)$ is known, the map $x\xmapsto{G}y$ is unknown to us. Therefore, we
only know “half” of the forward map $\omega\mapsto y(t)$. Since the forward
map is unknown, this problem does not fit into the usual definition of inverse
problems [3, 61]. Conversely, since we only observe a single experiment and do
not have a lot of data, it is not straightforwardly amenable to standard
machine learning methodology (see Sec. 6.1 for details).

Figure 2. The schematics of the proposed solution. For every $\omega\in\Omega$
a kernel $K^{x}(\omega)$ is computed. This Kernel is compared to the
observation kernel $K^{y}$, and the estimated $\hat{\omega}$ is chosen to
maximize their similarity score. The hypothesised $\omega$ values are either
predetermined (Algorithm 1) or dynamically determined using an optimization
scheme (Algorithm 2).

We propose a kernel-based approach to estimate the system’s parameters
$\omega^{*}$. We first study the case of a linear observation function $G$. A
maximum-likelihood estimation of $\omega^{*}$ then yields a maximization
problem for a normalized variant of the cross-covariance between the
observations and the model. To carry this idea to the general nonlinear case,
we “lift” both the observations and the model to an infinite-dimensional
Hilbert space (feature space [36, 55]). In the feature space, the two signals
are again linearly dependent. By constructing kernels for the observations
$y(t)$ and for the model $x(t;\omega)$, a covariance-like score in the feature
space is computed (see (15)) and maximized to estimate the system’s
parameters. By applying our method (Algorithms 1 and 2) to two examples of
chaotic dynamical systems - the double pendulum and the Lorenz system - we
demonstrate empirically that maximizing the kernel-based score indeed yields
an accurate estimate for $\omega^{*}$.

The application of the so-called kernel trick to generalize the linear notion
of covariance has been used for various statistical tasks such as kernel
principle component analysis (PCA), kernel canonical-correlation analysis
(CCA), and the Hilbert-Schmidt Independence criteria [6, 30, 32, 47, 54].
Kernels were also used in this context of kernel density estimators [63] or
for extracting latent variables from multi-modal observations as in [44, 42,
53, 72]. While resembling to some nonlinear kernel statistical problems on the
one hand, and to some machine learning and model discovery problems on the
other hand [4, 8, 12, 25, 23, 15, 18, 69], we note that the problem of
parameter learning under an unknown observation function is stated here, to
the best of our knowledge, for the first time. Consequently, our kernel-based
score does not seem to appear in the kernel methods literature; see discussion
in Section 3.

The remainder of this paper is organized as follows. Sec. 2 presents the
problem in formal terms. In Sec. 3 we derive our method initially for the
linear case, and then to the general nonlinear case. Sec. 4 presents the main
algorithms of this paper - the search-based Algorithm 1 and the optimization-
based Algorithm 2. The applications of our approach to the double pendulum and
to the Lorenz system are presented in Sec. 5. Finally, we discuss potential
applications of the method and its relationship to previous studies on model
discovery, inverse problems, and kernel methods in Sec. 6.

##  2\. Problem Formulation

Consider a parametric family of autonomous ordinary differential equations
(ODE)

(1)
$\left\\{\begin{array}[]{l}\dot{x}(t;\omega)=f(x;\omega)\,,\qquad\omega\in\Omega\subseteq\mathbb{R}^{m}\,,\\\
x(0;\omega)=x_{0}(\omega)\in\mathbb{R}^{d}\,,\end{array}\right.$

where $\Omega\subseteq\mathbb{R}^{m}$ is a convex set of possible parameters
and $f$ is sufficiently smooth such solutions are unique and exist globally.
The dynamics $x(t;\omega)$ are therefore completely determined by a fixed
vector of parameters $\omega^{*}\in\Omega$. Assume that $\omega^{*}$ is
unknown and that we do not observe $x(t;\omega^{*})$, but only a measurement
$y(t)\in\mathbb{R}^{D}$ for some dimension $D$. This observation can be viewed
as a noisy lifting of $x(t;\omega^{*})$ from the latent space $\mathbb{R}^{d}$
to the ambient observation space $\mathbb{R}^{D}$ by an unknown and possibly
noisy map, i.e.,

(2) $y(t)=G(x(t;\omega^{*});\zeta)\,,\qquad
G:\mathbb{R}^{d}\times\mathcal{Z}\to\mathbb{R}^{D}\,,$

where $\mathcal{Z}$ is some manifold in which $\zeta(t)$ is a stationary
random process with $\delta$ auto-correlation and the observation function
$G(\cdot,0)$ guarantees identifiability of $\omega$; see [64] and Section 4.3
for details .111It is not essential that $G$ is defined on all of
$\mathbb{R}^{d}$, but just on
$\cup_{\omega\in\Omega}\,\\{~{}x(t;\omega)~{}|~{}t\geq 0\\}$. Furthermore, in
our analysis and examples usually $\mathcal{Z}$ is either the space of the
model $\mathbb{R}^{d}$ or the observation space $\mathbb{R}^{D}$. For example,
if $x(t;\omega)$ describes the trajectory of a ballistic projectile in
$\mathbb{R}^{3}$, its video will embed this trajectory in $\mathbb{R}^{D}$,
where $D$ is the number of pixels in each video frame. As a practical matter,
we will further assume that $y(t)$ is measured in discrete times
$\\{t_{j}=(j-1)\Delta t\\}_{j=1}^{N}$ for some $\Delta t>0$. The main problem
of this paper can now be formally stated:

######  Problem.

Given

  1. (1)

A single observed time series $\\{y(t_{j})\\}_{j=1}^{N}$, defined by (2) with
unknown $G$ and $\omega^{*}$, and

  2. (2)

A solution $x(t;\omega)$ to the ODE (1) for all $t\geq 0$ and
$\omega\in\Omega$,

find the vector of underlying parameters $\omega^{*}$.

######  Remark 1.

Uncertainty in the observations may stem from several sources – modelling
misspecification, numerical errors, measurement noise, nuisance variables in
the experiment etc. The introduction of randomness in (2) is a modelling
decision aimed to capture all of these uncertainty sources.

We start by considering the simplified linear variant of (2), in which we can
derive a maximum likelihood estimator of $\omega^{*}$. In the general
nonlinear case, we use the kernel trick to map the phase-space coordinates
$x(t;\omega)$ and the observations $y(t)$ into a Hilbert space (feature space)
where the linear approach can be employed again.

##  3\. Derivation

###  3.1. The linear case - a maximum likelihood approach

It is instructive to first consider (2) where $G$ is linear in $x$ and
additive with respect to a Gaussian noise term, i.e.,

(3) $\bar{y}(t_{j})=A\bar{x}(t_{j};\omega^{*})+\zeta_{j}\,,\qquad
j=1,...,N\,,$

where $\bar{x}(t;\omega)=x(t;\omega)-\mathbb{E}_{\tau}x(\tau;\omega)$ and
$\bar{y}(t)=y(t)-\mathbb{E}_{\tau}y(\tau)$ are centered, $A\in
M_{D,d}(\mathbb{R})$, and for each $1\leq j\leq N$ the term $\zeta_{j}$ is
drawn iid from $\mathcal{N}(0,\sigma^{2}{\rm I})$ for some $\sigma>0$.

######  Remark 2.

We center the observations $y$ and model coordinates $x$ since the choice of
origin in either $\mathbb{R}^{D}$ and $\mathbb{R}^{d}$ is arbitrary from a
modelling/physics perspective. See more on the role of such invariances in
Section 4.3. In practice, the time-averages should be replaced by their
empirical estimates, e.g., $\mathbb{E}_{\tau}y(\tau)\approx
N^{-1}\sum_{j}y(t_{j})$.

To estimate $\omega^{*}$ from the observations $\\{y(t_{n})\\}_{n=1}^{N}$, we
use a maximum likelihood (ML) estimator [1]. The ML estimator is defined as

$\hat{\omega}_{\rm ml}=\arg\max\limits_{\begin{subarray}{c}\omega\in\Omega\\\
{\bf A}\in M_{D,d}(\mathbb{R})\end{subarray}}{\rm
Prob}_{\zeta_{1},\ldots,\zeta_{N}}(\bar{y}_{1},\ldots,\bar{y}_{N}|\omega,A)\,.$

Since the $\rm log$ function is monotonic increasing, we can replace the
likelihood function by log-likelihood to exploit the independence of the
normal $\zeta_{j}$’s to get

$\displaystyle\hat{\omega}_{\rm ml}$
$\displaystyle=\arg\max\limits_{\begin{subarray}{c}\omega\in\Omega\\\ {\bf
A}\in M_{D,d}(\mathbb{R})\end{subarray}}\log{\rm
Prob}_{\zeta_{1},\ldots,\zeta_{N}}(\bar{y}_{1},\ldots,\bar{y}_{N}|\omega,A)$
(4) $\displaystyle=\arg\max\limits_{\begin{subarray}{c}\omega\in\Omega\\\ {\bf
A}\in
M_{D,d}(\mathbb{R})\end{subarray}}-\frac{\sum\limits^{N}_{n=1}\|\bar{y}(t_{n})-A\bar{x}(t_{n},\omega)\|_{2}^{2}}{2\sigma^{2}}-\frac{DN}{2}\log(2\pi\sigma^{2})\,.$

Since $D$, $N$, $\sigma$, and the observations $y(t_{j})$ are independent of
$\omega$ and $A$, we can simplify the objective function on the right-hand-
side as follows:

$\displaystyle\hat{\omega}_{\rm ml}$
$\displaystyle=\arg\min_{\begin{subarray}{c}\omega\in\Omega\\\ A\in
M_{D,d}(\mathbb{R})\end{subarray}}\sum^{N}_{n=1}\|\bar{y}(t_{n})-A\bar{x}(t_{n},\omega)\|_{2}^{2}$
$\displaystyle=\arg\min_{\begin{subarray}{c}\omega\in\Omega\\\ A\in
M_{D,d}(\mathbb{R})\end{subarray}}\|AX(\omega)\|_{F}^{2}-2\langle
AX(\omega),Y\rangle+\|Y\|_{F}^{2}$ (5)
$\displaystyle=\arg\min_{\begin{subarray}{c}\omega\in\Omega\\\ A\in
M_{D,d}(\mathbb{R})\end{subarray}}\|AX(\omega)\|_{F}^{2}-2\langle
AX(\omega),Y\rangle\,,$

where $X(\omega)$ and $Y$ are matrices whose $j$-th columns are
$\bar{x}(t_{j};\omega)$ and $\bar{y}(t_{j})$, respectively, $\langle
B,C\rangle=\sum_{i,j}B_{i,j}C_{i,j}$ is the Frobenius inner product on
$M_{D,d}(\mathbb{R})$ and $\|B\|_{F}=\langle B,B\rangle^{1/2}$ is the
Frobenius norm. Next, we show that one can also normalize (5) by $\|AX\|_{F}$.
To see that, fix $\omega\in\Omega$ and $O\in M_{D,d}(\mathbb{R})$ such that
$A=\lambda O$, $\|OX\|_{F}=1$, and $\lambda\in\mathbb{R}$. Then, by direct
differentiation in $\lambda$

$\|AX(\omega)\|_{F}^{2}-2\langle
AX(\omega),Y\rangle=\lambda^{2}-2\lambda\langle OX(\omega),Y\rangle\,,$

is minimized when $\lambda=\langle OX,Y\rangle$. Hence,

$\hat{\omega}_{\rm ml}=\arg\min_{\begin{subarray}{c}\omega\in\Omega\\\
\|OX\|_{F}=1\end{subarray}}\langle OX(\omega),Y\rangle^{2}-2\langle\langle
OX(\omega),Y\rangle OX(\omega),Y\rangle\,.$

By setting $O=A\|AX\|_{F}^{-1}$, we get that the maximum likelihood estimator
is

(6) $\hat{\omega}_{ml}=\arg\max_{\begin{subarray}{c}\omega\in\Omega\\\ A\in
M_{D,d}(\mathbb{R})\end{subarray}}~{}\frac{{\Large{<}}AX(\omega),Y{\Large{>}_{F}}^{2}}{\|AX(\omega)\|_{F}^{2}\cdot\|Y\|_{F}^{2}}\,,$

where, since $Y$ is a constant matrix, we divide by $\|Y\|_{F}^{2}$ so that
the argument on the right hand side is always $\leq 1$.

###  3.2. Maximum likelihood in the nonlinear settings - a kernel approach

In the general model (2), the observations $y(t)$ do not depend linearly on
the model coordinates $x(t)$ as in (3). Rather, the two depend nonlinearly via
$G$, see (2). If we knew $G$, the linear maximum-likelihood approach (6) could
be applied to $y$ and $G(x(\cdot;\omega^{*});\zeta)$ in $\mathbb{R}^{D}$. Even
though we do not know $G$, the relation $\eqref{eq:Gdef}$ implies that $x$ and
$y$ are linearly dependent under nonlinear transformations
$\psi:\mathbb{R}^{d}\to\mathbb{R}^{D}$ and
$\phi:\mathbb{R}^{D}\to\mathbb{R}^{D}$ (feature maps), respectively, i.e.,

(7) $\phi(y(t))=\psi(x(t;\omega^{*});\zeta)\,.$

In what follows, we assume that the noise term $\zeta$ is again Gaussian and
additive, i.e., $\phi(y(t))=\psi(x(t;\omega^{*}))+\zeta$ where
$\zeta\sim\mathcal{N}(0,\sigma^{2}{\rm I})$. In this case, one possible set of
feature maps is simply $\phi(y)=y$ and $\psi(x)=G(x)$. However, we show in the
numerical experiments that our algorithm can estimate $\omega^{*}$ even for
non-Gaussian and non-additive noise sources, see Sec. 5.3.

######  Remark 3.

Seemingly, the nonlinear model (7) is more restrictive than its linear
counterpart (3), absent of the freedom to choose the linear transformation
$A$. Given the maps $\phi$ and $\psi$, however, $A$ is “absorbed” into the
definitions of $\phi$ and $\psi$.

Using the same maximum-likelihood argument of Sec. 3.1, the nonlinear model
(7) yields the following estimator of $\omega^{*}$ (compare with (6))

(8a)
$\hat{\omega}=\arg\max_{\begin{subarray}{c}\omega\in\Omega\end{subarray}}~{}\frac{\langle\bar{\Psi}(\omega),\bar{\Phi}\rangle^{2}_{F}}{\|\bar{\Psi}(\omega)\|_{F}^{2}\cdot\|\bar{\Phi}\|_{F}^{2}}\,,$
where (8b)
$\bar{\Psi}_{\cdot,j}(\omega):\,=\psi(x(t_{j};\omega))-\frac{1}{N}\sum\limits_{i=1}^{N}\psi(x(t_{i};\omega))\,,\qquad\bar{\Phi}_{\cdot,j}:\,=\phi(y(t_{j}))-\frac{1}{N}\sum\limits_{i=1}^{N}\phi(y(t_{i}))\,.$

We further note that

(9) $\|\bar{\Phi}\|_{F}^{2}=\rm
Tr(K^{y}H)\,,\qquad\|\bar{\Psi}(\omega)\|_{F}^{2}=\rm
Tr(K^{x}(\omega)H)\,,\qquad H_{ij}:\,=\delta_{ij}-\frac{1}{N}\,,$

where the Gram matrix $K^{y}$ is defined by

(10) $K^{y}_{ij}:\,=k^{y}(y(t_{i}),y(t_{j}))\,,\qquad
k^{y}(y,y^{\prime}):\,=\langle\phi(y),\phi(y^{\prime})\rangle\,,$

for every $1\leq i,j\leq N$, and $K^{x}$ and $k^{x}$ are defined analogously
to (10); see the details in Appendix A.

We can therefore rewrite (8a) as

(11)
$\hat{\omega}=\arg\max_{\omega\in\Omega}\frac{\langle\bar{\Psi}(\omega),\bar{\Phi}\rangle_{F}^{2}}{\|K^{x}(\omega)H\|_{F}^{2}\|K^{y}(\omega)H\|_{F}^{2}}\,.$

In practice, we do not know what the the maps $\phi$ and $\psi$ are, and
consequently cannot compute their corresponding kernels and Gram matrices
(10). However, the kernels $k^{x}$ and $k^{y}$ are both Mercer kernels, i.e.,
their Gram matrices (10) are always positive semi-definite [52], they define
an inner product on some (infinite-dimensional) reproducing kernel Hilbert
Space $\mathcal{H}$.222The kernels $k^{x}$ and $k^{y}$ are generally not inner
products in the original spaces $\mathbb{R}^{d}$ and $\mathbb{R}^{D}$,
respectively, since they are not linear in either of their components. This
observation suggests that one should apply a well known heuristic known as the
“kernel trick”, where rather than estimating the feature map $\psi$ from an
infinite-dimensional function space, one chooses the kernels $k^{x}$ and
$k^{y}$. The choice of kernels reflects our understanding and prior knowledge
of the data and what makes two samples $x(t_{i})$ and $x(t_{j})$ similar (and
similarly for $y$). The kernel trick here seems necessary, since we do not
know the observation function $G$. If we were to reconstruct $G$, then we
could have used $\psi=G$, $\phi={\rm Id}$, and $\mathcal{H}=\mathbb{R}^{D}$,
see e.g. [11]. In our setting, however, recovering $G$ might be nearly
impossible, since we have only a single output time-series $y(t)$ and do not
know the parameter vector $\omega^{*}$.

We revisit (11) with the kernel trick in mind. The numerator
$\langle\bar{\Psi}(\omega),\bar{\Phi}\rangle_{F}^{2}$ is equal to
$\text{Tr}^{2}(C_{xy}(\omega))$, where
$C_{xy}(\omega)=\bar{\Psi}(\omega)\bar{\Phi}^{\rm T}$ is the cross-covariance
operator. This operator norm cannot be computed, since we only choose the
kernels $k^{x}$ and $k^{y}$ and not the feature maps $\phi$ and $\psi$.
Nonetheless, we can use the kernels $k^{x}$ and $k^{y}$ to estimate ${\rm
Tr}(C_{xy}^{\rm T}C_{xy})=\|C_{xy}(\omega)\|_{F}^{2}$ as a surrogate to ${\rm
Tr}^{2}(C_{xy})$. The norm $\|C_{xy}(\omega)\|_{F}^{2}$ is known as the
Hilbert Schmidt independence criterion (HSIC), due to Gretton et al. [32, 33],
and can be estimated empirically by ${\rm Tr}(K^{x}(\omega)HK^{y}H)$. We
therefore use the ${\rm HSIC}(\omega)$ to define the following realizable
proxy of objective (8):

(12) $\hat{\omega}=\arg\max_{\omega\in\Omega}~{}\frac{{\rm
HSIC}(\omega)}{\|\bar{\Psi}(\omega)\|_{F}^{2}\cdot\|\bar{\Phi}\|_{F}^{2}}=\arg\max_{\omega\in\Omega}\frac{{\rm
Tr}(K^{x}(\omega)HK^{y}H)}{\|K^{x}(\omega)H\|_{F}^{2}\|K^{y}(\omega)H\|_{F}^{2}}\,.$

The HSIC is an indicator for the dependence of $\phi$ and $\psi$ (or $x$ and
$y$) as random variables. As $y(t)$ is determined by $x(t;\omega^{*})$ up to a
noise term, we expect that ${\rm HSIC}(\omega^{*})$ would express a high
statistical dependency. However, HSIC is also maximized as the covariance-
matrix of $x$ (or of $\Psi$ in the nonlinear settings) is maximized,
regardless of $C_{xy}$. We therefore normalize the HSIC estimator by the norm
of the standard deviation matrix
$\|\left(\bar{\Psi}(\omega)\bar{\Psi}(\omega)^{\rm T}\right)^{1/2}\|_{F}^{2}$.
This latter matrix is estimated in the nonlinear Hilbert space settings by
$\|K^{x}(\omega)\|_{F}^{2}$, which leads to the right-hand side of (12). The
empirical estimator of HSIC in the nominator on the right side is known to
have a bias of ${\mathcal{O}}(1/N)$ [32].

Different perspective: The estimator $\omega_{\rm ml}$ in (12) can also be
viewed in terms of kernel density estimators (KDE). In [63], the authors study
non-rigid shape correspondence in three-dimensional objects. Given $N$ points
$x_{i}$ and $y_{i}$ on two respective deformations of the same shape, the goal
is to find a permutation $\pi$ on $1,\ldots N$ such that each $y_{i}$
corresponds to $x_{\pi(i)}$ in the deformed shape. Letting $K^{x}(\pi)$ be the
Gram matrix of the permutated $x_{i}$’s, the term ${\rm Tr}(K^{x}(\pi)K^{y})$
can be understood as the KDE estimator of the joint probability of the points
under the permutation $\pi$. Indeed, in these settings the maximum likelihood
estimator of $\pi$ among all permutations maximizes ${\rm
Tr}(K^{x}(\pi)K^{y})$, much as in Section 3 of this paper. Note that since the
parameter estimated in [63] is a permutation, $\|K^{x}(\pi)\|_{F}^{2}$ is
constant (independent of $\pi$) and therefore the normalization is not needed.

In [47], the authors use the kernel trick to solve the nonparametric CCA
problem, i.e., identify nonlinear mappings $\psi$ and $\phi$ that maximize the
correlation between $\psi(x)$ and $\phi(y)$ (as in the numerator of (8a)).
They show that the solution can be expressed using the singular values of the
joint probability density of $X$ and $Y$, and use kernels to estimate this
joint density in a nonparametric fashion. In the context of our problem, this
solution can be used after $\omega^{*}$ is estimated for comparing the
trajectories in a low dimensional representation.

##  4\. Proposed method

###  4.1. Exhaustive search approach

In light of the analysis of Sec. 3, our first proposed method, Algorithm 1, is
a straightforward numerical application of (12). Assume for simplicity that
$\Omega$ is a box in $\mathbb{R}^{m}$, i.e.,
$\Omega=\prod_{i=1,\ldots,m}[a_{i},b_{i}]$ where $b_{i}>a_{i}$ for all
$i=1,\ldots,m$. We define the Gram matrices for $x$ and $y$ by

(13)
$K^{y}_{i,j}:\,=\exp\left(-\frac{\|y(t_{i})-y(t_{j})\|_{2}^{2}}{\varepsilon_{y}}\right)\,,\qquad\varepsilon_{y}>0\,.$

(14)
$K^{x}_{i,j}(\omega):\,=\exp\left(-\frac{\|x(t_{i};\omega)-x(t_{j};\omega)\|_{2}^{2}}{\varepsilon_{x}}\right)\,,\qquad\varepsilon_{x}>0\,.$

Then, following (12), we propose the score

(15) $s(\omega)=\frac{\mathop{\rm
Tr}\nolimits(K^{x}(\omega)HK^{y}H)}{\|K^{x}(\omega)H\|_{F}\cdot\|K^{y}H\|_{F}}\,,$

where $H_{ij}=\delta_{ij}-N^{-1}$, which Algorithm 1 maximizes on a pre-
determined set of grid-points in $\Omega_{\rm search}\subset\Omega$.

As noted in Sec. 3.2, since we do not know the feature maps $\phi$ and $\psi$,
one needs to choose the kernel $k^{x}$ and $k^{y}$ (or their Gram matrices
$K^{x}$ and $K^{y}$) to use (12). There are many possible choices of kernels
which reflect different notions of affinities between samples of $x(t)$ (and
of $y$), many of which might have worked well for estimating $\omega^{*}$, see
e.g., [36]. In this work we choose the widely popular Gaussian kernel
$k^{x}(x(t_{i}),x(t_{j}))=\exp(-\|x(t_{i})-x(t_{j})\|^{2}/\varepsilon_{x})$
(and respectively for $y$), for two main reasons: first, the Gaussian kernel
is translation invariant, i.e., $k^{x}(x,x^{\prime})=k(x-x^{\prime})$, which
ensures we capture only relative changes in the data, and not absolute values.
Second, the exponential decay of the Gaussian kernel attenuates the effect of
large distances. For the model coordinates $x(t)$, this makes intuitive sense
because the governing ODE (1) is local. For the observation coordinates,
attenuating large distances counteracts the spuriously large pairwise
distances which tend to appear in $\ell^{2}(\mathbb{R}^{D})$ for $D\gg 1$, see
e.g., [19, 35]. Finally, the Gaussian kernel has expressivity properties that
ensure its ability to capture polynomial-order nonlinear dependencies between
different times, but it is by no means the only possible choice. For detailed
analysis and alternative kernel choices, see [59] and the references therein.
A key takeaway from this work is that even though other tailored kernels can
be designed for better performance, the standard choice of the Gaussian kernel
yields good results in our numerical experiments.

Given $\\{y(t_{n})\\}_{n=1,\ldots,N}$.

1: Compute the $y$ kernel (13).

2: Choose $\Omega_{\rm search}\subset\Omega$ to be a finite Cartesian grid in
$\Omega\subseteq\mathbb{R}^{n}$.

3: for each $\omega\in\Omega_{\rm search}$ do

4: Solve (1) for $x(t;\omega)$.

5: Compute the $x$ kernel (14).

6: Compute the score $s(\omega)$, see (15).

7: end for

8: return $\hat{\omega}=\arg\max\limits_{\omega\in\Omega_{\rm
search}}s(\omega)$.

Algorithm 1 Kernel search-based method for estimating $\omega^{*}$

The effectiveness of Gaussian kernels strongly relies on proper tuning of the
kernels’ bandwidths $\varepsilon_{x}$ and $\varepsilon_{y}$. These parameters
directly affect the feature maps induced by the kernels. At one extreme,
setting $\varepsilon_{x}$ or $\varepsilon_{y}$ too large would result in
kernels that approach the all-ones matrices, i.e., where all samples are
equally affine. At the other extreme, as $\varepsilon_{x}\rightarrow 0$ (resp.
$\varepsilon_{y}$), the kernel $K$ approaches the identity matrix, i.e., all
affinities between different samples are neglected. Here, we use a max-min
measure suggested in [39] where the scale is set to

(16) $\varepsilon_{y}=\underset{j}{\max}[\underset{i,i\neq
j}{\min}(||y(t_{i})-y(t_{j})||^{2})]\,,\qquad i,j=1,...N\,,$

and analogously for $\varepsilon_{x}$. The max-min approach guarantees that
for each data-point, $K$ expresses a non-negligible affinity with at least one
other point. Moreover, this scheme generates kernels $K^{x}$ and $K^{y}$ that
are invariant to scaling in the ambient observation space $\mathbb{R}^{D}$.
Several other methods have been proposed for tuning $\varepsilon$ and can
certainly be used, see for example [39, 71, 43, 58]. The third tunable
parameter in Algorithm 1 is the grid size $|\Omega_{\rm search}|$. The choice
of a predeteremined grid $\Omega_{\rm search}$ affects Algorithm 1 in two
ways:

  1. (1)

Accuracy: For any estimator $\hat{\omega}$ of $\omega^{*}$, define the
estimation error

(17) ${\rm
Err}\,\omega_{j}=\frac{\|\hat{\omega}_{j}-\omega_{j}^{*}\|_{2}}{\omega_{j}^{*}}\,,\qquad
j=1,..,m\,.$

Generally in Algorithm 1, the grid does not necessarily include $\omega^{*}$,
i.e., $\omega^{*}\not\in~{}\Omega_{\rm search}$, and so
$\hat{\omega}\neq\omega^{*}$. Therefore, even in the best case scenario where
Algorithm 1 returns the closest grid point to $\omega^{*}$, the error is
typically bounded from below in terms of $\Delta\omega$. This error estimation
applies also to the case of a single parameter, as the average estimation
error (17) over many experiments scales like $\Delta\omega$, where

$\Delta\omega:\,=\min\\{|\omega_{i}-\omega_{j}|~{}~{}{\rm
s.t.}~{}\omega_{i},\omega_{j}\in\Omega_{\rm search}\,,~{}i\neq j\\}\,,$

is the spacing of the grid. To see that, suppose for simplicity that
$\Omega=[0,\Delta\omega]$ and that $\omega^{*}$ is drawn uniformly at random
from $\Omega$, i.e., the probability density function of $\omega^{*}=y$ is
$p(y)=\Delta\omega^{-1}$. Then by partitioning $\Omega$ to two halves, one
closer to the grid point $\omega=0$ and the other to the grid point
$\omega=\Delta\omega$, we have that

$\mathbb{E}_{\omega^{*}}\,{\rm dist}(\omega^{*},{\rm
grid})=\frac{1}{\Delta\omega}\left(\int\limits_{0}^{\Delta\omega/2}\omega^{*}\,d\omega^{*}+\int\limits_{\Delta\omega/2}^{\Delta\omega}(\Delta\omega-\omega^{*})\,d\omega^{*}\right)=\frac{\Delta\omega}{4}\,.$

A similar estimate holds for any dimension $m\geq 1$ of $\Omega$.

  2. (2)

Efficiency: In the multi-dimensional case $\Omega\subseteq\mathbb{R}^{m}$,
fine grids are computationally prohibitive since their size scales
exponentially with $m$. Large grids are especially an issue when either (i)
solving the underlying dynamical system (1) is computationally expensive, (ii)
the length of the time series $N$ is large. In the latter case, the
computation of the kernel (14), which requires evaluating all pairwise
distances, requires $\mathcal{O}(N^{2})$ operations. The cost of the kernel
computation can be reduced using methods such as $k$-sparse graph [66] which
enjoys a reduced complexity of ${\mathcal{O}}(NlogN+Nk)$, where $k$ is the
number of nearest neighbors used for building the graph.

###  4.2. Optimization approach

To overcome both of the accuracy and the efficiency problems of Algorithm 1,
we would like to replace the exhaustive grid search with an optimization
scheme to solve the following problem

(18) $\begin{array}[]{l}\text{maximize }s(\omega)=\frac{\mathop{\rm
Tr}\nolimits(K^{x}(\omega)HK^{y}H)}{\|K^{x}H(\omega)\|_{F}\cdot\|K^{y}H\|_{F}}\\\
\text{over }\omega\in\mathbb{R}^{m}\\\ \text{subject to }R(\omega)\leq
0\,,\end{array}$

where $R$ can be chosen to be any function such that $R(\omega)\leq 0$ if and
only if $\omega\in\Omega$.333It can often be the case that for $\omega$ values
outside of $\Omega$, the underlying ODE does not yield a well-posed solution.
For example, the harmonic oscillator $\ddot{x}(t)+\omega x=0$ for $\omega<0$
yields exponentially growing and ill-posed solutions. Critically, the
optimization problem (18) is in general non-convex. This is an inherent
feature of our problem and should not depend on the specific solution
strategy. Indeed, let $\tilde{s}(\omega)$ be any convex cost function for
which $\omega^{*}=\arg\max\tilde{s}(\omega)$. Since $\tilde{s}$ depends on
$\omega$ indirectly through $f(x;\omega)$ (see (1)) and since there is no
unique way to express the dependence of $f$ on its parameters, one can find an
equivalent parametrization of the ODE (1) such that the resulting
$\tilde{s}(\omega)$ is no longer convex. Therefore, the standard form of our
ODE of interest need not result in a convex parametrization of
$\tilde{s}(\omega)$. We note that this non-convexity property is already true
for the linear model (6).

Given $\\{y(t_{n})\\}_{n=1,\ldots,N}$.

1: Compute the kernel for $y$ (13).

2: Choose $\Omega_{\rm init}\subset\Omega$ at random.

3: for each $\omega_{j}\in\Omega_{\rm init}$ do

4: Solve the optimization problem (18) to find $\hat{\omega}_{j}$ using
interior point optimization initialized at $\omega_{j}$, where the kernel
$K^{x}(\omega)$ is given by (14) and $x(t;\omega)$ are the solutions of (1).

5: end for

6: return $\hat{\omega}=\arg\max\limits_{j}s(\hat{\omega}_{j})$.

Algorithm 2 Kernel optimization-based method for estimating $\omega^{*}$

To solve the constrained, nonlinear and non-convex problem (18), our proposed
method, Algorithm 2, has “two layers” of optimization. At the heart of
Algorithm 2, we use the Interior-point Algorithm (IPA), a solver for nonlinear
constrained optimization problems [16, 17], as the optimization scheme in Step
4 of Algorithm 2.444In our simulations, we used MATLAB’s implementation of the
IPA method. The IPA method first takes a Newton step by attempting to solve a
linear approximation of the problem (18), then, a gradient step is performed
using a trust region [65]. Since the problem is not convex (see details in the
next paragraph), we repeatedly initialize the IPA method at random points
$\Omega_{\rm init}\subset\Omega$, and then chooses the optimal result over all
iterations. In our experience, the multiple initialization mechanism improves
our chances of finding the global maximum; see experimental results in Sec. 5.

Since the IPA routine in Algorithm 2 dynamically samples $\omega\in\Omega$
values, the overall number of samples does not scale exponentially with the
dimension, and the accuracy is not limited by the grid spacing. Algorithm 2 is
therefore computationally cheaper than Algorithm 1. The key parameter in
comparing the two is the number of times the ODE (1) is solved and $K^{x}$
(14) is computed. In Algorithm 1, this is exactly the grid size $|\Omega_{\rm
search}|$. In Algorithm 2, the number of evaluations of $x(t;\omega)$ is the
number of optimization initializations $|\Omega_{\rm init}|$ multiplied by the
number of iterations in each optimization process. In the examples considered
in this study, $|\Omega_{\rm init}|$ was kept orders of magnitude smaller than
the $|\Omega_{\rm search}|$ without much loss of accuracy, and the number of
iterations in each optimization process is usually below $20$; see e.g., Fig.
6. Algorithm 2 therefore suggests an avenue to solve (18) efficiently and
accurately even as the dimension $m$ of the parameter space $\Omega$
increases.

######  Remark 4.

Another practical implementation aspect of Algorithms 1 and 2 is the numerical
solution method of the ODE (1). Throughout this paper we used the standard
fourth order Runge-Kutta method [37]. Other numerical methods for ODEs can be
used; see [38] for a more thorough discussion.

###  4.3. Degeneracies and Identifiability

When is the problem of estimating $\omega^{*}$ unsolvable? If
$x(t;\omega)=x(t;\omega^{*})$ for some $\omega\neq\omega^{*}$, these two
parameters are indistinguishable in terms of the resulting dynamics. Moreover,
if $G(x(t;\omega))=G(x(t;\omega^{*}))$, then the experiment/observation of the
dynamical system cannot distinguish between the parameters. These are extreme
cases for obstructions of identifiability and observability, topics which have
been studied in both the statistics and control literature, see e.g., [48, 64]
and the references therein. The system (1) combined with the observation
function (2) is said to be unidentifiable if the problem of estimating
$\omega^{*}$ is not solvable, regardless of the method of solution. In loose
terms, such $G$ then corresponds to a flawed experiment which is not designed
to estimate $\omega$.

Even for an identifiable and observable system, the fact that the observation
function $G$ is unknown limits our use of standard inference techniques. In
what follows we wish to discuss the limitation of our method even for
identifiable systems: suppose that neither $x$ nor $G$ are degenerate (as
functions of $\omega$ and $x$, respectively), but that the Gaussian kernels in
(13) and (14) are degenerate. To explore the effect of these degeneracies, we
will consider the case where $G$ is an
$\ell^{2}(\mathbb{R}^{d}\to\mathbb{R}^{D})$ isometry and noiseless, i.e.,
$\|G(x_{1})-G(x_{2})\|_{2}=\|x_{1}-x_{2}\|_{2}$ for any
$x_{1},x_{2}\in\mathbb{R}^{d}$. This discussion highlights some of the
considerations that go into designing $K^{x}$ and $K^{y}$.

######  Lemma 1.

Let $y(t)=G(t)$ where $G$ is an $\ell^{2}(\mathbb{R}^{d}\to\mathbb{R}^{D})$
isometry. Then

$s(\omega^{*})=1=\max_{\omega\in\Omega}s(\omega)\,.$

The other $\omega\in\Omega$ values where $s(\omega)=1$ are precisely those for
which $x(\omega)=Tx(\omega^{*})$ for some $\ell^{2}(\mathbb{R^{d}})$ isometry
$T$.

See the proof in Appendix B. Intuitively, Lemma 1 implies that while
$s(\omega)$ is not uniquely maximized in this case, it is “reasonably
degenerate”, in the sense that its only other maximizers (when $G(x)=x$) are
$\omega\in\Omega$ for which the trajectory $x(t;\omega)$ is isometric to
$x(t;\omega^{*})$. For example, if $x$ is a scalar, it would mean that such
degeneracies are only translations and reflections.

This is a fundamental consideration in the kernel design - if
$\omega,\omega^{*}\in\Omega$ manifest in isometric observations/trajectories,
then one cannot distinguish between them using the observations. The
assumption that underlies our choice of the $\ell^{2}$ norm in the kernel (14)
now becomes apparent - it expresses the kind of degeneracies we wish to allow.
Indeed, different choices of norms would yield different equivalence classes
of parameters in $\Omega$.

Finally, we make a crucial note regarding our choice of the Gaussian kernel.
In the proof of Lemma 1 we show that, in ideal settings, $s(\omega)$ is
maximized only when $\Delta=K^{x}(\omega)-K^{x}(\omega^{*})=0$. However, since
the kernel (14) is exponentially decaying with $\|x_{i}-x_{j}\|_{2}^{2}$, the
entries $\Delta_{ij}$ are practically null for most far-away indices (times)
$i$ and $j$. Hence, even local in time degeneracies are sufficient to cause
estimation error. As noted above, the choice of $\varepsilon_{x}$ and
$\varepsilon_{y}$ determines how “local” the respective kernels are.

##  5\. Experimental Results

To test Algorithms 1 and 2, we apply them to two classical chaotic dynamical
systems - the double pendulum and the Lorenz system.

###  5.1. Double Pendulum

The double pendulum consists of two pendulums, one attached to the end of the
other. The Lagrangian of this system is given by [40]

$L=\frac{1}{2}\left(m_{1}+m_{2}\right)l_{1}^{2}\dot{\theta}_{1}^{2}+\frac{1}{2}m_{2}l_{2}^{2}\dot{\theta}_{2}^{2}+m_{2}l_{1}l_{2}\dot{\theta}_{1}\dot{\theta}_{2}\cos(\theta_{1}-\theta_{2})\,,$

where $l_{j},m_{j}$, and $\theta_{j}(t)$ are the length, mass, and clock-wise
angle from the negative $y$ direction of the $j$-th pendulum, for $j=1,2$.
From this Lagrangian, a fourth-order system of Euler-Lagrange ODEs for
$(\theta_{1},\theta_{2},\dot{\theta}_{1},\dot{\theta}_{2})$ can be derived;
see e.g., [56] and Appendix C.

Figure 3. Left: A frame of an artificial double pendulum video. Right: The
chaotic trajectory $(x_{2}(t),y_{2}(t))$ of the bottom bob. Figure 4.
Application of Algorithms 1 to estimate the double pendulum parameters
$\omega=(l_{1},l_{2},m_{2})$. The horizontal and vertical coordinates $x_{i}$
and $y_{i}$ of the two bobs ($i=1,2)$ for the true parameter $\omega^{*}$
(black) and using $\hat{\omega}$ when estimated using Algorithm 1 (dashed
blue).

The parameters of the system are $\omega=(m_{2},l_{1},l_{2})$, where we set
$m_{1}=1$ since the motion of the double pendulum depends only on the ratio
$m_{2}/m_{1}$; see Appendix C. We describe the dynamics of the double pendulum
using the Cartesian coordinates of the two bobs $x=(x_{1},y_{1},x_{2},y_{2})$
(with a slight abuse of notations). The dynamics of the pendulum are observed
through a synthetic video with a frame rate of $\Delta t=~{}0.01$. Each frame
embeds the model in a high dimensional space $y(t_{i})\in\mathbb{R}^{D}$,
where $D=~{}171\times 217=37,107$ is the number of pixels in each frame; see
Fig. 3.

Figure 5. Application of Algorithms 1 and 2 to estimate the double pendulum
parameters $\omega=(l_{1},l_{2},m_{2})$. Top: Box plots of normalized errors
(17) of estimated parameters $l_{1},l_{2}$ and $m_{2}$ for the double pendulum
based on $20$ test cases. Results base on Algorithm 1 (left) and Algorithm 2
(right). Bottom: Box plots of normalized errors for the linear algorithms 1
and 2, (19) and (20), respectively. In all box plots, red lines represent the
medians, the blue box represents the $25$-th and $75$-th quantiles, the black
whiskers represent $\pm 2.7$ standard deviations, and any data point beyond
this distance is considered as an outlier and marked by a red plus.

To evaluate the proposed approach, we generate $20$ instances of a double
pendulum with the parameters $m_{2},l_{1},l_{2}$ each drawn iid from a uniform
distribution in the interval $[1,10]$ and with initial conditions randomly
drawn independently from $\mathcal{N}(0,0.5)$, where $N=1,000$. We then
generate from each instance a movie of total run-time $T=N\cdot\Delta t=10$.
First, we apply Algorithm 1 to each movie and search over $20$ values of each
parameter (where the initial conditions are known), i.e., a search-grid of
size $|\Omega_{search}|=20^{3}=8,000$ and $\Delta\omega=0.45$. In Fig. 4, we
see that the estimated parameters yield pendulum trajectories nearly
indistinguishable from the true ones (Fig. 4). Overall, the median normalized
parameter estimation error (17) of Algorithm 1 is $7\%$, $8\%$, and $5\%$ for
$m_{2}$, $l_{1}$, and $l_{2}$, respectively. We repeat the same experiment for
the optimization-based Algorithm 2 with $|\Omega_{\rm init}|=1,000$, where the
overall median normalized errors are $3.6\%$, $2.3\%$, and $2.4\%$ for
$l_{1}$, $l_{2}$, and $m_{2}$, respectively; see box-plots in Fig. 5. In this
experiment, as well as in the Lorenz system (Sec. 5.3), we assume that the
initial conditions are known. This is a realistic assumption if the experiment
is controlled by the observer, or if the initial conditions can be estimated
independently. If the initial conditions are not known, one option is to
estimate them as additional parameters using Algorithm 1 or 2.

To provide baselines for the proposed algorithms we use two linear variants of
Algorithm 1. Both variants are based on the same grid search procedure as in
Algorithm 1 but using a different score. The score for the first linear
estimator (Linear 1) is defined by

(19) $s^{\rm
lin_{1}}(\omega)=\frac{\|\bar{X}\bar{Y}^{T}\|_{F}}{\|\bar{X}\|_{F}\|\bar{Y}\|_{F}},$

where $\bar{X}$ and $\bar{Y}$ are the centered versions of $X$ and $Y$. For
the second linear estimator, we use Gram matrices instead of Gaussian kernels
for $k^{y}$ and $k^{x}$ , thus the score is defined by

(20) $s^{\rm lin_{2}}(\omega)=\frac{\mathop{\rm
Tr}\nolimits(G^{x}(\omega)HG^{y}H)}{\|G^{x}H(\omega)\|_{F}\cdot\|G^{y}H\|_{F}},$

where $G^{x}=X^{T}X$ and $G^{y}=Y^{T}Y$. As evident from the box-plots of
Linear 1 and 2 (Fig. 5), the performance of the linear variants of Algorithm 1
are inferior compared with their kernel based counterparts. Specifically, the
overall median error of Linear 1 is higher than Algorithm 1 by a factor of
$6.5$. For Linear 2 this factor is $9.7$.

Next, to demonstrate the efficacy of the optimization-based Algorithm 2, we
record the values attained by the IPA method throughout its iterations. As can
be seen in Fig. 6, most of the IPA runs converge after 10 iterations, and all
of them converge in less than 30 iterations (results not shown). Some runs
converge to parameters with a relatively low score. Nevertheless, since we
choose the IPA run of the highest score, we find that the overall error of
Algorithm 2 is low.

Figure 6. $10$ runs of the IPA optimization method in Algorithm 2 for the case
of the double pendulum; see Sec. 5.1. The score (15) vs. the optimization
step. Here the initial parameters for the IPA are drawn uniformly at random
from $[0.9,1.1]\cdot\omega^{*}$. The solid line indicates the best IPA run in
terms of the score $s(\omega)$.

In light of the discussion in Sec. 4.3, it is worth asking whether measuring
the estimation error (17) is the right way to estimate one’s performance in
estimating a dynamical system’s parameters. Certainly, it might be that the
problem is inherently ill-posed and that $\omega^{*}$ is non-identifiable if
another parameter $\omega\in\Omega$ produce observations similar to $y(t)$.
Along this reasoning, we propose another metric to measure our performance,
the prediction error. This metric compares the normalized mean square error
(MSE) of the predicted trajectory $x(t;\hat{\omega})$ from the true one
$x(t;\omega^{*})$

(21) $\frac{\int\limits_{t=0}^{t_{\rm
f}}\|x(\tau;\hat{\omega})-x(\tau;\omega^{*})\|_{2}^{2}\,d\tau}{\int\limits_{t=0}^{t_{\rm
f}}\|x(\tau;\hat{\omega})\|_{2}^{2}\,d\tau}.$

Figure 7. Box plots of normalized prediction error (21) for the double
pendulum experiment, as in Fig. 5, for Algorithms 1 and 2.

We evaluate the prediction error of Algorithms 1 and 2 by comparing the true
pendulum trajectory to the estimated trajectory using the same $20$ synthetic
videos used for the box-plots in Fig. 5. The results (see Fig. 7) demonstrate
that the median prediction errors for Algorithm 1 and 2 are $0.27\%$ and
$0.2\%$ respectively. Moreover, both methods attain a prediction error smaller
than $9\%$ in all of the simulations. In this experiment, the prediction error
(21) is comparable to the estimation error (17) squared, as could be expected
from their respective definitions.

###  5.2. Estimation error decreases with signal length

Figure 8. Application of Algorithm 1 to the double pendulum systems with
varying total run-time $T_{\rm f}$; see Sec. 5.2. Mean normalized estimation
error (17) for the three parameters $l_{1}$, $l_{2}$ and $m_{2}$ as a function
of the total length $T_{\rm f}$, with error-bars of one standard deviation.

In many experiments, the overall run time $T_{\rm f}=N\cdot\Delta t$ of the
experiment is a tunable parameter. How does $T_{\rm f}$ affects the estimation
error (17)? In particular, we want to test the intuition according to which
longer signals provide more information. For each $T_{\rm f}=1,2,\ldots,10$ we
draw 100 samples of $l_{1},l_{2},m_{2}$, each with the uniform distribution on
the set $\\{1.5,2,2.5,\ldots,7\\}$. To each set of these parameters we apply
Algorithm 1 as in the previous section. The mean estimation error (17) decays
as a function of $T_{\rm f}$, see Fig. 8. Furthermore, it is evident that the
standard deviation generally decreases with $T_{\rm f}$, which reinforces the
intuition of more information in longer signals.

We remark, however, that we do not expect the accuracy to increase with
$T_{\rm f}$ for all dynamical systems, especially not for systems with
attractors or limiting cycles. Consider for example $\dot{x}(t)=-x$ with
$x(0)=\omega\in\mathbb{R}_{+}$. Since $x(t;\omega)=~{}\omega e^{-t}$, then
$x(t;\omega)\approx 0$ for $t\gg 1$ independently of $\omega$. Since the
kernels we use in Algorithms 1 and 2 weight all times equally, the longer the
signal is, the more weight that is given to times $t\gg 1$ where $\omega^{*}$
is practically non-identifiable; see discussion in Sec. 4.3.

###  5.3. Lorenz ’63 system

Figure 9. (a) Lorenz system (22) with $\rho=60$, $\sigma=20$, $\beta=~{}8/3$.
(b) The observed $y(t)$, see (23). (c) The coordinates of $x(t;\omega^{*})$
(solid, black) and its noisy variant $x(t;\omega^{*})+\zeta$ (grey) vs. time
for the true parameter. For each coordinate we present the trajectory using
the estimated parameters based on the clean signal $\hat{\omega}_{c}$ (dots,
blue) and noisy signal $\hat{\omega}_{n}$. (dashes, red)

Consider the Lorenz (or Lorenz ’63 [46]) System of ODEs

(22) $\begin{split}\dot{x}_{1}(t)&=\sigma(x_{2}-x_{1})\,,\\\
\dot{x}_{2}(t)&=x_{1}(\rho-x_{3})-x_{2}\,,\\\ \dot{x}_{3}(t)&=x_{1}x_{2}-\beta
x_{3}\,,\end{split}$

where $\sigma,\rho,\beta\in\mathbb{R}_{+}^{3}$ are the model parameters. The
Lorenz System is nonlinear, and for certain parameters and initial conditions
it is chaotic; see Fig. 9(a). To embed this system in a high-dimensional
space, we follow the nonlinear transformation introduced in [15, 18]. Let
$u_{j}\in\mathbb{R}^{128}$ be the $j$-th order Legendre polynomial evaluated
on $128$ uniformly-spaced points in $[-1,1]$, and let555This precise form of
$G$ is not particularly important, only that it is truly nonlinear and
incorporates all of $x$’s coordinates. Other observation functions were tested
to yield similar result with our algorithm (results not shown).

(23)
$y(t)=G(x(t)):\,=u_{1}{x}_{1}(t)+u_{2}{x}_{2}(t)+u_{3}{x}_{3}(t)+u_{4}{x}_{1}(t)^{2}+u_{5}{x}_{2}(t)^{2}+u_{6}{x}_{3}(t)^{2}\,,$

where the $\omega$ notations were omitted for brevity; see Fig. 9(b). We apply
Algorithm 2 to $y(t)$ to estimate $\omega^{*}$. Our algorithm’s estimation of
these parameters leads to nearly indistinguishable low-dimensional
trajectories $x(t)$; see Fig. 9(c).

We then consider a noisy observation function

(24) $\mathcal{G}(x(t);\zeta)=G(x(t)+\zeta)\,,$

where $G$ is given by (23) and $\zeta\sim\mathcal{N}(0,\sigma^{2}I_{3})$ is
three-dimensional normally distributed with $\sigma=15$. Note that in this
case, the resulting noise in the observation is neither additive nor Gaussian.
The estimated parameters produce comparable trajectories; see Fig. 9(d).

Figure 10. Box plots of normalized errors (17) for estimated parameters
$\sigma$ and $\rho$ for the Lorenz ’63 model (22) given the observation
function (23). Top left: Algorithm 1. Top right: Algorithm 2. Red lines
represent the medians and whisker bars indicate the $25$-th and $75$-th
quantiles. Bottom left: Oracle estimator using the unknown observation
function $G$. Bottom right: estimation of $3$ parameters using Algorithm 2.

To systematically study the performance of our proposed approach, we repeat
the following experiment $20$ times: we set $\beta=~{}8/3$, draw
$\sigma\in[15,25]$ and $\rho\in[40,80]$, the initial conditions $x_{1}(0)$,
$x_{2}(0)$, and $x_{3}[0]$ from $[0,1]$, all uniformly at random. For each set
of parameters and initial conditions we solve the Lorenz ODE (22) and sample
the solution with $\Delta t=0.01$ and $T_{\rm f}=N\cdot\Delta t=10$, where
$N=1,000$. For each of these $20$ instances we estimate the parameters
$\sigma$, and $\rho$ using Algorithms 1 and 2. For Algorithm 1 we use a grid
with $50$ values of each parameter, i.e., a search-grid of size $|\Omega_{\rm
search}|=50^{2}=2,500$, with $\Delta\sigma=0.2$ and $\Delta\rho=0.8$. The
median normalized estimation error (17) of Algorithm 1 is $12\%$, and $1.2\%$
for $\sigma$, and $\rho$, respectively. We repeat the same experiment for the
optimization-based Algorithm 2 with $|\Omega_{\rm init}|=100$, where the
overall median normalized errors are $7.5\%$, and $1.6\%$ for $\sigma$, and
$\rho$, respectively; see box-plots in Fig. 10. To provide a baseline, we
further evaluate the performance of an “Oracle” estimator. We define the
Oracle estimator as the linear maximum likelihood estimator (4) when $G$ is
known. Using this estimator we find the parameter $\omega\in~{}\Omega_{\rm
search}$ that minimizes the square difference between $y$ and $G(x(\omega))$.
The box-plots of the oracle estimation appears in the right side of Fig. 10,
the median normalized errors are $6.4\%$ and $0.3\%$ for $\sigma$ and $\rho$,
respectively. Next, we use Algorithm 2 to estimate all three parameters of the
Lorenz system $\sigma,\rho$ and $\beta$. To ensure observability, we use
$\beta=8/3$ and draw $\sigma\in[15,25]$ and $\rho\in[40,80]$, the initial
conditions $x_{1}(0)$, $x_{2}(0)$, and $x_{3}(0)$ from $[0,1]$, all uniformly
at random.666We observe that for $\beta\neq 8/3$, many different sets of
parameters produce nearly indistinguishable trajectories, thus obstructing
observability, see discussion in Sec. 4.3. We use $|\Omega_{init}|=100$
($\beta$ is initialized by drawing from $[2,3]$), and obtain an overall median
normalized of $7.4\%$, $5.3\%$ and $5.1\%$ for $\sigma$, $\rho$ and $\beta$,
respectively; see bottom right box-plots in Fig. 10.

Figure 11. Estimating $\sigma$ in the Lorenz equation (22) from the
observations (23) where all other parameters are fixed and known. Parameter
estimation error (17) vs. number of grid points $|\Omega_{\rm grid}|$; see
Alg. 1.

The accuracy of Algorithm 1 is limited by the grid resolution of $\Omega_{\rm
search}$. The introduction of an optimization scheme in Algorithm 2 removes
this limitation, but brings a host of other accuracy issues, e.g., the non-
convexity of the optimization problem (12). But, practical considerations
aside, what is the inherent accuracy of the score (15), independently of the
search or optimization method? To answer this question, we repeated the test
of Algorithm 1 on the noiseless observation (23) of the Lorenz system (22),
but when only the parameter $\sigma$ is unknown. This allows us to
substantially refine the grid and increase the size of $\Omega_{\rm grid}$
from $10$ to $10^{4}$. For each grid-resolution we average over $100$
simulations and observe that the median estimation error (17) decreases by
more than an order of magnitude; see Fig. 11. Whether this trend saturates at
some point, or conversely the error vanishes as $|\Omega_{\rm
grid}|\to\infty$, remains an interesting open question.

###  5.4. Model coordinates and phase space coordinates

The results above highlight an important distinction between the model
coordinates $x(t;\omega)$ and phase space coordinates. From a dynamical
systems perspective, it might seem as if our model coordinates $x(t;\omega)$
in the double pendulum example consist of only partial data - the double
pendulum system (Appendix C) is a fourth-order ODE, with two angles and two
angular velocities, while we use only the position coordinates. Therefore, our
model coordinates $x(t;\omega)$ do not specify a single point in phase space.
The key observation here is that we observe a time series, and not a single
point in time. We claim that temporal derivatives are implied by the time
series. This can be understood heuristically, as the differences between
subsequent times are indicative of the velocities. In the spirit of Taken’s
theorem, delay coordinates $(x(t_{1};\omega),\ldots,x(t_{n};\omega))$ can
reconstruct the phase space, especially in this instance where we only “drop”
two coordinates, see e.g., [26].

Figure 12. Box plots of the normalized errors (17) for estimating $\sigma$ and
$\rho$ in the Lorenz system (22) based on the partial observations function
$y_{\text{partial}}$, see (25).

To test the hypothesis that time-series of partially observed data can be
sufficient for parameter estimation purposes, we applied Algorithm 2 for the
case where only two out of the three coordinates of the Lorenz model (22) are
available; we use the following partial observation function

(25)
$y_{\text{partial}}(t)=G(x(t)):\,=u_{1}{x}_{1}(t)+u_{2}{x}_{2}(t)+u_{4}{x}_{1}(t)^{2}+u_{4}{x}_{2}(t)^{2},$

where ${x}_{1}(t)$ and ${x}_{2}(t)$ are coordinates of the Lorenz ’63 system
(22) and $u_{i}\in\mathbb{R}^{128}$ be the $i$-th Legendre polynomial. Next,
we simulate $20$ trajectories based on the same scheme described in section
5.3. In this experiment, since $y_{\text{partial}}(t)$ only depends on
${x}_{1}(t)$ and ${x}_{2}(t)$, we only used these two coordinates to generate
the kernel $K^{x}$ (see (14)). We apply Algorithm 2 with
$|\Omega_{init}|=100$, and obtain an overall median normalized of $10.65\%$
and $0.74\%$ for $\sigma$ and $\rho$, respectively; see box-plots in Fig. 12.
This is evidence that a time-series of partially observed data is sufficient
to estimate the underlying parameters.

##  6\. Discussion

###  6.1. Relevant literature - learning and dynamics

In what follows, we discuss the relation between the main problem of this
paper (parameter estimation with an unknown observation function) to notable
questions at the interface of dynamical systems, statistics, and machine
learning.

This study is related to the fast-growing field of model discovery and machine
learning of physical systems in general [4, 8, 10, 12, 25, 23, 15, 18, 69].
Generally (notwithstanding their many differences), these works aim to
discover governing equations from data using various machine learning
techniques. There are three features which distinguish our study from model
discovery studies. First, our approach and settings are not agnostic to the
physical modelling, and this study does not aim at “physics discovery” [4,
15]. Rather, we use the laws of physics and known models to estimate the
parameters. The second distinction is that the unknown/unspecified portion of
our settings is the correspondence between the model and the observations
(denoted by $G$, see (2)), which in model discovery studies is usually assumed
to be known.

The third distinction between this study and machine learning problems in
general is that we do not have ample training data, i.e., many pairs of a
parameter $\omega$ and the resulting observation $y(t)$. Nor do we even
observe many different signals $y(t)$ which correspond to different (unknown)
parameters $\omega$ (as in e.g., [45]). Critically, even though one can
generate many trajectories $x(t;\omega)$ by solving the underlying ODE, our
data includes only a single experiment with a single observed $y(t)$.

Another increasingly popular application of machine learning to dynamical
systems is learning implicit propagation models. By observing many instances
of a system’s evolution, one learns the time-evolution or Koopman operator to
propagate the observations in time without learning an underlying ODE or PDE,
either in a model-free fashion [27, 51, 68], or using partial knowledge of the
underlying system [2, 5, 67, 70]. In this study, propagating the observations
$y(t)$ (e.g., a video) in time does not seem to advance the estimation of the
system’s parameters $\omega^{*}$.

Our problem can be considered as a novel variant of standard inverse problems
[3, 61]. Broadly speaking, in these problems a known forward (perhaps noisy)
map $\omega\mapsto y(t)=G(x(t;\omega);\zeta)$ is inverted in some sense to
recover $x$ or $\omega^{*}$. This is typically done using the observations $y$
and some a-priori knowledge on the model. The key difference between standard
inverse problems and this study is that $G$ in our case is completely unknown,
and we do not attempt to recover it. A particularly relevant type of inverse
problems is that where $G$ is only known approximately, due to e.g., modelling
errors or numerical approximation [20, 57].

Also related to this study is the general problem of nonlinear dimensionality
reduction and manifold learning. Since $y=G(x)$ (suppressing noise), then the
$\\{y(t;\omega)~{}|~{}t\geq 0,~{}\omega\in\Omega\\}$ can be viewed as a $m+1$
dimensional manifold in the ambient space $\mathbb{R}^{D}$. Manifold learning
techniques can therefore be applied to recover this low-dimensional structure
[7, 21, 24, 25, 60, 73]. These techniques, however, do not provide a
straightforward way to compare the sub-manifold resulting from $y(t)$ with the
model coordinates $x(t;\omega)$. Even diffusion based techniques which allow
one to identify the modality of $x(t;\omega^{*})$ in the observation space
[44, 42] do not lead directly to ways to identify $\omega^{*}$ from the full
parameter space $\Omega$. We consider diffusion-based solutions of our problem
an interesting direction of future studies.

###  6.2. Future works

In this work, we presented two examples (the double pendulum and the Lorenz
’63 system) where the number of estimated parameters is relatively small. As
in other inverse problems, many issues might arise from considering a high-
dimensional parameter, e.g., identifiability, computational efficiency, and
convergence of the optimization scheme. Adjusting and extending our scheme to
high-dimensional parameters is therefore an interesting remaining challenge.

One class of potential applications is the deduction of homogenized constants
from reduced models. Consider, for concreteness, a Hamiltonian describing
optical propagation of laser beams in a waveguides array, or of a quantum-
mechanical electron in two-dimensional graphene. The dynamics of the
corresponding system can be approximately reduced to an effective Dirac
equation, where short-time dynamics and the micro-structure of the lattice are
homogenized to a single constant wave speed. It might be possible, using our
method, to observe the full dynamics, and by indirect comparison to the
reduced model deduce the homogenized velocity and other constants of the
lattice such as topological charges. Even though this problem involves complex
PDEs, our proposal is to estimate a low dimensional parameter. This is only an
example to a broad class of homogenization schemes in optics, radio-frequency
arrays, fluid dynamics, and continuum mechanics. If successful, such an
approach could allow measurement from complex dynamical systems by indirect
and non-explicit comparison to cheaper and simpler reduced models.

Another potential application of our method is predicting hematoma expansion
after intracerebral hemorrhage based on non-contrast computed tomography (CT).
Over the years pathological observations have led to the development of
several parametric models of the propagation of hemorrhage, e.g., [29, 13,
14]. This is another case where low dimensional parameters underlies high-
dimensional data, and requires the ability to estimate the model parameters of
a dynamical system from noisy observations consisting of only a single
trajectory. Here, the parameter estimation problem is to identify the
dynamical regime. One can investigate the ability of our method to estimate
the model parameters from a small number of CT scans from a single patient,
thereby allowing us to determine whether the hematoma is likely to expand,
which in turn necessitates a life-saving, yet risky medical intervention, or
whether the hematoma contracts and only monitoring is required. We believe
that there exists a large number of such applications from the realm of
medical data analysis and related fields, which can benefit from the presented
method.

The approach presented in this paper can be viewed as a general scheme (see
Fig. 2) of which Algorithms 1 and 2 are two effective representatives. Our
choice of Gaussian kernels in (14) and (13) is judicious, standard, and
effective, but might be improved. One avenue for improvement would be
considering kernels that take into account the temporal progression of the
samples, as in e.g., [25, 62]. One of the main constraints in this paper is
that we observe only a single time-series $y(t)$. Suppose the problem is
extended to measure many such time-series, either by altering $\omega^{*}$ or
by changing the initial conditions, then an appropriate kernel might be
learned, see e.g., [41, 22, 50].

##  7\. Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Derivation of (9)

We write here the details for $\bar{\Phi}$ and $K^{y}$.

$\displaystyle\|\bar{\Phi}\|_{F}^{2}$ $\displaystyle={\rm
Tr}(\bar{\Phi}^{\top}\bar{\Phi})$
$\displaystyle=\sum\limits_{j=1}^{N}\langle\bar{\Phi}_{\cdot,j},\bar{\Phi}_{\cdot,j}\rangle$
$\displaystyle=\sum\limits_{j=1}^{N}\langle\phi(y(t_{j}))-\frac{1}{N}\sum\limits_{i=1}^{N}\phi(y(t_{i})),\phi(y(t_{j}))-\frac{1}{N}\sum\limits_{i^{\prime}=1}^{N}\phi(y(t_{i^{\prime}}))\rangle$
$\displaystyle=\sum\limits_{j=1}^{N}\left[\langle\phi(y(t_{j})),\phi(y(t_{j}))\rangle-\frac{2}{N}\sum\limits_{i}\langle\phi(y(t_{j})),\phi(y(t_{i}))\rangle+\frac{1}{N^{2}}\sum\limits_{i,i^{\prime}}\langle\phi(y(t_{i})),\phi(y(t_{i^{\prime}}))\rangle\right]$
$\displaystyle=\sum\limits_{j=1}^{N}\langle\phi(y(t_{j})),\phi(y(t_{j}))\rangle-\frac{2}{N}\sum\limits_{i,j}\langle\phi(y(t_{j})),\phi(y(t_{i}))\rangle+\frac{1}{N^{2}}\sum\limits_{j,i,i^{\prime}}\langle\phi(y(t_{i})),\phi(y(t_{i^{\prime}}))\rangle$
$\displaystyle=\sum\limits_{j=1}^{N}\langle\phi(y(t_{j})),\phi(y(t_{j}))\rangle-\frac{2}{N}\sum\limits_{i,j}\langle\phi(y(t_{j})),\phi(y(t_{i}))\rangle+\frac{N}{N^{2}}\sum\limits_{i,i^{\prime}}\langle\phi(y(t_{i})),\phi(y(t_{i^{\prime}}))\rangle$
$\displaystyle=\sum\limits_{j=1}^{N}k^{y}(y(t_{j}),y(t_{j}))-\frac{1}{N}\sum\limits_{i,\ell}k^{y}(y(t_{i}),y(t_{\ell}))$
$\displaystyle={\rm Tr}(K^{y}H)\,,\qquad H_{ij}=\delta_{ij}-\frac{1}{N}\,,$

where we have used the fact that $\phi$ is real valued and therefore the inner
product is symmetric.

##  Appendix B Proof of Lemma 1

Since $\alpha\cdot\rm Tr(A)=\rm Tr(\alpha A)$ for every scalar $\alpha$ and
square matrix $A$, then

$s(\omega)={\rm
Tr}\left(\frac{K^{x}(\omega)H}{\|K^{x}(\omega)H\|_{F}}\frac{K^{y}H}{\|K^{y}H\|_{F}}\right)=\langle\frac{K^{x}(\omega)H}{\|K^{x}(\omega)H\|_{F}}\,,\,\frac{K^{y}H}{\|K^{y}H\|_{F}}\rangle_{F}\,,$

where, as before, $\langle\cdot,\cdot\rangle_{F}$ is the Frobenius inner
product. Hence we can restrict the analysis to the case of
$\|K_{x}H\|_{F}=\|K_{y}H\|_{F}=1$. We first prove that

$A=\arg\max\limits_{\begin{subarray}{c}B\in M_{n}(\mathbb{R})\\\
\|B\|_{F}=1\end{subarray}}\langle A,B\rangle_{F}\,,$

for any $A\in M_{n}(\mathbb{R})$ with $\|A\|_{F}=1$. Define $\Delta:\,=B-A$.
Then

$1=\|B\|_{F}^{2}=\|A\|_{F}^{2}+\langle
A,\Delta\rangle_{F}+\langle\Delta,A\rangle_{F}+\|\Delta\|_{F}^{2}\,.$

Since $A$ and $B$ are real, $\langle
A,\Delta\rangle_{F}=\langle\Delta,A\rangle_{F}$ and so $\langle
A,\Delta\rangle_{F}=-\|\Delta\|_{F}^{2}/2$. Therefore

$\displaystyle\langle A,B\rangle_{F}$ $\displaystyle=\langle
A,A\rangle_{F}+\langle A,\Delta\rangle_{F}$
$\displaystyle=\|A\|_{F}^{2}-\frac{1}{2}\|\Delta\|_{F}^{2}$
$\displaystyle=1-\frac{1}{2}\|\Delta\|_{F}^{2}\,.$

Therefore, the maximum of this expression is attained when $\|\Delta\|_{F}=0$,
i.e., when $\Delta=0$ and $A=B$.

In the case where $A=K^{y}H$ and $B=K^{x}(\omega)H$, the two centered kernels
are equal if and only if we have the pairwise equalities
$\|y_{i}(\omega^{*})-y_{j}(\omega^{*})\|_{2}=\|x_{i}(\omega)-x_{j}(\omega)\|_{2}$
for all times $t_{i}$ and $t_{j}$. Since $G$ is an
$\ell^{2}(\mathbb{R}^{d}\to\mathbb{R}^{d})$ isometry, this inequality occurs
if and only if
$\|x_{i}(\omega^{*})-x_{j}(\omega^{*})\|_{2}=\|x_{i}(\omega)-x_{j}(\omega)\|_{2}$
for all times $t_{i}$ and $t_{j}$, i.e., where
$x_{i}(\omega)=Tx_{i}(\omega^{*})$ for an $\ell^{2}(\mathbb{R}^{d})$ isometry.

##  Appendix C Explicit ODEs for the double pendulum

For completeness, we include here the Euler-Lagrange ODEs that govern the
double pednulum system; see [56] for derivation and details. Denote by
$\theta_{1}$ and $\theta_{2}$ the angles of the respective pendulums from the
negative $y$ axis, then

$\frac{d}{dt}(\vec{\theta})=\frac{d}{dt}\left(\begin{array}[]{l}\theta_{1}\\\
\theta_{2}\\\ \theta_{3}\\\
\theta_{4}\end{array}\right)=\left(\begin{array}[]{l}\theta_{3}\\\
\theta_{4}\\\ g_{1}(\vec{\theta})\\\ g_{1}(\vec{\theta})\end{array}\right)\,,$

where

$g_{1}(\vec{\theta})=\frac{g(\sin\theta_{2}\cos\Delta\theta-\mu\sin\theta_{1})-(l_{2}\dot{\theta_{2}}^{2}+l_{1}\dot{\theta_{1}}^{2}\cos\Delta\theta)\sin\Delta\theta}{l_{1}(\mu-\cos^{2}\Delta\theta)}\,,$
$g_{2}(\vec{\theta})=\frac{g\mu(\sin\theta_{1}\cos\Delta\theta-\sin\theta_{2})+(\mu
l_{1}\dot{\theta_{1}}^{2}+l_{2}\dot{\theta_{2}}^{2}\cos\Delta\theta)\sin\Delta\theta}{l_{2}(\mu-\cos^{2}\Delta\theta)}\,,$

and where

$\Delta\theta=\theta_{1}-\theta_{2}\,,\qquad\mu=1+\frac{m_{1}}{m_{2}}\,.$

## References

  * [1] Abramovich F and Ritov Y,  Statistical Theory: a Concise Introduction,  CRC Press, 2013.
  * [2] Arcomano T, Szunyogh I, Pathak J, Wikner A, Hunt BR, and Ott E,  A machine learning-based global atmospheric forecast model,  Geo. Res. Lett., 47, e2020GL087776, 2020.
  * [3] Aster RC, Borchers B, and Thurber CH,  Parameter Estimation and Inverse Problems, Elsevier, 2018.
  * [4] Atkinson S, Subber W, Wang L, Khan G, Hawi P, and Ghanem R,  Data-driven discovery of free-form governing differential equations,  arXiv preprint, arXiv:1910.05117, 2019.
  * [5] Azencot O, Erichson NB, Lin V, and Mahoney MW,  Forecasting sequential data using consistent Koopman autoencoders,  arXiv preprint, arXiv:2003.02236, 2020.
  * [6] Bach FT and Jordan MI,  Kernel independent component analysis,  J. Mach. Learn. Res., 3, pp. 1–48, 2002.
  * [7] Belkin M and Niyogi P,  Laplacian eigenmaps for dimensionality reduction and data representation,  Nueral Comput., 15, pp. 1373–1396, 2003.
  * [8] Bertalan T, Dietrich F, Mezic I, and Kevrekidis IG.  On learning Hamiltonian systems from data,  Chaos, 29, 121107, 2020.
  * [9] Blaschko MB and Gretton A. Taxonomy Inference Using Kernel Dependence Measures, Max-Planck-Insitut Technical Report No. 181, 2008.
  * [10] Bongard J and Lipson H,  Automated reverse engineering of nonlinear dynamical systems,  Proc. Natl. Acad. Sci., 104 pp. 9943–9948, 2007.
  * [11] Boots B and Gordon GJ,  Two-manifold problems with applications to nonlinear system identification,  Proc. International Conference on Machine Learning, 2012\.
  * [12] Bramburger JJ, and Kutz JN,  Poincré maps for multiscale physics discovery and nonlinear Floquet theory,  Phys. D, 408, 132479, 2020.
  * [13] Brouwers B, Biffi A, Ayres A M, , Schwab K, Cortellini L, Romero J M, and Goldstein J N.  Apolipoprotein E genotype predicts hematoma expansion in lobar intracerebral hemorrhage.  Stroke 43(6), 1490-1495, 2012
  * [14] Brouwers B, Chang Y, Falcone J, Cai X, Ayres M, Battey T, Vashkevich A, McNamara K, Valant V and Schwab K.  Predicting hematoma expansion after primary intracerebral hemorrhage.  JAMA neurology, pp. 158–164, 2014.
  * [15] Brunton SL, Proctor JL, and Kutz JN.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems  Proc. Natl. Acad. Sci., 113:3932–3937, 2016.
  * [16] Byrd RH, Hribar ME, Nocedal J,  An interior point algorithm for large-scale nonlinear programming,  SIAM J. Optim., 9, 877-900, 1999.
  * [17] Byrd RH, Gilbert JC, and Nocedal J.  A trust region method based on interior point techniques for nonlinear programming,  Math. programming, 89, 149-185, 2000.
  * [18] Champion K, Lusch B, Kutz JN, and Brunton SL.  Data-driven discovery of coordinates and governing equations,  Proc. Natl. Acad. Sci., 116:22445–22451, 2019.
  * [19] Charu A, Hinneburg A, and Keim D,  On the surprising behavior of distance metrics in high dimensional space. _International Conference on Database Theory,_ Springer, Berlin, Germany, 2001.
  * [20] Clearey E, Garbuno-Inigo A, Lan S, Schneider T, and Stuart AM,  Calibrate, emulate, sample,  arXiv preprint, arXiv:200.036989.
  * [21] Coifman RR and Lafon S,  Diffusion maps  Appl. Comput. Harm. Anal., 21, 5–30, 2006.
  * [22] Cortes C, Mehryar M, and Afshin R,  Algorithms for learning kernels based on centered alignment,  The Journal of Machine Learning Research, 2019, 13.1, pp. 795-828.
  * [23] Cranmer M, Sanchez-Gonzalez A, Battaglia P, Xu R, Cranmer K, Spergel D, and Ho S,  Discovering symbolic models from deep learning with inductive biases,  arXiv preprint, arXiv:2006.11287, 2020.
  * [24] Donoho DL and Grimes C,  Hessian eigenmaps: locally linear embedding techniques for high-dimensional data,  Proc. Natl. Acad. Sci., 100, pp. 5591–5596, 2003.
  * [25] Dsilva CJ, Talmon R, Coifman RR, and Kevrekidis IG,  Parsimonious representation of nonlinear dynamical systems through manifold learning: A chemotaxis case study,  Appl. Comput. Harmon. Anal., 44, pp. 759–773, 2018.
  * [26] Dsilva CJ, Talmon R, Gear CW, Coifman RR, and Kevrekidis IG,  Data-driven reduction for a class of multiscale fast-slow stochastic dynamical systems,  SIAM J. Appl. Dyn. Sys., 15, pp. 1327–1351, 2016.
  * [27] Dubois P, Gomez T, Planckaert L, and Perret L.  Data-driven predictions of the Lorenz system,  Phys. D, 132495, 2020.
  * [28] Dunford N and Schwartz JT,  Linear Operators, Part I,  Wiley, New York NY, 1958.
  * [29] Fisher, C Miller.  Pathological observations in hypertensive cerebral hemorrhage.  Journal of Neuropathology & Experimental Neurology pp. 536-550, 1971.
  * [30] Fukumizu K, Bach FR, and Gretton A.  Statistical Consistency of Kernel Canonical Correlation Analysis.  J. Mach. Learn. Res., 8:361–383, 2007.
  * [31] Fukumizu K, Gretton A, Sun X, and Schölkopf B.  Kernel measures of conditional dependence.  In Advances in Neural Information Processing, 2008, pp. 489–496.
  * [32] Gretton A, Bousquet O, Smola A, and Schölkopf B.  Measuring statistical dependence with Hilbert-Schmidt norms.  In International Conference on Algorithmic Learning Theory, 63–77. Springer, Berlin, 2005.
  * [33] Gretton A, Herbrich R, Smola AJ, Bousquet O and Schölkpf B.  Kernel Methods for Measuring Independence,  J. Mach. Learn. Res. 6, 2005, pp. 2075–2129.
  * [34] Gretton A, Fukumizo K, Teo CH, Song K, Schölkopf B, and Smola AJ.  A kernel statistical test of independence,  In Advances in Neural Information Processing, 2008, pp. 585–592.
  * [35] Hinneburg A, Charu A, and Keim D,  What is the nearest neighbor in high dimensional spaces?  26th Internat. Conference on Very Large Databases. 2000\.
  * [36] Hofmann T, Schölkopf B, and Smola AJ.  Kernel Methods in Machine Learning,  Ann. Statist. 36, 2008, pp 1171–1220.
  * [37] Iserles A,  A First Course in the Numerical Analysis of Differential Equations, Cambridge University, UK, 2009.
  * [38] Keller R and Du Q,  Discovery of Dynamics Using Linear Multistep Methods,  arXiv preprint, arXiv:1912.12728, 2020.
  * [39] S. Lafon, Y. Keller, and R. Coifman, “Data fusion and multicue data matching by diffusion maps,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 28 no. 11, p. 1784–1797, 2006.
  * [40] Landau LD and Lifshitz EM,  Course of Theoretical Physics, Vol. 1, Mechanics, Pergamon, Oxford, UK, 1980.
  * [41] Le L, Hao J, Xie Y, and Priestley J.  Deep kernel: Learning kernel function from data using deep neural network,  Proceedings of the 3rd IEEE/ACM International Conference on Big Data Computing, Applications and Technologies. 2016\.
  * [42] Lederman RR and Talmon R,  Learning the geometry of common latent variables using alternating-diffusion,  Appl. Comput. Harm. Anal., 44, pp. 509–536, 2018.
  * [43] Lindenbaum O, Salhov M, Yeredor A, and Averbuch A,  Gaussian bandwidth selection for manifold learning and classification, _Data Mining and Knowledge Discovery_ , 1-37, 2020.
  * [44] Lindenbaum O, Yeredor A, Salhov M, and Averbuch A,  Multi-view diffusion maps,  Inf. Fusion, 55, pp. 127–149, 2020.
  * [45] Liu C-H, Tao Y, Hsu D, Du Q, and Billinge SJL,  Using a machine learning approach to determine the space group of a structure from the atomic pair distribution function  Acta Cryst., A75, PP. 633–643, 2019.
  * [46] Lorenz EN,  Deterministic nonperiodic flow,  J. Atmos. Sci., 20, pp. 130–141, 1963.
  * [47] Michaeli T, Wang W, and Livescu K,  Nonparametric canonical correlation analysis  International Conference on Machine Learning, pp. 1967–1976, 2016.
  * [48] Nguyen VV and Wood EF,  Review and unification of linear identifiability concepts,  SIAM Rev., 24, pp. 34–51, 1982.
  * [49] Nocedal J and Wright SJ,  Numerical Optimization Springer, New York NY, 2006.
  * [50] Owhadi H and Yoo GR,  Kernel flows: from learning kernels from data into the abyss,  J. Comput. Phys., 389, pp. 22–47, 2019.
  * [51] Pathak J, Hunt B, Girvan M, Lu Z, and Ott E,  Model-free prediction of large spatiotemporally chaotic systems from data: a reservoir computing approach,  Phys. Rev. Lett., 120, 024102, 2018.
  * [52] Saitoh S,  Theory of Reproducing Kernels and its Applications, Longman Scientific and Technical, Harlow, UK, 1988.
  * [53] Salhov M, Lindenbaum O, Aizenbud Y, Silberschatz, A, Shkolnisky, Y, and Averbuch, A.  Multi-view kernel consensus for data analysis,  Applied and Computational Harmonic Analysis, 49(1), 208-228.
  * [54] Schölkopf B, Smola A, and Müller K-S,  Nonlinear component analysis as a kernel eigenvalue problem,  Neural Comput., 10, pp. 1299–1319, 1998.
  * [55] Shawe-Taylor J and Cristianini N,  Kernel Methods for Pattern Analysis, Cambridge University, UK, 2004.
  * [56] Shinbrot T, Grebogi C, Wisdom J, and Yorke JA,  Chaos in a double pendulum,  Amer. J. Phys., 60, pp. 491–499. 1992.
  * [57] Shulkind G, Horesh L, and Avron H,  Experimental design for nonparametric correction of misspecified dynamical models,  SIAM/ASA J. Uncertain. Quant., 6, pp. 880–906, 2018.
  * [58] A. Singer, R. Erban, I. Kevrekidis, and R. R. Coifman,  Detecting intrinsic slow variables in stochastic dynamical systems by anisotropic diffusion maps,  Proc. Nat. Acad. Sci., 106, pp. 16090–16095, 2009.
  * [59] Smola AJ and Schölkopf B,  Learning with Kernels, Vol. 4, GMD-Forschungszentrum Informationstechnik, 1998.
  * [60] Sober B and Levin D,  Manifold approximation by moving least-squares projection (MMLS),  Constr. Approx., 2019\.
  * [61] Stuart AM,  Inverse problems: a Bayesian perspective,  Acta Numerica, 19, pp. 451–559, 2010.
  * [62] Talmon R and Coifman RR,  Empirical intrinsic geometry for nonlinear modeling and time series filtering,  Proc. Natl. Acad. Sci., 110, pp. 12535–12540, 2013.
  * [63] Vestner M, Litman R, Rodola A, Bronstein M, Cremers D,  Product manifold filter: Non-rigid shape correspondence via kernel density estimation in the product space,  Proc. Computer Vision and Pattern Recognition, 2017\.
  * [64] Villaverde AF,  Observability and structural identifiability of nonlinear biological systems,  Complexity, 2019\.
  * [65] Waltz RA, Morales JL, Nocedal J, and Orban D,  An interior algorithm for nonlinear optimization that combines line search and trust region steps,  Math. Programming, 107:391–408, 2006.
  * [66] Wang D, Shi L, and Cao J,  Fast algorithm for approximate k-nearest neighbor graph construction.  In In IEEE 13th International Conference on Data Mining Workshops pp. 349-356, 2013.
  * [67] Wikner A, Pathak J, Hunt BR, Girvan M, Arcomano T, Szunyogh I, Pomerance A, and Ott E,  Combining machine learning with knowledge-based modeling for scalable forecasting and subgrid-scale closure of large, complex, spatiotemporal systems,  Chaos, MACL2020, 053111, 2020.
  * [68] Williams MO, Kevrekidis IG, and Rowley CW,  A data-driven approximation of the Koopman operator: extending dynamic mode decomposition,  J. Nonline. Sci., 25, pp. 1307–1346, 2015.
  * [69] Yair O, Talmon R, Coifman RR, and Kevrekidis IG,  Reconstruction of normal forms by learning informed observation geometries from data,  Proc. Natl. Acad. Sci., 114, E7865-E7874, 2017.
  * [70] Yu J and Hesthaven JS,  Flowfield reconstruction method using artificial neural networks,  AIAA J., 57, pp. 482–498, 2019.
  * [71] Zelnik-Manor L and Perona P,  Self-tuning spectral clustering,  in _Advances in Neural Information Processing Systems_ , 2004, pp. 1601–1608.
  * [72] Yair O, Dietrich F, Mulayoff R, Talmon R, Kevrekidis IG,  Spectral Discovery of Jointly Smooth Features for Multimodal Data,  in _arXiv preprint arXiv:2004.04386_ , 2020.
  * [73] Holiday A, Kooshkbaghi M, Bello-Rivas JM, Gear CW, Zagaris A, Kevrekidis IG,  Manifold learning for parameter reduction,  in _Journal of computational physics_ , 2019, 392, pp.419-431.