# Foundational Inference Models for Dynamical Systems

Patrick Seifner  Kostadin Cvejoski  Ramsés J. Sánchez

###### Abstract

Ordinary differential equations (ODEs) underlie dynamical systems which serve
as models for a vast number of natural and social phenomena. Yet inferring the
ODE that best describes a set of noisy observations on one such phenomenon can
be remarkably challenging, and the models available to achieve it tend to be
highly specialized and complex too. In this work we propose a novel supervised
learning framework for zero-shot inference of ODEs from noisy data. We first
generate large datasets of one-dimensional ODEs, by sampling distributions
over the space of initial conditions, and the space of vector fields defining
them. We then learn neural maps between noisy observations on the solutions of
these equations, and their corresponding initial condition and vector fields.
The resulting models, which we call foundational inference models (FIM), can
be (i) copied and matched along the time dimension to increase their
resolution; and (ii) copied and composed to build inference models of any
dimensionality, without the need of any finetuning. We use FIM to model both
ground-truth dynamical systems of different dimensionalities and empirical
time series data in a zero-shot fashion, and outperform state-of-the-art
models which are finetuned to these systems. Our (pretrained) FIMs are
available online.

Machine Learning, ICML

##  1 Introduction

Dynamical systems are mathematical systems that change with time according to
a fixed evolution rule, and serve as representational and analytical tools for
phenomena, both natural and social, which generate patterns that change over
time. Very often, the recorded changes of these empirical patterns are such
that they can be viewed as occurring continuously in time, and thus can be
represented mathematically by systems whose evolution rule is defined through
differential equations. Dynamical systems governed by ordinary differential
equations (ODEs) correspond to an important subset of these models, and
describe the rate of change of a single, generally vector-valued function
(that is, the state of the system) as time evolves, by means of a vector field
(that is, the time derivative of the function). These deceptively simple
systems have had a fundamental role in our understanding of many natural
processes across nearly every scientific discipline — from their very
introduction and application to celestial mechanics in the late seventeenth
century (Newton, 1687; Bernoulli, 1712), to their function as models of
concentration changes in molecular reaction networks (Hoff, 1986);
oscillations in natural populations in biology (Lotka, 1925; Volterra, 1927);
the process of atmospheric convection and its chaotic features (Lorenz, 1963);
and the dynamics of coherent, high energy modes within turbulent flows (Noack
et al., 2003), just to name a few — and continue to be the go-to mathematical
objects for representing dynamic phenomena today.

This work deals with the general problem of finding the ODE (dynamical) system
that best describes a set of noisy observations on a given empirical process.
The current machine learning paradigm tackles this problem by (implicitly)
constraining the model to handle a single process only. Indeed, practitioners
typically encode their inductive biases either into the model architecture or
the training objective, but optimize the model parameters to fit a single
empirical distribution. One problem with this paradigm is that inference
models trained to fit a single process tend to be overly specific to its
distribution, and thus can rarely be reused to infer the governing dynamics of
a second one, even when both processes are assumed to be described by e.g.
similar ODE systems. Another problem is that, to succeed, the paradigm
requires (i) practitioners to have access to enough observations on the
process they study, to train and test their models from scratch; and that (ii)
they also have the experience and expertise to face the trials and
tribulations of their intricate training procedures. Note that in our context,
the inductive bias is nothing but the assumption that an ODE system underlies
the observations, and two representative examples that encode this assumption
into the architecture and the training objective, respectively, are neural
ODEs (Chen et al., 2018) and physics-informed neural networks (Karniadakis et
al., 2021).

A radical departure from this paradigm took place recently — first within the
natural language processing community, with the advent of large language
models, and second with the introduction of neural operators — and seeks
instead to encode (some of) the inductive biases of interest into the training
dataset itself. In fact, large language models have been shown to perform
remarkably well in zero-shot down-stream tasks, after being trained on
datasets which, opposite to the standard single-task paradigm, are built to
contain wildly diverse demonstrations of natural language tasks, in as many
domains and contexts as possible (Radford et al., 2019; Brown et al., 2020).
Similarly, neural operators are able to solve differential equations in zero-
shot settings, even at unseen resolutions, after being trained on carefully
constructed datasets of input functions and their images under the (fixed,
differential) operator being encoded (Lu et al., 2021; Kovachki et al., 2023).

We aim to set this departure onto a more concrete program with the
introduction of Foundational Inference Models (FIM) for dynamical systems. We
define FIMs as pretrained, neural models that perform zero-shot inference of
ODE systems, whose solution encapsulate sets of noisy observations on a given
empirical process. A satisfactory FIM should be able to (I) cope with dynamic
phenomena of very diverse nature, and thus capture various secular and
seasonal patterns, as well as non-periodic fluctuations. It should also be
able to (II) handle observations from empirical processes of different
dimensionalities; (III) process any number of such observations, which might
be recorded regularly in time or not, and adapt its predictions and confidence
accordingly; and (IV) deal with noise signals of different kinds.

The most obvious technical challenges imposed by requirements (II) and (III)
above are that one does not know a priori (that is, at the time of the FIM
construction) either the process dimensionality or how frequently it will be
observed. Similarly, and as regards requirement (IV) above, one does not know
a priori the strength of the to-be-measured noise signal, even after one has
invoked maximum entropy arguments. A naive way of dealing with these
challenges would be to specify ahead a maximum number of dimensions, sequence
length and noise level that the model can process, but it would necessarily
constrain the model’s applicability. We face these challenges by instead
incorporating the notion of compositionality into FIM. Indeed, one of our main
observations is that, without loss of generality, one can treat the components
of the ODE system that fits a given deterministic process as independent, and
hence focus only on the inference of the one-dimensional ODEs describing each
dimension separately. We will empirically show that once such a one-
dimensional FIM is available, it can be used to compose FIMs of any
dimensionality and increase their time-context (that is, their resolution)
without the need of any finetuning.

How can one train a model to infer one-dimensional ODEs from noisy data in a
zero-shot setting? Our key idea is to frame the zero-shot inference problem as
a supervised learning problem, and to encode our inductive biases — viz. ODEs
as dynamical models of 1D (possibly) irregularly-sampled data — directly into
a synthetic training dataset of ODEs and their solutions. Note that this
dataset encodes our beliefs about the class of dynamic phenomena we expect to
encounter in practice, and therefore addresses our pending requirement (I)
from above. We train FIMs to map noisy and sparse observations on the
solutions of these equations to their corresponding ODEs.

In what follows we first briefly review previous work on the inference of ODE
dynamical systems in Section 2, which mostly follow the classical inference
paradigm. Section 3 introduces our main ideas, the data generation process
encoding our assumptions, and the supervised inference algorithm. We test our
methodology on a benchmark of 63 ground-truth dynamical system models, 2 real-
world datasets and one large scale simulation in Section 4. Finally, Section 5
closes the paper with some concluding remarks about future work, while Section
6 comments on the main limitations of our approach.

##  2 Related work

The problem of fitting an ODE dynamical system model to noisy data has many
instances and, consequently, a long history in machine learning. In some
applications, the functional form of the vector field can be constructed from
first principles, and only a few physically interpretable parameters need to
be estimated from data. Earlier solutions to this (instance of the) problem
involved iteratively and numerically solving the corresponding ODE system, for
some initial condition and initial guess of the unknown parameters, and then
updating the latter to maximize the likelihood of the data (Bard, 1974;
Biegler et al., 1986). An alternative approach, commonly known as gradient
matching, avoids solving ODEs altogether and attempts instead to interpolate
the noisy observations, either via splines (Varah, 1982; Ramsay et al., 2007)
or Gaussian processes (Solak et al., 2002; Macdonald et al., 2015), to compute
a finite difference approximation of the vector field. Gradient matching then
proceeds by optimizing the model parameters to minimize the discrepancy
between this finite difference approximation and the functional vector field.

In most applications one deals instead with empirical processes for which no
functional form of the vector field is a priori known; it is this instance of
the problem what we shall focus on. Some early proposals were close in spirit
to gradient matching and simply replaced the unknown vector field with either
linear function models (De Hoon et al., 2002) or Gaussian processes (Äijö &
Lähdesmäki, 2009). More recently, researchers have abandoned data
interpolations and returned to numerically solving the ODE system forward in
time, while leveraging sensitivity methods (Pontryagin, 1987) and modeling the
unknown vector field again with Gaussian processes (Heinonen et al., 2018) or
with neural networks (Chen et al., 2018). Note that these last models encode
their inductive bias into the model architecture whereas (some versions of)
gradient matching encode it into the training objective. Neural ODEs have been
widely employed to model many continuous-time dynamical systems in the last
few years, with examples ranging from direct applications to irregularly-
sampled time series data (Rubanova et al., 2019; Yildiz et al., 2019;
Norcliffe et al., 2021) to parameterization and inference of master equations
(Seifner & Sanchez, 2023). Despite their popularity, the training of neural
ODEs is known to be a formidable hurdle which many have attempted to overcome,
either through regularization (Dupont et al., 2019; Finlay et al., 2020;
Choromanski et al., 2020; Pal et al., 2021) or through the learning of simpler
dynamics (Zhi et al., 2022), but remains a severe limitation of the model.
Perpendicular to these developments runs a large family of models that use
symbolic regression to infer the unknown vector field in symbolic form. We
refer the reader to La Cava et al. (2021) and Makke & Chawla (2024) for
contemporary reviews thereof.

All these inference models find themselves under the umbrella of the classical
paradigm, insofar as they are optimized with respect to a single empirical
process. There are however two recent exceptions. Similar to our approach,
Becker et al. (2023) and d’Ascoli et al. (2023) generate large datasets of ODE
systems and their (noisy) solutions. Opposite to us, they generate their
vector fields as symbolic expressions, which limits expressivity in favor of
interpretability, and tokenise the numerical solution of the corresponding
ODEs. Both works frame the problem of inferring the symbolic vector field from
the ODE solution as a sequence-to-sequence problem, and leverage transformer
models to perform the map. The work of Becker et al. (2023) is however limited
to one dimensional ODEs, whereas that of d’Ascoli et al. (2023) is limited to
six-dimensional ones. Below we shall compare our methodology to that of
d’Ascoli et al. (2023).

Let us finish this section by briefly commenting on the classical alternative
to foundational models, namely transfer learning, which normally involves
freezing a subset of weights in a pretrained model, and finetuning a smaller
subset to solve the new task. Transfer learning has been recently and
successfully applied to ODE inference problems via physically-informed neural
networks (Desai et al., 2022; Pellegrin et al., 2022) and neural operators
(Goswami et al., 2022). Yet, transfer learning inherits the problems of the
classical paradigm, viz. it requires additional datasets and expert knowledge.

##  3 Foundational Inference Models

In this section we introduce a novel methodology to the problem of inferring
ODE systems that model sets of noisy observations on $D$-dimensional empirical
processes, in zero-shot settings.

A system of $D$ coupled ODEs describes the rate of change of a vector-value
function (of time) $\mathbf{x}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{D}$ as
time evolves, by means of a vector field
$\mathbf{f}:\mathbb{R}^{+}\times\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$. Let
us write these equations in the form

$\mathbf{\dot{x}}(t)=\mathbf{f}(t,\mathbf{x}(t)),\,\,\text{where}\,\,\,\mathbf{\dot{x}}(t)=\frac{d\mathbf{x}(t)}{dt}.$
(1)

To model a $D$-dimensional empirical process with an ODE entails assuming that
the available observations correspond to the values taken by the solution of
the ODE, at the observation times, potentially corrupted by some noise signal
of which only a few statistics are known. Now, let us recall that in order to
completely specify a unique ODE solution, it is necessary to specify a set of
auxiliary conditions, like the initial state (condition) of the system.
Accordingly, to infer an ODE system model not only entails inferring the
vector field that defines it, but also the initial condition that specifies
the relevant solution.

Let us now assume that we have found an ODE system defined by a vector field
$\mathbf{f}^{*}$, whose solution $\mathbf{x}^{*}(t)$ with respect to the
initial condition $\mathbf{x}^{*}(0)$ captures the observations made on our
empirical process, up to some additive noise factors. Since the dynamics are
assumed to be deterministic, the process is tied up with that unique solution
$\mathbf{x}^{*}(t)$. Therefore we can think of the ODE system that underlies
our empirical process as defined by the vector field evaluated on the solution
associated to that process, that is $\mathbf{f}^{*}(t,\mathbf{x}^{*}(t))$.
Such a process-dependent vector field is now a function of time only, and its
correspondent (process-dependent) ODE system is now a set of $D$ uncoupled
(that is, independent) albeit complex one-dimensional ODEs.

This simple observation indicates that one does not need to explicitly deal
with systems of coupled equations; inferring one-dimensional ODEs — one per
each dimension of the empirical process in question — suffices. It also
suggests approaching the inference of such 1D ODEs from noisy data in a
different (zero-shot) key. Indeed, if to each (component of our) empirical
process there is an unknown one-dimensional function (of time only) $f^{*}(t)$
that approximately captures its dynamics, one can frame the zero-shot
inference problem as a supervised learning problem, in which one seeks to map
the set of values $x^{*}(\tau_{1}),\dots x^{*}(\tau_{l})$ taken by the ODE
solution $x^{*}(t)$ at some observation times $\tau_{1},\dots,\tau_{l}$ into
both the function $f^{*}(t)$ and the initial condition $x^{*}(0)$ that
specified $x^{*}(t)$ in the first place.

Since we do not know a priori the class of empirical processes we might
encounter in practice — and hence the class of vector fields and initial
conditions that explain them — we opt for (i) choosing probability
distributions over them, which encode our beliefs about the family of
processes we aim to model; and (ii) generating synthetic datasets by sampling
from these distributions and simulating the resulting ODEs. It thus follows
that a model trained to map sets of values of these solutions
$(x(\tau_{1}),\dots x(\tau_{l}))$, on some time grid
$(\tau_{1},\dots,\tau_{l})$, to their corresponding vector field $(f)$ and
initial condition $(x(0))$ will by construction perform zero-shot inference of
1D ODEs on any sequence of empirical observations. This procedure defines our
FIMs for dynamical systems.

The FIM framework consist then of two main components. First, a synthetic
dataset which encodes our beliefs about the empirical phenomena we aim to
model. Second, an inference model that leverages ideas from neural operators
(Lu et al., 2021; Kovachki et al., 2023) to map sets of observations onto
continuous functions. In what follows we shall delve into the details of these
two components.

###  3.1 Synthetic data generation process

In this section we describe the synthetic data generation process, through
which we encode our inductive biases and address the requirements (I), (III)
and (IV) proposed in the introduction. As motivated in the previous section,
the FIM training dataset should consist of a set of one-dimensional ODEs
(defined through their vector field) and their solutions (obtained wrt. a
given initial condition). We define all these functions on the interval
$[0,T]$ for some time horizon $T\in\mathbb{R}^{+}$. The $j$th instance in the
dataset is given by

$\displaystyle f_{j}$ $\displaystyle\sim$
$\displaystyle\phantom{\int}F(\eta_{j}),\quad
f_{j}:[0,T]\rightarrow\mathbb{R},$ (2) $\displaystyle x_{j}(0)$
$\displaystyle\sim$
$\displaystyle\phantom{\int}\mathcal{N}(0,\sigma^{2}_{j}),$ (3) $\displaystyle
x_{j}(t)$ $\displaystyle=$ $\displaystyle
x_{j}(0)+\int_{0}^{t}f_{j}(\tau)d\tau,$ (4)

where $F$ is a distribution over the space of vector fields, whose
hyperparameter $\eta_{j}$ is sampled, together with the variance
$\sigma^{2}_{j}$ over possible initial condition, from a fixed
(meta)distribution $P$. Before giving any details about the distributions $F$
and $P$, we note that there are yet two further assumptions we need to
incorporate into the data generation process, namely our lack of knowledge
about the measurement process (in other words, requirement (III)) and the
noise signal (requirement (IV)). We encode the former into a distribution $G$
over possible observations grids, and the latter into a Gaussian distribution.
To wit

$\displaystyle\tau_{1},\dots,\tau_{L_{j}}$ $\displaystyle\sim$ $\displaystyle
G(\gamma_{j}),\quad\tau_{i}\in(0,T),$ (5)
$\displaystyle\varepsilon_{j1},\dots,\varepsilon_{jL_{j}}$ $\displaystyle\sim$
$\displaystyle\mathcal{N}(0,\rho^{2}_{j}),$ (6) $\displaystyle
y_{j}(\tau_{i})$ $\displaystyle=$ $\displaystyle
x_{j}(\tau_{i})+\varepsilon_{ji},$ (7)

where the hyperparameter $\gamma_{j}$ and the variance $\rho^{2}_{j}$ over
possible noise signals are sampled from the same fixed (meta)distribution $P$.

Distribution over the space of vector fields. What class of vector fields give
rise to ODE solutions that capture observations on naturally occurring
dynamical processes? Or, put another way, what are the most likely vector
fields to be found when approximating empirical processes with ODE dynamical
systems? We do not attempt to answer these questions here. Instead we take a
pragmatic route and define the distribution $F$ over vector fields either via
Gaussian processes, with Radial Basis Function (RBF) kernels (Williams &
Rasmussen, 1995; Vert et al., 2004), or as truncated Chebyshev expansions,
whose coefficients and degree are both randomly sampled. Both alternative have
different hyperparameters ($\eta$ in our context).

On the one hand, RBF kernels depend on a single free parameter $\eta$, which
controls the scale at which variations take place. Thus a smaller $\eta$
results in vector fields with short-range fluctuations, whereas larger $\eta$
produce smoother vector fields that capture broader trends. When using these
kernels we define the distribution $P$ over $\eta$ as a Beta distribution. On
the other hand, truncated Chebyshev expansions can be written as

$f_{j}(t)=\sum_{m=1}^{M_{j}}a^{(j)}_{m}T_{m}(t),$ (8)

where $T_{m}(t)$ is the $m$th Chebyshev polynomial with real coefficient
$a^{(j)}_{m}$. To generate vector fields with Eq. 8 we sample both the degree
$M_{j}$ and the set of coefficients $(a^{(j)}_{1},\dots,a^{(j)}_{M_{j}})$ from
$P$. Intuitions says that we want high order polynomials to occur rarely, and
that the scale of their coefficient should also be small. We thus choose the
distribution $P$ over the degree of the polynomial to follow a Zipf law, and
that over the coefficients to follow a Gaussian distribution whose variance
decreases with the polynomial order like $\frac{1}{m}$.

In practice, we generate a dataset half of which consists of vector fields
sampled from Chebyshev expansions, while the other half consists of vector
fields sampled via Gaussian processes. We refer the reader to Appendix A for
details on hyperparameters and the sizes of our datasets.

Distribution over observation grids. The FIM dataset has an inherent, fine
resolution $\Delta t$ which is specified by the fine grid on $[0,T]$ on which
we evaluate the vector fields and the ODE solutions (that is, Eqs 2 and 4
above). We choose the observation grid to be random because we do not know a
priori whether the empirical processes we aim to study were recorded regularly
or irregularly in time. We also do not know the recording frequency. We thus
consider both regular and irregular observation grid instances. The regular
instances are defined by their instance-dependent resolution
$\Delta\tau^{j}\geq\Delta t$, where $\Delta t$ is the finest resolution within
the dataset. For concreteness, we choose $\Delta\tau^{j}$ uniformly from the
set $(\Delta t,2\Delta t,3\Delta t,4\Delta t)$. The irregular instances are
instead defined by sub-sampling the fine grid with a Bernoulli mask of
probability $p=\frac{1}{2}$ and $\frac{1}{4}$. Both subsampling procedures
generate observation grids with at least $\frac{L}{4}$ points, where $L\Delta
t=T$. This number defines the minimum number of context points needed by FIM
to function. Additional information is provided in Appendix A.

Noise model. We choose to model the noise signal with a Gaussian distribution,
simply because when observing a given physical process, one typically only has
access to the mean square error in those observations. According to the
maximum entropy principle, a Gaussian distributions is the best guess one can
make (actually, the most likely distribution) given the available information
(that is, given those first two moments) (Jaynes, 2003).

###  3.2 Supervised inference models

Given a set of noisy observations $(y_{1},\tau_{1}),\dots,(y_{l},\tau_{l})$ on
the ODE solution $x(t)$ as defined in the data generation process (Eq. 7
above), our goal is to estimate both the vector field $f(t)$ and initial
condition $x(0)$ that specified $x(t)$. In other words, we want to reverse the
data generation process. Note that we used $y_{i}$ to label $y(\tau_{i})$.
Also note that we have at our disposition the set of clean observations
$(x_{1},\tau_{1}),\dots,(x_{l},\tau_{l})$ too. Let us now label the set of
noisy observations with $\mathcal{O}$, the set of clean observations with
$\tilde{\mathcal{O}}$, the space of vector fields with $\mathcal{F}$ and that
of initial conditions with $\mathcal{X}_{0}$. We can then split the inference
procedure into two steps: (i) a denoising step that maps $\mathcal{O}$ onto
$\tilde{\mathcal{O}}$; and (ii) an inference step that maps
$\tilde{\mathcal{O}}$ onto $\mathcal{F}$ and $\mathcal{X}_{0}$.

Let us, for the time being, assume that the denoising process has already been
done. Learning a neural map
$h_{\theta}:\tilde{\mathcal{O}}\rightarrow\mathcal{F}$ can be done in the
spirit of neural operators (Lu et al., 2021; Kovachki et al., 2023). The main
difference is that the input space is not a space of functions; rather it is
the space of all possible observations on the ODE solution $x(t)$.

Let us use $\phi$ and $\psi$ to denote feedforward neural networks (FFN) and
sequence processing networks, respectively, where the latter can be e.g. LSTM
(Hochreiter & Schmidhuber, 1997) or transformer networks (Vaswani et al.,
2017). Let us also denote the networks’ parameters with $\theta$. We define a
function
$h_{\theta}:\tilde{\mathcal{O}}\times\mathbb{R}^{+}\rightarrow\mathbb{R}$ as

$h_{\theta}(t)=\phi_{2}(\phi_{1}(t,\theta),\psi_{1}(x_{1},\tau_{1},\dots,x_{l},\tau_{l},\theta),\theta),$
(9)

where the networks $\phi_{1}$ and $\psi_{1}$ can be interpret as the trunk and
branch nets of DeepONets (Lu et al., 2021), albeit the inputs to the latter
(aka sensors) are not fixed in our case. We then define the mean and variance
of a Gaussian distribution over the values that our estimated vector field can
take as

$\hat{f}(t)=\phi_{3}(h(t),\theta),$ (10)
$\log\text{Var}(\hat{f})(t)=\phi_{4}(h(t),\theta),$ (11)

the variance $\text{Var}(\hat{f})$ can be used to estimate the confidence of
the model when inferring ODE from empirical processes. Furthermore, we model
the initial condition $\hat{x}(0)$ as a Gaussian random variable, whose mean
and variance are given by the composition of FFNs with $\psi_{1}$, that is

$\displaystyle\mu_{0}$ $\displaystyle=$
$\displaystyle\phi_{5}\circ\psi_{1}(x_{1},\tau_{1},\dots,x_{l},\tau_{l},\theta),$
(12) $\displaystyle\sigma^{2}_{0}$ $\displaystyle=$
$\displaystyle\phi_{6}\circ\psi_{1}(x_{1},\tau_{1},\dots,x_{l},\tau_{l},\theta).$
(13)

We can then integrate the mean function $\hat{f}(t)$ to get an estimate of the
ODE solution as follows

$\hat{x}(t)=\hat{x}(0)+\int_{0}^{t}\hat{f}(\tau)d\tau.$ (14)

We refer the reader to Appendix B where we provide an alternative formulation
of our neural operator (Eq. 9) in terms of the kernel framework of Kovachki et
al. (2023).

Training objective. We train to model to maximize the log-likelihood of the
values taken by the vector field on a fine target grid of $L$ points

$\mathcal{L}_{1}=-\underset{f_{1:L}\sim
p_{{\text{data}}}}{\mathbb{E}}\Big{\\{}\sum_{i=1}^{L}\frac{(f_{i}-\hat{f}_{i})^{2}}{2\text{Var}(\hat{f}_{i})}+\frac{1}{2}\log(\text{Var}(\hat{f}_{i}))\Big{\\}},$
(15)

and minimize the one-step-ahead reconstruction error of the integrated
solution

$\mathcal{L}_{2}=\underset{x_{1:L}\sim
p_{{\text{data}}}}{\mathbb{E}}\sum_{i=1}^{L}(x_{i}-x_{i-1}-\hat{f}_{i-1}\Delta
t)^{2}.$ (16)

Note that $\mathcal{L}_{1}$ corresponds to a supervised objective, whereas
$\mathcal{L}_{1}$ is an unsupervised one.

Supervised denoising model. Our supervised learning framework allows us to
treat the inference and denoising problems separately. This separation opens
up the door to the composition of FIMs with different families of denoising
models. In our experiments below we shall use the popular Savitzky–Golay
filter (Savitzky & Golay, 1964). However in Appendix C, we present a
supervised model which can denoise time-series data in zero-shot settings,
too.

Table 1: Reconstruction on ODEBench. We report the percentage ($\%$) of predictions for which the $R^{2}$ score exceeds a threshold of 0.9. Results for ODEFormer, ODEFormer-opt and PySR are extracted from (d’Ascoli et al., 2023). Our results are averaged over 10 noise and subsampling realizations. Model | $\rho=0,\sigma=0$ | $\rho=0,\sigma=0.05$ | $\rho=0.5,\sigma=0$ | $\rho=0.5,\sigma=0.05$
---|---|---|---|---
ODEFormer-opt | 75.9 | 55.7 | 74.7 | 66.6
ODEFormer | 71.2 | 52.2 | 69.9 | 60.3
PySR | 82.3 | 63.2 | 77.0 | 38.2
1-5 FIM |  84.9 $\pm$ 0 |  76.5 $\pm$ 1.2 | 85.3 $\pm$ 1.5 | 80.2 $\pm$ 1.7
FIM-N | 85.7 $\pm$ 0 | 82.1 $\pm$ 0.9 |  82.9 $\pm$ 1.4 |  79.6 $\pm$ 1.5

###  3.3 Compositionality of FIM

Any foundational inference model should generally fulfill requirements (II)
and (III) from the Introduction. Requirement (II) can be achieved with our
formulation by simply leveraging different copies (instances) of our (1D)FIMs
to process each component of a $D$-dimensional process separately. We
empirically demonstrate in the sections below that we can compose high-
dimensional FIMs that outperform state-of-the-art models finetuned to the
target datasets.

However, there are some aspects of requirement (III) that we have not fully
addressed yet. Although our FIMs can (by construction) process observations
recorded both regularly and irregularly in time, they have been trained to
model time series data with at most $L$ points, and at least $\frac{L}{4}$
points (see Section 3.1).

The upper bounds could limit the applicability of FIMs to time series with
more than $L$ points. Nevertheless, we can split such time series into
successive, overlapping windows and use copies (instances) of FIMs to get
local estimates $\hat{x}(t)$. By interpolating these estimates on the window
overlaps, we can combine them to a global estimate for the whole time series.
We refer to Appendix D for a detailed discussion.

Table 2: Reconstruction error on the Lorenz system data. Results are averaged over the $128$ realizations in the test set Model |  MAE$\downarrow$ |  MSE$\downarrow$ |  R2$\uparrow$
---|---|---|---
LatentODE | 1.93 $\pm$ 0.69 |  9.07 $\pm$ 7.64 |  0.88 $\pm$ 0.11
1-4 FIM |  2.04 $\pm$ 0.20 | 7.91 $\pm$ 1.67 | 0.90 $\pm$ 0.02
FIM-N |  2.06 $\pm$ 0.18 |  8.77 $\pm$ 1.58 |  0.88 $\pm$ 0.02

##  4 Experiments

In this section we test our methodology on three different task which involve
datasets of varying complexity, different dimensionalities, and corrupted by
noise signals of very different nature. The first task consists of inferring a
set of 63 ground-truth dynamical systems, four of which exhibit chaotic
behaviour, and comparing against the state-of-the-art symbolic regression
methods of d’Ascoli et al. (2023) and Cranmer (2023), as well as against
latent ODEs (Chen et al., 2018). The second task deals with dynamical systems
evolving in low-dimensional manifolds and the imputation of missing values.
Specifically, we perform imputation of $20\%$ of (i) the data from the
Carnegie Mellon University motion capture database (as pre-processed by
Heinonen et al. (2018)) and (ii) a Navier-Stokes simulation. Finally, we infer
the hidden dynamics of the ETT (Electricity Transformer Temperature) datasets
collected from electricity transformers (Zhou et al., 2021), which is widely
used in forecasting tasks.

We report the results obtained from two pretrained FIMs in zero-shot settings
only. We label these models FIM and FIM-N. They both need a minimum of 32
(context) points to function and differ only in that FIM-N was trained on
noisy input data. Furthermore, we compose both our FIM with Savitzky–Golay
filters.

###  4.1 Inferring ground-truth dynamical systems

In this first experiment we show the ability of FIMs to infer ground-truth ODE
(dynamical) systems. We consider the recently introduced ODEBench dataset
(d’Ascoli et al., 2023) that consists of 63 ODEs (specifically 23 1D, 28 2D,
10 3D and 2 4D equations) and their solutions, evaluated on a 512 grid. Many
of these ODE systems have been used to model physical phenomena of widely
different nature, and examples include the damped harmonic oscillator, the
Lotka-Volterra model, the Lorentz system and the Rössler attractor — the last
two exhibiting chaotic behaviour.

Figure 1:  Inference on the Lorenz attractor time series on a irregular grid
from the test set, split into their individual components. Noisy observations
(black dots) of the system (black line) are interpolated by LatentODE (magenta
dashed line) and FIM-N (blue line).

Given some set of observations $\\{y_{1},\dots,y_{l}\\}$ on the solutions of
these equations, the task consists of inferring the hidden vector field $f(t)$
and comparing the integrated trajectory $\hat{x}(t)$, with the ground-truth
(noiseless and dense) ODE solution $x(t)$. We recreate the settings in
d’Ascoli et al. (2023) and consider two types of data corruption, namely (i)
multiplicative noise ($y_{i}=(1+\epsilon)x_{i}$,
$\epsilon\sim\mathcal{N}(0,\sigma)$), and (ii) the subsampling of a fraction
$\rho$ from the observations. Furthermore, we choose $\rho$ to be either $0$
or $\frac{1}{2}$ and $\sigma$ to be either $0$ or $0.05$. We compare FIM
against ODEFormer (d’Ascoli et al., 2023) and PySR (Cranmer, 2023) both of
which rely on symbolic regression. Note that these models have already been
shown (d’Ascoli et al., 2023) to be superior to many common symbolic
regression methods (like e.g. SINDy), which we then do not include in our
baselines. ODEFormer can be used in zero-shot mode or can be optimized to the
target dataset (ODEFormer-opt). PySR is fitted to each equation separately.
Details about these fitting procedures can be found in (d’Ascoli et al.,
2023). We compose four FIM along the time axis to process the time series in
all the experiments of this section.

Table 1 shows the percentage ($\%$) of predictions for which the $R^{2}$ score
exceeds a threshold of 0.9, averaged over all 63 ODE solutions. Note that our
models clearly outperform all baselines, even the ones that are finetuned to
the data.

We now set to test FIM against a black-box ODE model which does not output a
symbolic expression and hence has higher expressivity. Let us consider the
Latent ODE framework of Chen et al. (2018), and apply it to the chaotic
Lorentz system, as defined in ODEBench (ID 56). We simulate a set of $1280$
time series, which we split into training, validation and test sets following
a ratio of $(8:1:1)$, for a set of initial conditions sampled from a Gaussian
distribution with mean $\mu=(2.3,8.1,12.4)$ (one of the initial conditions in
ODEBench ID 56) and diagonal standard deviation $\sigma=(1,1,1)$. Again, we
corrupt the data by randomly subsampling a fraction $\rho=\frac{1}{2}$ from
it, and perturbing it with a noise signal of $\sigma=0.05$. Table 2 shows our
(reconstruction) results with respect to the mean absolute error (MAE), mean
square error (MSE) and $R^{2}$ score, averaged over all 128 time series in our
test set. Remarkably, both our FIMs are comparable to (and even better than)
Latent ODE. Figure 1 further demonstrates the behaviour of our models when
applied to an instance of the Lorentz system: Our FIMs nicely fit the ground-
truth solutions.

Ablation study. For completeness we report in Appendix H the reconstruction
errors of FIMs when applied to the ODE systems of dimensions 1, 2, 3 and 4
separately. We also report the mean absolute error (MAE), mean square error
(MSE) and $R^{2}$ scores (that is, not only the percentages). Specifically,
Table 3 shows our results when no denoising model is applied to the input
data. These results should be compared with those in Table 4 which correspond
to FIM + Savitzky–Golay configurations. The main observation we can make is
that combining FIM with a denoising model significantly improves its
performance in all cases (that is, in all systems, regardless of their
dimension). Another interesting observations is that, besides the 3D cases
which contain three out of the four chaotic systems in ODEBench, FIMs perform
(in average) remarkably well (wrt. e.g. $R^{2}$ scores) regardless of the
dimensionality of the system. Figure 6 further shows the FIM reconstruction
for all 63 ODE systems.

Figure 2:  First 3 dimensions of the MC trajectory with imputation in data
space. The trajectory is split into an observed (black dots) and an unobserved
set (black triangles). Blue (Dashed magenta) lines correspond to FIM (ground-
truth) reconstructions.  Figure 3:  Zero-shot interpolation of the last
component of every ETT sub-dataset with FIM-N.

###  4.2 Low dimensional manifolds, amputations and FIM

Having demonstrated that FIMs can successfully infer ground-truth ODE systems
in zero-shot settings, and even outperform models (black-box and not)
finetuned on the empirical distributions under study, we would like to
investigate how our methodology handles more complex, high-dimensional and
experimental data. In particular we study the problem of zero-shot imputation.

One interesting application of FIM is to model systems whose dynamics mainly
evolve in low-dimensional manifolds. One such case is the motion capture (MC)
dataset, which consists of a 50D pose measurements from walking humans, but
which is often studied as evolving in a 3D manifold (Wang et al., 2007). We
use a single trajectory with 94 observations from the preprocessed dataset of
Heinonen et al. (2018), and project the 50D process onto 3D space via PCA. We
then remove $20\%$ of the observations from the middle of the trajectory in
latent space and ask FIM to recover them. FIM can perform imputation by
upsampling (zero-shot super-resolution) the regions outside the gap, and then
composing FIMs along time to interpolate the signal across the gap. Figure 2
displays the first three dimensions of the 50D MC trajectory, together with
the reconstruction from the full 3D dynamics and the reconstruction from our
imputation samples. The results show that zero-shot imputation is indeed
plausible. See also Figure 5 where we additionally show the first 9 components
of the 50D process. As a second experiment we also study imputations in low-
dimensional manifolds of incompressible fluids past a circular cylinder
(Navier-Stokes simulations) in Appendix F.

###  4.3 FIM in the wild

All experiments so far have dealt with physical systems which naturally admit
a description in terms of ODEs. In this last section we consider instead noisy
records of electricity consumption and show that FIM can seamlessly
interpolate the data. Indeed, Figure 3 shows the last of the seven dimensions
of the ETT time series dataset, together with the FIM zero-shot
interpolations. Note how FIM correctly captures the fine structures and
patterns in the process. Details about the ETT dataset can be found in
Appendix G, and a complete study wrt. all seven dimensions of this dataset is
shown in Figure 8.

##  5 Conclusions

In this work we introduced a novel methodology for zero-shot inference of
dynamical systems governed by systems of ordinary differencial equations. We
empirically demonstrated that our foundational models successfully infer
ground-truth dynamical systems of any dimensionality in zero-shot settings. We
also showed that our models can perform zero-shot data imputation of complex,
high-dimensional processes and seamlessly interpolate complex time series
data.

##  6 Limitations

The main limitations of FIM are imposed by our synthetic dataset. Specifically
by the class of vector fields and observation grids we include in it. On the
one hand, inference of high frequency ODE systems is difficult for a single
FIM, because of the lack of long, high frequency vector fields in our current
dataset. On the other hand, FIMs require at least about 32 context points to
function (already a sizeable number) regardless of the complexity of the
vector field. Finding a balance between high and low frequency functions, and
large and small context points within the synthetic data needs further
research.

##  7 Impact statements

Our proposal could have a strong impact in the scientific community.
Specifically from the point of view of accessibility to the more modern
machine learning learning techniques. Indeed, foundation models do not require
machine expert knowledge, which necessarily entails that they can be used by a
broader community of scientist and engineers.

## References

  * Äijö & Lähdesmäki (2009) Äijö, T. and Lähdesmäki, H.  Learning gene regulatory networks from gene expression measurements using non-parametric molecular kinetics.  _Bioinformatics_ , 25(22):2937–2944, 2009.
  * Bard (1974) Bard, Y.  _Nonlinear parameter estimation_.  Academic Press, New York, 1974.
  * Becker et al. (2023) Becker, S., Klein, M., Neitz, A., Parascandolo, G., and Kilbertus, N.  Predicting ordinary differential equations with transformers.  In _International Conference on Machine Learning_ , pp. 1978–2002. PMLR, 2023.
  * Bernoulli (1712) Bernoulli, J.  Extrait de la Réponse de M. Bernoulli à M. Herman, datée de Basle le 7. Octobre 1710. (Both in Mémoires de l’Ac. Royale des Sciences, Boudot, Paris, 1710), 1712.
  * Biegler et al. (1986) Biegler, L. T., Damiano, J. J., and Blau, G. E.  Nonlinear parameter estimation: A case study comparison.  _AIChE Journal_ , 32(1):29–45, 1986.  doi: https://doi.org/10.1002/aic.690320105.
  * Brown et al. (2020) Brown, T., Mann, B., Ryder, N., Subbiah, M., Kaplan, J. D., Dhariwal, P., Neelakantan, A., Shyam, P., Sastry, G., Askell, A., et al.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901, 2020.
  * Chen et al. (2018) Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K.  Neural ordinary differential equations.  _Advances in neural information processing systems_ , 31, 2018.
  * Choromanski et al. (2020) Choromanski, K. M., Davis, J. Q., Likhosherstov, V., Song, X., Slotine, J.-J., Varley, J., Lee, H., Weller, A., and Sindhwani, V.  Ode to an ode.  _Advances in Neural Information Processing Systems_ , 33:3338–3350, 2020.
  * Course & Nair (2023) Course, K. and Nair, P. B.  State estimation of a physical system with unknown governing equations.  _Nature_ , 622(7982):261–267, 2023.
  * Cranmer (2023) Cranmer, M.  Interpretable machine learning for science with pysr and symbolicregression. jl.  _arXiv preprint arXiv:2305.01582_ , 2023.
  * d’Ascoli et al. (2023) d’Ascoli, S., Becker, S., Mathis, A., Schwaller, P., and Kilbertus, N.  Odeformer: Symbolic regression of dynamical systems with transformers.  _arXiv preprint arXiv:2310.05573_ , 2023.
  * De Hoon et al. (2002) De Hoon, M. J., Imoto, S., Kobayashi, K., Ogasawara, N., and Miyano, S.  Inferring gene regulatory networks from time-ordered gene expression data of bacillus subtilis using differential equations.  In _Biocomputing 2003_ , pp. 17–28. World Scientific, 2002.
  * Desai et al. (2022) Desai, S., Mattheakis, M., Joy, H., Protopapas, P., and Roberts, S. J.  One-shot transfer learning of physics-informed neural networks.  In _ICML 2022 2nd AI for Science Workshop_ , 2022.
  * Dupont et al. (2019) Dupont, E., Doucet, A., and Teh, Y. W.  Augmented neural odes.  _Advances in neural information processing systems_ , 32, 2019.
  * Finlay et al. (2020) Finlay, C., Jacobsen, J.-H., Nurbekyan, L., and Oberman, A.  How to train your neural ode: the world of jacobian and kinetic regularization.  In _International conference on machine learning_ , pp. 3154–3164. PMLR, 2020.
  * Goswami et al. (2022) Goswami, S., Kontolati, K., Shields, M. D., and Karniadakis, G. E.  Deep transfer operator learning for partial differential equations under conditional shift.  _Nature Machine Intelligence_ , 4(12):1155–1164, 2022.
  * Heinonen et al. (2018) Heinonen, M., Yildiz, C., Mannerström, H., Intosalmi, J., and Lähdesmäki, H.  Learning unknown ode models with gaussian processes.  In _International conference on machine learning_ , pp. 1959–1968. PMLR, 2018.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long short-term memory.  _Neural Computation_ , 9(8):1735–1780, 1997.
  * Hoff (1986) Hoff, J. H. v. t.  _Studies in Chemical Dynamics_.  Amsterdam: F. Mueller and Co. London: Williams and Norgate, 1986.
  * Jaynes (2003) Jaynes, E. T.  _Probability Theory: the Logic of Science_.  Cambridge University Press, 2003.
  * Karniadakis et al. (2021) Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., and Yang, L.  Physics-informed machine learning.  _Nature Reviews Physics_ , 3(6):422–440, 2021.
  * Klambauer et al. (2017) Klambauer, G., Unterthiner, T., Mayr, A., and Hochreiter, S.  Self-normalizing neural networks.  _Advances in neural information processing systems_ , 30, 2017.
  * Kovachki et al. (2023) Kovachki, N., Li, Z., Liu, B., Azizzadenesheli, K., Bhattacharya, K., Stuart, A., and Anandkumar, A.  Neural operator: Learning maps between function spaces with applications to pdes.  _Journal of Machine Learning Research_ , 24(89):1–97, 2023.
  * La Cava et al. (2021) La Cava, W., Orzechowski, P., Burlacu, B., de França, F. O., Virgolin, M., Jin, Y., Kommenda, M., and Moore, J. H.  Contemporary symbolic regression methods and their relative performance.  _arXiv preprint arXiv:2107.14351_ , 2021.
  * Lorenz (1963) Lorenz, E. N.  Deterministic nonperiodic flow.  _Journal of atmospheric sciences_ , 20(2):130–141, 1963.
  * Loshchilov & Hutter (2017) Loshchilov, I. and Hutter, F.  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_ , 2017.
  * Lotka (1925) Lotka, A. J.  _Elements of physical biology_.  Williams & Wilkins, 1925.
  * Lu et al. (2021) Lu, L., Jin, P., Pang, G., Zhang, Z., and Karniadakis, G. E.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  _Nature machine intelligence_ , 3(3):218–229, 2021.
  * Macdonald et al. (2015) Macdonald, B., Higham, C., and Husmeier, D.  Controversy in mechanistic modelling with gaussian processes.  In Bach, F. and Blei, D. (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 1539–1547. PMLR, 2015.
  * Makke & Chawla (2024) Makke, N. and Chawla, S.  Interpretable scientific discovery with symbolic regression: a review.  _Artificial Intelligence Review_ , 57(1):2, 2024.
  * Newton (1687) Newton, I.  _Philosophiae Naturalis Principia Mathematica_.  cf. especially Liber de Motu I Sects. II, III, VII. London, 1687.
  * Noack et al. (2003) Noack, B. R., Afanasiev, K., MORZYŃSKI, M., Tadmor, G., and Thiele, F.  A hierarchy of low-dimensional models for the transient and post-transient cylinder wake.  _Journal of Fluid Mechanics_ , 497:335–363, 2003.
  * Norcliffe et al. (2021) Norcliffe, A., Bodnar, C., Day, B., Moss, J., and Liò, P.  Neural {ode} processes.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=27acGyyI1BY.
  * Pal et al. (2021) Pal, A., Ma, Y., Shah, V., and Rackauckas, C. V.  Opening the blackbox: Accelerating neural differential equations by regularizing internal solver heuristics.  In _International Conference on Machine Learning_ , pp. 8325–8335. PMLR, 2021.
  * Pellegrin et al. (2022) Pellegrin, R., Bullwinkel, B., Mattheakis, M., and Protopapas, P.  Transfer learning with physics-informed neural networks for efficient simulation of branched flows.  _NeurIPS Workshop on Machine Learning and the Physical Sciences_ , 2022.
  * Pontryagin (1987) Pontryagin, L. S.  _Mathematical theory of optimal processes_.  Routledge, 1987.
  * Radford et al. (2019) Radford, A., Wu, J., Child, R., Luan, D., Amodei, D., Sutskever, I., et al.  Language models are unsupervised multitask learners.  _OpenAI blog_ , 1(8):9, 2019.
  * Ramsay et al. (2007) Ramsay, J. O., Hooker, G., Campbell, D., and Cao, J.  Parameter estimation for differential equations: a generalized smoothing approach.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 69(5):741–796, 2007.
  * Rubanova et al. (2019) Rubanova, Y., Chen, T. Q., and Duvenaud, D.  Latent Ordinary Differential Equations for Irregularly-Sampled Time Series.  In _Advances in Neural Information Processing Systems 32_ , pp. 5321–5331, 2019.
  * Savitzky & Golay (1964) Savitzky, A. and Golay, M. J.  Smoothing and differentiation of data by simplified least squares procedures.  _Analytical chemistry_ , 36(8):1627–1639, 1964.
  * Seifner & Sanchez (2023) Seifner, P. and Sanchez, R. J.  Neural Markov jump processes.  In _Proceedings of the 40th International Conference on Machine Learning_ , volume 202. PMLR, 2023.
  * Shukla & Marlin (2020) Shukla, S. N. and Marlin, B.  Multi-time attention networks for irregularly sampled time series.  In _ICML Workshop on the Art of Learning with Missing Values (Artemiss)_ , 2020.
  * Solak et al. (2002) Solak, E., Murray-Smith, R., Leithead, W., Leith, D., and Rasmussen, C.  Derivative observations in gaussian process models of dynamic systems.  _Advances in neural information processing systems_ , 15, 2002.
  * Varah (1982) Varah, J. M.  A spline least squares method for numerical parameter estimation in differential equations.  _SIAM Journal on Scientific and Statistical Computing_ , 3(1):28–46, 1982.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Vert et al. (2004) Vert, J.-P., Tsuda, K., and Sch”olkopf, B.  A primer on kernel methods.  _Kernel methods in computational biology_ , 47:35–70, 2004.
  * Volterra (1927) Volterra, V.  _Variazioni e fluttuazioni del numero d’individui in specie animali conviventi_ , volume 2.  Societá anonima tipografica” Leonardo da Vinci”, 1927.
  * Wang et al. (2007) Wang, J. M., Fleet, D. J., and Hertzmann, A.  Gaussian process dynamical models for human motion.  _IEEE transactions on pattern analysis and machine intelligence_ , 30(2):283–298, 2007.
  * Williams & Rasmussen (1995) Williams, C. and Rasmussen, C.  Gaussian processes for regression.  _Advances in neural information processing systems_ , 8, 1995.
  * Yildiz et al. (2019) Yildiz, C., Heinonen, M., and Lahdesmaki, H.  Ode2vae: Deep generative second order odes with bayesian neural networks.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Zhi et al. (2022) Zhi, W., Lai, T., Ott, L., Bonilla, E. V., and Ramos, F.  Learning efficient and robust ordinary differential equations via invertible neural networks.  In _International Conference on Machine Learning_ , pp. 27060–27074. PMLR, 2022.
  * Zhou et al. (2021) Zhou, H., Zhang, S., Peng, J., Zhang, S., Li, J., Xiong, H., and Zhang, W.  Informer: Beyond efficient transformer for long sequence time-series forecasting.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 35, pp. 11106–11115, 2021.

##  Appendix A Synthetic data generation process: Additional Information

Sampling the space of vector fields. The synthetic dataset that we use for
training our models has in total 8M (800K for the test set) time series each
containing 128 time points. Half of the time series (4M train, 400K test) are
sampled from GP, where the covariance matrix is parameterized by summation of
two RBF kernels ($k_{1}$ and $k_{2}$). To generate a time series $j$, first we
sample from Beta distribution the lengthscale value ($\eta^{1}_{j}$ and
$\eta^{2}_{j}$) for the two kernels. Namely, for each kernel we have two
separate Beta distributions

$\displaystyle\eta^{1}_{j}\sim\mathcal{B}(2,7),$ (17)
$\displaystyle\eta^{2}_{j}\sim\mathcal{B}(5,2).$ (18)

The first kernel returns functions with faster change where as the second
kernel returns smoother functions.

The second half of the dataset is contained of Chebyshev polynomials where for
each time series $j$ we sample the degree $M_{j}$ from Zipf distribution (the
min degree is 1 and the max degree is 100) and the coefficients of the
polynomial $(a^{(j)}_{1},\dots,a^{(j)}_{M_{j}})$ are sampled from Gaussian
distribution.

$\displaystyle M_{j}\sim Zipf(2),$ (19)
$\displaystyle(a^{(j)}_{1},\dots,a^{(j)}_{M_{j}})\sim\mathcal{N}(0,\frac{1}{M_{j}}).$
(20)

Distribution of Observation Grids. Each generated time series is sub-sampled
to create observable data points for model input. We employ three sub-sampling
schemes: full, regular, and irregular. In the full scheme, the entire time
series is observed. For the regular sub-sampling scheme, observations are
obtained by regularly omitted, typically every second, third, or fourth point
from the fine grid. Lastly, the irregular sub-sampling scheme involves sub-
sampling the fine grid using a Bernoulli mask with probabilities $p=0.5$ and
$0.25$.

##  Appendix B Neural operator (kernel) version

Let’s consider an $M$-layer neural operator

$h_{\theta}(t)=\mathcal{P}\int k_{M}(t,\tau,\theta)u_{M}(\tau,\theta)d\tau,$
(21)

where all functions $u_{i}$ for $i>1$ are given by

$u_{i+1}(t,\theta)=\sigma\left(w_{i}\cdot u_{i}(t)+b_{i}(t)+\int
k_{i}(t,\tau,\theta)u_{i}(\tau,\theta)d\tau\right).$ (22)

Now the function $u_{1}$ has to be dealt with differently because we only have
access to the input function at the observation grid
$(x_{1},\tau_{1}),\dots,(x_{l},\tau_{l})$. Thus we define

$u_{1}(t,\theta)=\sigma\left(\sum_{i=1}^{l}k_{1}(t,\tau_{i},\theta)\cdot(x_{i}\cdot
w_{v})\right),$ (23)

with the kernel $k_{1}$ defined as

$k_{1}(t,\tau_{i})=\frac{\exp\big{\\{}(\phi(t)\cdot
w_{q})\cdot(\phi(\tau_{i})\cdot
w_{k})^{T}/D\big{\\}}}{\sum\limits_{n=1}^{l}\exp\big{\\{}(\phi(t)\cdot
w_{q})\cdot(\phi(\tau_{n})\cdot w_{k})^{T}/D\big{\\}}},$ (24)

where the $\phi:\mathbb{R}\rightarrow\mathbb{R}^{d}$ is the time embedding of
Shukla & Marlin (2020), and the weights $w_{q},w_{k}\in\mathbb{R}^{d\times
d_{\text{att}}}$ and $w_{v}\in\mathbb{R}^{1\times d_{\text{att}}}$ are
learnable.

##  Appendix C Supervised denoising model

Consider

$\hat{x}_{1},\dots,\hat{x}_{l},\hat{\sigma}_{x_{1}},\dots\hat{\sigma}_{x_{l}}=\psi_{2}(y_{1},\tau_{1},\dots,y_{l},\tau_{l},\theta).$
(25)

with $\psi_{2}$ a sequence-processing network like e.g. an LSTM or transformer
network.

We can train this simple model to optimize the supervised loss

$\mathcal{L}_{3}=\underset{(x_{1:l},\sigma)\sim
p_{{\text{data}}}}{\mathbb{E}}\Big{\\{}\sum_{i=1}^{l}\frac{(x_{i}-\hat{x}_{i})^{2}}{2\hat{\sigma}_{x_{i}}^{2}}+\frac{1}{2}\log(2\pi\hat{\sigma}_{x_{i}}^{2})\Big{\\}}$
(26)

##  Appendix D Composing FIMs for long time series

Time series in the training set of FIMs are bounded in length by a fixed upper
bound $L$, as described in Section 3.1. Therefore, FIMs can only (reasonably)
process time series which lengths do not exceed $L$.

To satisfy requirement (III), adaptability of FIMs to any number of
observations, we split the time series into successive, overlapping windows
that are still processable by the FIMs. Applying (instances of) FIMs to each
window individually yields local estimates of the solution $\hat{x}(t)$. To
get a global estimate out of the local estimates, one can combine them on the
overlaps in several way. Here we decided on a simple interpolation scheme,
which is described below for two windows.

Consider two successive, overlapping time windows A and B with, associated to
time windows $[t_{0}^{A},t_{1}^{A}]$ and $[t_{0}^{B},t_{1}^{B}]$ respectively.
By the overlapping assumption, we have

$t_{0}^{A}<t_{0}^{B}<t_{1}^{A}<t_{1}^{B}\quad.$

Applying a FIM to each window yields local solutions $\hat{x}_{A}(t)$ and
$\hat{x}_{B}(t)$ in the respective windows. By defining

$\hat{x}(t)=\begin{cases}\hat{x}_{A}(t),&\text{if }t<t_{0}^{B}\\\
\frac{t_{1}^{A}-t}{t_{1}^{A}-t_{0}^{B}}\hat{x}_{A}(t)+\frac{t-t_{0}^{B}}{{t_{1}^{A}-t_{0}^{B}}}\hat{x}_{B}(t),&\text{if
}t_{0}^{B}\leq t\leq t_{1}^{A}\\\ \hat{x}_{B}(t),&\text{if }t>t_{1}^{A}\\\
\end{cases}$

for $t\in[t_{0}^{A},t_{1}^{B}]$, we combine $\hat{x}_{A}(t)$ and
$\hat{x}_{B}(t)$ to a global solution $\hat{x}(t)$ with linear interpolation
on the overlap.

##  Appendix E Experimental setup

We introduce two models, denoted as FIM and FIM-N, both sharing identical
architectures. The key distinction lies in their input data: the FIM model is
trained using clean observations, while the FIM-N model is trained using noisy
observations.

Model Architectures: The model contains one recurrent sub-module $\psi_{1}$
and six feed-forward sub-modules $\phi_{1-6}$. The size of the hidden state of
the model is $\\*h_{\theta}\in\mathbb{R}^{512}$. The recurrent sub-module
$\psi_{1}$ which gives the input to the branch net is bi-directional LSTM with
hidden size of 512. The branch net is feed-forward neural network with four
hidden layers with size 1024. As trunk net $\phi_{1}$, we use the time
embedding from Shukla & Marlin (2020) in combination with feed-forward neural
network with 4x1024 hidden layers. The dimension of the time embedding is 512.
We combine the output of the trunk net and the branch net by using feed-
forward neural network $\phi_{2}$ with 4 hidden layers of size 1024. The
Gaussian distribution over the values that our estimated vector field is
parameterized by two linear layers, one for the mean $\phi_{3}$ and the second
for the variance $\phi_{4}$. Lastly, for inferring the initial condition we
use feed-forward neural network with 4x1024 layers. For all the feed-forward
networks we use SeLU (Klambauer et al., 2017) as activation and dropout rate
of 0.1.

Hyperparameter Tuning: For the tuning of the hyperparameters we used grid
search. We used AdamW (Loshchilov & Hutter, 2017) optimizer with learning rate
of $1e^{-6}$ and weight decay $1e^{-4}$. The batch size is 1024.

Training Procedure: All the models are trained on 4 A100 80Gb GPUs. They were
train approximately 500 epoch. Early stopping was used as stopping criteria.
The models were trained by maximizing the likelihood. The implementation is
done in Jax111https://jax.readthedocs.io/en/latest/index.html. The code and
trained models are provided also in the supplementary materials.

##  Appendix F Imputations and Two-Dimensional Navier-Stokes

As an application of FIMs to high-dimensional data, we consider the simulation
of a two-dimensional, incompressible Navier-Stokes equation from (Course &
Nair, 2023). The equation is simulated on a two-dimensional grid of size
$199\times 1499$ for a total of $596,602$ states.

For the imputation task of Figure 7, we take the test interval of Figure 3
from (Course & Nair, 2023), which consists of $40$ equidistant observations in
their time interval $[256.2,260.2]$. We impute $8$ points in the middle of the
time series. The frame of Figure 7 is recorded at time $258.2$, the center of
the imputation window.

##  Appendix G ETT Dataset – Electricity Transformer Temperature

The ETT222https://github.com/zhouhaoyi/ETDataset from (Zhou et al., 2021)
contains data of two electricity transformers, labeled $1$ and $2$, with two
resolutions, $15$ minutes (m) and $1$ hour intervals (h), for a total of $4$
sub-datasets. There are $7$ features recorded at every observation.

For the interpolation task of Figures 3 and 8, we chose time series out of the
last $20\%$ of each sub-dataset, which is usually used as a test set for
models trained on ETT data. The time series displayed are of length $512$.

##  Appendix H Additional Results

Table 3: ODEBench - Without denoising model. |  |  MAE$\downarrow$ |  MSE$\downarrow$ |  R2$\uparrow$ |  R2_ACC$\uparrow$ |  RMSE$\downarrow$
---|---|---|---|---|---|---
Dim | Model |  |  |  |  |
1 | FIM - N |  0.556 $\pm$ 0.052 |  1.808 $\pm$ 0.41 |  0.612 $\pm$ 0.245 |  0.828 $\pm$ 0.036 |  0.723 $\pm$ 0.059
FIM |  1.098 $\pm$ 0.141 |  6.395 $\pm$ 1.958 |  -2.155 $\pm$ 3.866 |  0.691 $\pm$ 0.041 |  1.384 $\pm$ 0.178
RBF-70-PR-NOISY |  0.381 $\pm$ 0.04 |  0.785 $\pm$ 0.201 |  0.583 $\pm$ 0.378 |  0.889 $\pm$ 0.03 |  0.493 $\pm$ 0.051
RBF-70-PR |  1.861 $\pm$ 0.25 |  23.099 $\pm$ 7.947 |  -3.651 $\pm$ 4.65 |  0.52 $\pm$ 0.033 |  2.323 $\pm$ 0.302
baseline |  0.32 $\pm$ 0.039 |  0.559 $\pm$ 0.167 |  0.672 $\pm$ 0.237 |  0.907 $\pm$ 0.018 |  0.418 $\pm$ 0.05
2 | FIM - N |  0.183 $\pm$ 0.014 |  0.297 $\pm$ 0.141 |  0.921 $\pm$ 0.024 |  0.85 $\pm$ 0.023 |  0.237 $\pm$ 0.017
FIM |  0.234 $\pm$ 0.028 |  0.456 $\pm$ 0.203 |  0.87 $\pm$ 0.046 |  0.845 $\pm$ 0.034 |  0.32 $\pm$ 0.037
3 | FIM - N |  3.774 $\pm$ 0.19 |  76.846 $\pm$ 11.648 |  0.155 $\pm$ 0.105 |  0.185 $\pm$ 0.024 |  4.953 $\pm$ 0.242
FIM |  4.451 $\pm$ 0.578 |  141.019 $\pm$ 55.513 |  -0.171 $\pm$ 0.227 |  0.175 $\pm$ 0.026 |  5.79 $\pm$ 0.843
RBF-70-PR-NOISY |  3.631 $\pm$ 0.221 |  73.977 $\pm$ 10.079 |  0.214 $\pm$ 0.124 |  0.255 $\pm$ 0.037 |  4.817 $\pm$ 0.282
RBF-70-PR |  4.528 $\pm$ 0.214 |  106.606 $\pm$ 8.225 |  -0.465 $\pm$ 0.211 |  0.135 $\pm$ 0.034 |  5.878 $\pm$ 0.233
baseline |  3.859 $\pm$ 0.301 |  88.525 $\pm$ 15.137 |  0.224 $\pm$ 0.067 |  0.28 $\pm$ 0.048 |  5.078 $\pm$ 0.349
4 | FIM - N |  0.025 $\pm$ 0.003 |  0.001 $\pm$ 0.0 |  0.935 $\pm$ 0.025 |  0.875 $\pm$ 0.177 |  0.031 $\pm$ 0.003
FIM |  0.031 $\pm$ 0.005 |  0.002 $\pm$ 0.001 |  0.857 $\pm$ 0.064 |  0.5 $\pm$ 0.167 |  0.04 $\pm$ 0.006
all | FIM - N |  0.884 $\pm$ 0.034 |  12.99 $\pm$ 1.828 |  0.687 $\pm$ 0.08 |  0.737 $\pm$ 0.022 |  1.156 $\pm$ 0.045
FIM |  1.212 $\pm$ 0.114 |  24.921 $\pm$ 8.79 |  -0.4 $\pm$ 1.417 |  0.671 $\pm$ 0.021 |  1.568 $\pm$ 0.165
2-7 |  ODEFormer (opt) $\rho=0.5,\sigma=0.05$ | - | - | - | 66.6 | -
|  ODEFormer $\rho=0.5,\sigma=0.05$ | - | - | - | 60.3 | -
|  PySR $\rho=0.5,\sigma=0.05$ | - | - | - | 38.2 | -
Table 4: ODEBench - With Savgol filter. |  |  MAE$\downarrow$ |  MSE$\downarrow$ |  R2$\uparrow$ |  R2_ACC$\uparrow$ |  RMSE$\downarrow$
---|---|---|---|---|---|---
Dim | Model |  |  |  |  |
1 | FIM - N |  0.259 $\pm$ 0.016 |  0.341 $\pm$ 0.063 |  0.708 $\pm$ 0.088 |  0.933 $\pm$ 0.016 |  0.335 $\pm$ 0.02
FIM |  0.384 $\pm$ 0.034 |  0.784 $\pm$ 0.181 |  0.537 $\pm$ 0.448 |  0.926 $\pm$ 0.021 |  0.495 $\pm$ 0.04
2 | FIM - N |  0.166 $\pm$ 0.018 |  0.316 $\pm$ 0.203 |  0.942 $\pm$ 0.012 |  0.871 $\pm$ 0.026 |  0.208 $\pm$ 0.023
FIM |  0.14 $\pm$ 0.031 |  0.222 $\pm$ 0.225 |  0.946 $\pm$ 0.025 |  0.9 $\pm$ 0.027 |  0.184 $\pm$ 0.041
3 | FIM - N |  3.661 $\pm$ 0.12 |  70.72 $\pm$ 5.195 |  0.234 $\pm$ 0.087 |  0.24 $\pm$ 0.032 |  4.829 $\pm$ 0.167
FIM |  4.085 $\pm$ 0.532 |  111.568 $\pm$ 39.827 |  0.129 $\pm$ 0.139 |  0.215 $\pm$ 0.024 |  5.26 $\pm$ 0.725
4 | FIM - N |  0.02 $\pm$ 0.002 |  0.001 $\pm$ 0.0 |  0.97 $\pm$ 0.011 |  0.95 $\pm$ 0.105 |  0.024 $\pm$ 0.003
FIM |  0.018 $\pm$ 0.003 |  0.001 $\pm$ 0.0 |  0.964 $\pm$ 0.022 |  0.925 $\pm$ 0.121 |  0.021 $\pm$ 0.004
all | FIM - N |  0.75 $\pm$ 0.022 |  11.49 $\pm$ 0.84 |  0.745 $\pm$ 0.034 |  0.796 $\pm$ 0.015 |  0.982 $\pm$ 0.029
FIM |  0.851 $\pm$ 0.082 |  18.094 $\pm$ 6.297 |  0.668 $\pm$ 0.17 |  0.802 $\pm$ 0.017 |  1.098 $\pm$ 0.115
2-7 |  ODEFormer (opt) $\rho=0.5,\sigma=0.05$ | - | - | - | 66.6 | -
|  ODEFormer $\rho=0.5,\sigma=0.05$ | - | - | - | 60.3 | -
|  PySR $\rho=0.5,\sigma=0.05$ | - | - | - | 38.2 | -
Figure 4:  Inference of the lorenz attractor time series from Figure 1 by our model. Results from our model are shown in three dimensions.  Figure 5:  First 9 dimensions of a motion capture trajectory with imputation in data space. Imputation targets (black circles) are not part of the inputs (black dots). Reconstructions from our model (blue line) are based on interpolating the first 3 PCA dimensions of the trajectory. Reconstructions from the first 3 PCA dimensions of the data (magenta dashed line) are provided for comparison. Figure 6: Reconstructions of noisy observations from of all equations in the ODEBench dataset, using the first set of initial conditions.  Table 5: Mean and standard deviation for reconstruction metrics on a test set of the Lorenz attractor of size $128$.  Model |  MAE$\downarrow$ |  MSE$\downarrow$ |  R2$\uparrow$ |  RMSE$\downarrow$
---|---|---|---|---
LatentODE | 1.93 $\pm$ 0.69 |  9.07 $\pm$ 7.64 |  0.88 $\pm$ 0.11 |  2.83 $\pm$ 1.04
1-5 FIM |  2.04 $\pm$ 0.20 | 7.91 $\pm$ 1.67 | 0.90 $\pm$ 0.02 | 2.80 $\pm$ 0.29
FIM - N |  2.06 $\pm$ 0.18 |  8.77 $\pm$ 1.58 |  0.88 $\pm$ 0.02 |  3.00 $\pm$ 0.27
Figure 7:  Imputation task on Navier Stokes trajectory test trajectory of
length $40$ from (Course & Nair, 2023). We drop $20\%$ out of the center of
the trajectory and impute the missing values with FIM-N. Top: observation.
Mid: Reconstruction with the first $38$ PCA components. Bottom: Reconstruction
via imputation of FIM-N on $38$ PCA components.  Figure 8:  Interpolation of
all $7$ dimensions for one time series out of every ETT sub-dataset with
FIM-N. The columns correspond to the sub-datasets, the rows to the individual
dimensions.