Neural Context Flows for Meta-Learning of Dynamical Systems

## Abstract

Abstract Neural Ordinary Differential Equations (NODEs) often struggle to adapt to new dynamic behaviours caused by parameter changes in the underlying system, even when these dynamics are similar to previously observed behaviours. This problem becomes more challenging when the changing parameters are unobserved, meaning their value or influence cannot be directly measured when collecting data. To address this issue, we introduce Neural Context Flow (NCF), a robust and interpretable Meta-Learning framework that includes uncertainty estimation. NCF uses higher-order Taylor expansion to enable contextual self-modulation, allowing context vectors to influence dynamics from other domains while also modulating themselves. After establishing convergence guarantees, we empirically test NCF and compare it to related adaptation methods. Our results show that NCF achieves state-of-the-art Out-of-Distribution performance on 5 out of 6 linear and non-linear benchmark problems. Through extensive experiments, we explore the flexible model architecture of NCF and the encoded representations within the learned context vectors. Our findings highlight the potential implications of NCF for foundational models in the physical sciences, offering a promising approach to improving the adaptability and generalization of NODEs in various scientific applications. Our code is openly available at https://github.com/ddrous/ncflow .

## 1 Introduction

Figure: Figure 1: Predicted angle $\alpha$ for the simple pendulum $\frac{\text{d}^{2}\alpha}{\text{d}t^{2}}+g\sin\alpha=0$. The OFA Neural ODE that disregards context fails to generalize its oscillation frequency $\nicefrac{{\sqrt{g}}}{{2\pi}}$ to unseen environments, due in part to merged inhomogeneous training data. Our work investigates Neural Context Flows and related Meta-Learning methods to overcome this issue.
Refer to caption: x1.png

A prototypical autonomous dynamical system describes the continuous change, through time $t\in\mathbb{R}^{+}$, of a quantity $x(t)\in\mathbb{R}^{d}$. Its dynamics are heavily dependent on its parameters $c\in\mathbb{R}^{d_{c}}$, according to the (ordinary) differential equation

$$ $\displaystyle\frac{\text{d}x}{\text{d}t}(t)=f(x(t),c),$ (1) $$

where $f:\mathbb{R}^{d}\times\mathbb{R}^{d_{c}}\rightarrow\mathbb{R}^{d}$ is the vector field. Learning a dynamical system from data is synonymous with approximating $f$, a task neural networks have been remarkably good at in recent years (Chen et al., ; Kidger, ).

Consider the problem of reconstructing the mechanical motion of an undamped pendulum, with limited data collected from two domains a.k.a. environments: Earth’s surface, and outer space. Unaware of the physical discrepancies between these environments (gravitational constant, ambient temperature, etc.), we could learn one single environment-agnostic vector field using all available data in a strategy termed One-For-All (OFA) Yin et al. (). However, our model would struggle not only to meaningfully fit such inhomogeneous data, but also to generalize when presented with data from a different environment, e.g., the Moon’s surface. This problem is illustrated in Figure 1 when gravity is the parameter that changes. Similarly, learning one vector field per environment would fail to capture the commonalities between the training environments, which hinders fast adaptation to new ones Yin et al. (). In these conditions, it is imperative to *learn what to learn* from all the data while acknowledging each environment’s unique properties.

Neural Ordinary Differential Equations (Neural ODEs) Chen et al. () have emerged as a versatile class of methods for learning ordinary, stochastic, and controlled differential equations Kidger et al. (). They can be viewed as continuous-time limits of Residual Neural Networks He et al. (); Haber and Ruthotto (); Weinan (), and were popularized in the Dynamical Systems community through Universal Differential Equations Rackauckas et al. (). Trained using Differentiable Programming techniques Nzoyem et al. (); Ma et al. (), they have shown extensive applications in areas like chemical engineering Owoyele and Pal (), geosciences Shen et al. (), turbulence forecasting Portwood et al. (), etc., where the parameters are typically known or the environments fixed.

In Scientific Machine Learning Cuomo et al. () the problem of generalization has largely been tackled by injecting domain knowledge. It is commonly understood that adding a term in Equation 1 that captures as much of the dynamics as possible leads to lower evaluation losses Yin et al. (). For such terms to be added, however, it is essential to have knowledge of the parameters that change, which may then either be directly estimated, or predicted by a neural network within the vector field Rackauckas et al. (). We are naturally left to wonder how to efficiently learn a generalizable dynamical system when such physics are absent.

Under constantly changing experimental conditions, two major obstacles to learning the parameter dependence of a vector field arise:

- (P1)
Limited data: models like Neural ODEs are known to be data-hungry Yin et al. (), and limited data from each environment may not be enough to learn a vector field suitable for all environments.
- (P2)
Unobserved parameters: during the data collection process, one might be unfamiliar with the basic physics of the system, or may not know *which* or *how many* parameters might be worth paying attention to.

Solving these two problems would contribute to the efficient generalization of the learned models, an urgent concern given the tendency of the Machine Learning community to prioritize Out-of-Distribution (OoD) *detection* at the expense of *adaptation* to new settings Arjovsky (). OoD adaptation would massively reduce cost and benefit recent efforts towards foundational models in the physical sciences Subramanian et al. (); McCabe et al. ().

This work presents Neural Context Flows (NCFs), a novel approach for multi-environment generalization of dynamical systems. By leveraging the regularity of the vector field with respect to unobserved parameters, NCFs employ a neural network vector field shared between all environments. This vector field is subsequently modulated by environment-specific latent context vectors. Our contributions are threefold:

- (1)
We introduce a robust Meta-Learning methodology for enhancing the generalizability of dynamical systems. This approach effectively addresses problems (P1) and (P2), challenging the prevailing notion that standard Deep Learning techniques are inherently inefficient for Out-of-Distribution (OoD) adaptation Mouli et al. ().
- (2)
We present an interpretable framework for Multi-Task representation learning, incorporating a straightforward method for uncertainty quantification. This work extends the emerging trend of explainable linearly parameterized physical systems Blanke and Lelarge () to non-linear settings, thus broadening its applicability.
- (3)
We provide a curated set of benchmark problems specifically designed for Meta-Learning of dynamical systems. This collection encompasses a diverse range of problems frequently encountered in the physical science literature, all accessible through a unified interface: [https://github.com/ddrous/gen-dynamics](https://github.com/ddrous/gen-dynamics).

## 2 Related Work

Learning data-driven Neural ODEs that generalize across parameters is only a recent endeavor. To the best of the authors’ knowledge, all attempts to solve (P1) and (P2)
have relied on Multi-Task and Meta-Learning to efficiently adapt to new parameter values, thus producing methods with varying levels of interpolation and extrapolation capabilities.

#### Multi-Task Learning.

Multi-Task Learning (MTL) describes a family of techniques where a model is trained to jointly perform multiple tasks. In Scientific Machine Learning, one of the earliest methods to attack this generalization problem is LEADS Yin et al. (). In LEADS, the vector field is decomposed into shared dynamics $f_{\phi}$ and environment-specific $g^{e}_{\psi}$ components

$$ $\displaystyle\frac{\text{d}x^{e}}{\text{d}t}=f_{\phi}(x)+g^{e}_{\psi}(x),$ (2) $$

where the superscript $e$ identifies the environment in which the dynamical system evolves, and $\{\phi,\psi\}$ are learnable neural network weights.

While LEADS excels at interpolation tasks, it performs poorly during extrapolation Kirchmeyer et al. (). Furthermore, it requires retraining a new network $g^{e}_{\psi}$ each time a new environment is encountered, which can be constraining in scenarios where adaptation is frequently required. Before LEADS, other MTL approaches had been proposed outside the context of dynamical systems Caruana (); Rebuffi et al. (, ); Lee and Parish (). Still, they do not address the crucial adaptation to new tasks which is our focus.

#### Meta-Learning.

Another influential body of work looked at Meta-Learning: a framework in which, in addition to the MTL joint training scheme, shared representation is learned for rapid adaptation to unseen tasks with only minimum data Wang et al. (). The seminal MAML Finn et al. () popularized Gradient-Based Meta-learning (GBML) by nesting an inner gradient loop in the typical training process. Since then, several variants aimed at avoiding over-fitting and reducing cost have been proposed, e.g. ANIL Raghu et al. () and CAVIA Zintgraf et al. (). The latter is a contextual learning approach Garnelo et al. () that partitions learnable parameters into some that are optimized on each environment, and others shared across environments, i.e., *meta-trained*.

DyAd Wang et al. () is one of the earliest Meta-Learning approaches to target generalizable dynamical systems. It learns to represent time-invariant features of a trajectory by an *encoder* network, followed by a *forecaster* network to learn shared dynamics across the different environments.
DyAd only performs well under weak supervision, that is when the underlying (observed) parameters are made known to the loss function via penalization.

Arguably the most successful Meta-Learning method for physical systems is CoDA Kirchmeyer et al. (), which assumes that the underlying system is described by parametrized differential equations whose form is shared by all environments. However, these equations differ by the values of the vector field’s weights, which are produced by a (linear) *hypernetwork*. For the environment indexed by $e$, these weights are computed by(^1^11The general CoDA formulation encompasses the GMBL adaptation rule Kirchmeyer et al. (). Furthermore, MTL models can be identified to Equation 3 with $\theta^{c}=\mathbf{0}$.)

$$ $\displaystyle\theta^{e}=\theta^{c}+W\xi^{e},$ (3) $$

where $\theta^{c}$ and $W$ are shared across environments, and $\xi^{e}\in\mathbb{R}^{d_{\xi}}$ is an environment-specific latent context vector (or simply context). While it achieves state-of-the-art performance on many physical systems, the main limitation of CoDA is its hypernetwork approach, which might hinder parallelism and scaling to large root or target networks. In practice, since they involve two networks instead of one, methods based on hypernetworks require more computational resources to train, and exhibit a more complex optimization landscape Chauhan et al. ().

#### Taylor expansion.

Recently, Blanke and Lelarge () modeled the variability of linearly parameterized dynamical systems with an affine function of low-dimensional environment-specific context vectors. They empirically showed that this improved interpretability, generalization abilities, and computation speed. In our work, we similarly explain non-linear systems, thus generalizing existing work with higher-order Taylor expansion. Additionally, we extract additional benefits such as parallelizability and uncertainty quantification.

## 3 Neural Context Flow (NCF)

A training dataset $\operatorname{\mathcal{D}_{\text{tr}}}:=\big{\{}x^{e}_{i}(\cdot)\big{\}}_{i\in
\operatorname{[\![}1,S\operatorname{]\!]}}^{e\in\operatorname{[\![}1,m
\operatorname{]\!]}}$ is defined as a set of trajectories collected
from $m$ related environments, with $S$ trajectories per environment, each of length $N\in\mathbb{N}^{*}$ over a time horizon $T>0$. Given $\operatorname{\mathcal{D}_{\text{tr}}}$, we aim to find the neural network weights $\theta$ that parameterize a vector field $f_{\theta}$, along with several context vectors $\{\xi^{e}\}_{e\in\operatorname{[\![}1,m\operatorname{]\!]}}$ that modulate its behavior such that

$$ $\displaystyle\frac{\text{d}x_{i}^{e}}{\text{d}t}(t)=f_{\theta}(x_{i}^{e}(t), \xi^{e}),\qquad\forall t\in\left[0,T\right],\quad\forall i\in\operatorname{[\! [}1,S\operatorname{]\!]},\quad\forall e\in\operatorname{[\![}1,m\operatorname{ ]\!]}.$ (4) $$

We learn a single vector field for all environments in our training set $\operatorname{\mathcal{D}_{\text{tr}}}$. The same vector field will be reused, unchanged, for future testing and adaptation to environments in $\operatorname{\mathcal{D}_{\text{ad}}}$.

The vector field $f_{\theta}$ is assumed to not only be *continuous*, but also *smooth* in its second argument $\xi$. Exploiting this inductive bias, we “collect” information from other environments by Taylor-expanding $f_{\theta}$ around any other $\{\xi^{j}\}_{j\in\mathrm{P}}$, where $\mathrm{P}\subseteq\operatorname{[\![}1,m\operatorname{]\!]}$ is a context pool(^2^22We note that $\mathrm{P}$ might include $e$ itself, and is reconstituted at each evaluation of loss function (9). Only its size $p$ is constant.) containing $p\geq 1$ distinct environment indices. This gives rise, for fixed $e$ and $i$, to $p$ Neural ODEs

$$ $\displaystyle\begin{dcases}\frac{\text{d}x^{e,j}_{i}}{\text{d}t}(t)=T^{k}_{f_{ \theta}}(x^{e,j}_{i}(t),\xi^{e},\xi^{j}),\\ x^{e,j}_{i}(0)=x^{e}_{i}(0),\end{dcases}\qquad\forall j\in\mathrm{P},$ (5) $$

where $x^{e,j}_{i}(\cdot)\in\mathbb{R}^{d}$, and $\xi^{e},\xi^{j}\in\mathbb{R}^{d_{\xi}}$. $T^{k}_{f_{\theta}}(\cdot,\xi^{e},\xi^{j})$ denotes the $k$-th order Taylor expansion of $f_{\theta}$ at $\xi^{e}$ around $\xi^{j}$.
In particular, $T^{1}_{f_{\theta}}$ can be written as

$$ $\displaystyle T^{1}_{f_{\theta}}(x^{e,j}_{i},\xi^{e},\xi^{j})$ $\displaystyle=f_{\theta}(x^{e,j}_{i},\xi^{j})+\nabla_{\xi}f(x^{e,j}_{i},\xi^{j })(\xi^{e}-\xi^{j})+o(\|\xi^{e}-\xi^{j}\|),$ (6) $$

where $o(\cdot)$ captures negligible residuals. Equation 6 directly consists of a Jacobian-Vector Product (JVP), making its implementation memory-efficient. Since higher-order Taylor expansions of vector-valued functions do not readily display the same property, we provide the following proposition to facilitate the second-order approximation.

###### Proposition 3.1 (Second-order Taylor expansion with JVPs) .

Assume $f:\mathbb{R}^{d}\times\mathbb{R}^{d_{\xi}}\rightarrow\mathbb{R}^{d}$ is $\mathcal{C}^{2}$ wrt its second argument. Let $x\in\mathbb{R}^{d},\xi\in\mathbb{R}^{d_{\xi}}$, and define $g:\bar{\xi}\mapsto\nabla_{\xi}f(x,\bar{\xi})(\xi-\bar{\xi})$. The second-order Taylor expansion of $f$ around any $\tilde{\xi}\in\mathbb{R}^{\xi}$ is then expressed as

$$ $\displaystyle f(x,\xi)$ $\displaystyle=f(x,\tilde{\xi})+\frac{3}{2}g(\tilde{\xi})+\frac{1}{2}\nabla g( \tilde{\xi})(\xi-\tilde{\xi})+o(\|\xi-\tilde{\xi}\|^{2}).$ (7) $$

###### Proof.

We refer the reader to Appendix A.
∎

By setting $\xi^{e}:=\xi$, and $\xi^{j}:=\tilde{\xi}$, Proposition 3.1 yields an expression for $T^{2}_{f_{\theta}}(\cdot,\xi^{e},\xi^{j})$ in terms of JVPs, an implementation of which we detail in Appendix E. During training (as described below) trajectories from all $p$ Neural ODEs are used within the loss function.

Figure: Figure 2: Illustration of the Neural Context Flow (NCF). (a) Given an initial condition $x_{i}^{e}(0)$ for a training trajectory $i$ in the environment $e$, NCF predicts in parallel, several candidate trajectories $\{\hat{x}^{e,j}_{i}\}_{j\in\mathrm{P}}$ that are all compared to the ground truth $x_{i}^{e}$, upon which $\{\theta,\xi^{e}\}$ is updated. (b) Illustration of the 3-networks architecture for $f_{\theta}$ where the state vector and the arbitrary context are projected into the same representational space before they can interact inside a main network (in purple).
Refer to caption: x2.png

This new framework is called Neural Context Flow (NCF) as is the resulting model. The “flow” term refers to the capability of the context from one environment, i.e., $j$, to influence predictions in another environment, i.e., $e$, by means of Taylor expansion. It allows the various contexts to not only modulate the behavior of the vector field, but to equally modulate theirs, since they are forced to remain close for the Taylor approximations to be accurate. This happens while the same contexts are pushed apart by the diversity in the data. This self-modulation process and the beneficial friction it creates is notably absent from other contextual Meta-Learning approaches like CAVIA Zintgraf et al. () and CoDA Kirchmeyer et al. ().

We depict the NCF framework in Figure 2, along with a 3-networks architecture that projects the contexts and the state vectors into the same representational space before they can interact. Concretely, $x_{i}^{e}(t)$ and the context $\xi$ are first processed independently into $\tilde{x}_{i}^{e}(t)$ (by a state network in blue) and $\tilde{\xi}$ (by a context network in green) respectively, before they are concatenated and fed to a main network (in purple) to produce $f_{\theta}(x_{i}^{e}(t),\xi)$. This allows the framework to automatically balance the potentially nonlinear influence of the context with that of the state vector at each evaluation of the vector field.

###### Proposition 3.1 (Second-order Taylor expansion with JVPs) .

###### Proof.

### 3.1 Meta-Training

Starting from the same initial state $\hat{x}^{e,j}_{i}(0):=x^{e}_{i}(0)$, the Neural ODEs (5) are integrated using a differentiable numerical solver Kidger (); Nzoyem et al. (); ; Chen ():

$$ $\displaystyle\hat{x}^{e,j}_{i}(t)=\hat{x}^{e,j}_{i}(0)+\int_{0}^{t}T^{k}_{f_{ \theta}}(\hat{x}^{e,j}_{i}(\tau),\xi^{e},\xi^{j})\,\text{d}\tau,\quad\forall j \in\mathrm{P}.$ (8) $$

The resulting candidate trajectories are evaluated at time steps $\left\{t_{n}\right\}_{n\in\operatorname{[\![}1,N\operatorname{]\!]}}$ such that $t_{1}=0\text{ and }t_{N}=T$. We feed this to a supervised (inner) loss function

$$ $\displaystyle\ell(\theta,\xi^{e},\xi^{j},\hat{x}^{e,j}_{i},x^{e}_{i}):=\frac{1 }{N\times d}\sum_{n=1}^{N}{\|\hat{x}^{e,j}_{i}(t_{n})-x^{e}_{i}(t_{n})\|_{2}^{ 2}}+\frac{\lambda_{1}}{d_{\xi}}\|\xi^{e}\|_{1}+\frac{\lambda_{2}}{d_{\theta}} \|\theta\|_{2}^{2},$ (9) $$

where
$d,d_{\xi},$ and $d_{\theta}$ are the dimensions of the state space, context vectors, and flattened network weights, respectively. $\|\cdot\|_{1}$ and $\|\cdot\|_{2}$ denote the $L^{1}$ and $L^{2}$ norms, with $\lambda_{1}$ and $\lambda_{2}$ their penalty coefficients, respectively. The $L^{1}$ regularization term is particularly vital in promoting sparsity in a possibly high-dimensional context vector. To ensure its independence from environment count, context pool size, and trajectory count per environment, the overall MSE loss function $\mathcal{L}$ is expressed as in Equation 10; after which it is minimized wrt both weights $\theta$ and contexts $\xi^{1:m}:=\{\xi^{e}\}_{e\in\operatorname{[\![}1,m\operatorname{]\!]}}$ via gradient descent, alternating between updates:

$$ $\displaystyle\mathcal{L}(\theta,\xi^{1:m},\operatorname{\mathcal{D}_{\text{tr} }}):=\frac{1}{m\times S\times p}\sum_{e=1}^{m}\sum_{i=1}^{S}\sum_{j=1}^{p}\ell (\theta,\xi^{e},\xi^{j},\hat{x}^{e,j}_{i},x^{e}_{i}).$ (10) $$

The most effective way to train NCFs is via proximal alternating minimization as described in Algorithm 1. Although more computationally demanding compared to ordinary alternating minimization (see Algorithm 3), it is adept at dealing with non-smooth loss terms like the $L^{1}$ norm in Equation 9 Parikh et al. (). Not to mention that ordinary alternating minimization (Algorithm 3) can easily lead to sub-optimal convergence Attouch et al. (); Li et al. (), while its proximal counterpart converges, with random initialization, almost surely to second-order stationary points provided mild assumptions are satisfied (see Theorem 3.2). We provide a short discussion on those assumptions in Appendix A.

The above comparison demands the definition of two variants of Neural Context Flows. NCF-$\mathit{t}_{2}$ uses *second*-order Taylor expansion in (5), and is trained via *proximal* alternating minimization (Algorithm 1). NCF-$\mathit{t}_{1}$ on the other hand, is implemented using a *first*-order Taylor expansion, and trained using *ordinary* alternating minimization (Algorithm 3). Although less robust than NCF-$\mathit{t}_{2}$, NCF-$\mathit{t}_{1}$ is faster and serves as a powerful baseline in our experiments.

###### Theorem 3.2 (Convergence to second-order stationary points) .

Assume that $\mathcal{L}(\cdot,\cdot,\operatorname{\mathcal{D}_{\text{tr}}})$ satisfies the Kurdyka-Lojasiewicz (KL) property, is $L$ bi-smooth, and $\nabla\mathcal{L}(\cdot,\cdot,\operatorname{\mathcal{D}_{\text{tr}}})$ is Lipschitz continuous on any bounded subset of domain $\mathbb{R}^{d_{\theta}}\times\mathbb{R}^{d_{\xi}\times m}$. Under those assumptions, let $(\theta_{0},\xi^{1:m}_{0})$ be a random initialization and $(\theta_{q},\xi^{1:m}_{q})$ be the sequence generated by Algorithm 1. If the sequence $(\theta_{q},\xi^{1:m}_{q})$ is bounded, then it converges to a second-order stationary point of $\mathcal{L}(\cdot,\cdot,\operatorname{\mathcal{D}_{\text{tr}}})$ almost surely.

Figure: Algorithm 1 Proximal Alternating Minimization

Figure: Algorithm 2 Sequential Adaptation of NCF

###### Proof.

The proof of Theorem 3.2 is straightforward by adapting Assumptions 1 and 4, then Theorem 2 of Li et al. ().
∎

###### Theorem 3.2 (Convergence to second-order stationary points) .

###### Proof.

### 3.2 Adaptation (or Meta-Testing)

Few-shot adaptation to a new environment $e^{\prime}\in\operatorname{[\![}a,b\operatorname{]\!]}$ in or out of the meta-training distribution requires relatively less data, i.e., its trajectory count $S^{\prime}\ll S$. Here, the network weights are frozen, and the goal is to find a context $\xi^{e^{\prime}}$ such that

$$ $\displaystyle\frac{\text{d}x_{i}^{e^{\prime}}}{\text{d}t}(t)=f_{\theta}(x_{i}^ {e^{\prime}}(t),\xi^{e^{\prime}}),\qquad\forall i\in\operatorname{[\![}1,S^{ \prime}\operatorname{]\!]}.$ (11) $$

Our adaptation rule, as outlined in Algorithm 2 is extremely fast, converging in seconds for trainings that took hours. In scenarios where we want to adapt to more than one environment, we outline a bulk version in Section A.5 (see Algorithm 4). Although this algorithm is parallelizable and framed in the same way as the meta-training algorithm, it does not allow for “flow” of contextual information, since this causes significant accuracy degradation. Most importantly, disabling the Taylor expansion at this stage limits the size of the context pool to $p=1$ and significantly improves memory efficiency, a crucial requirement when adapting to a large number of environments.

### 3.3 What’s in a context pool P P \mathrm{P} roman_P ?

The content of the context pool $\mathrm{P}$ not only defines the candidate trajectories we get from Equation 8, but also the speed and memory cost of the meta-training process. We outline 3 pool-filling strategies for selecting $p$ neighboring contexts:

- •
Random-All (RA): all $p$ distinct contexts can be selected by randomly drawing their indices from $\operatorname{[\![}1,m\operatorname{]\!]}$. By repeatedly doing so, we maximize variance and reduce overfitting by ensuring the broadest form of self-modulation – since information can always (in the stochastic limit) flow from any environment into $e$.
- •
Nearest-First (NF): only the $p$ closest contexts to $\xi^{e}$ are selected. Each index is drawn from the discrete uniform distribution $\mathcal{U}\{1,m\}$. This effective strategy encourages environments to form clusters while maintaining intra-cluster self-modulation. (Note, however, that if $p=1$, then $\mathrm{P}=\{e\}$ itself, and no self-modulation occurs.) In our experiments, we observe that this strategy is the most balanced with regards to the context pool size, InD and OoD performance.
- •
Smallest-First (SF): the lowest contexts in $L^{1}$ norm are selected first. Since an environment with context close to 0 can be interpreted as an environment-agnostic feature like in Kirchmeyer et al. (); Blanke and Lelarge (), this strategy prioritizes the flow of information from that base or canonical environment to the one of interest $e$.

### 3.4 InD and OoD Testing

We distinguish two forms of testing. For In-Domain (InD) testing, the environments are the same as the ones used during training. In InD testing data, the underlying parameters of the dynamical system we aim to reconstruct are unchanged, and so the meta-learned context vectors are reused. Out-of-Distribution (OoD) testing considers environments encountered during adaptation. The data is from the same dynamical system, but defined by different parameter values (either interpolated or extrapolated).

In all forms of testing, only the main predicted trajectory corresponding to $j=e$ or $j=e^{\prime}$ is used in the MSE and MAPE metrics computation (see Section C.3 for definitions), thus returning to a standard Neural ODE. Other candidate trajectories are aggregated to provide a measure of uncertainty. The initial conditions that define the trajectories are always unseen, although their distribution never changes from meta-training to meta-testing.

## 4 Main Results

In this section, we evaluate the effectiveness of our proposed framework by investigating two main questions. (i) How good are NCFs at resolving (P1) and (P2) for interpolation tasks (Section 4.2)? (ii) How does our framework compare to state-of-the-art Meta-Learning baselines (Section 4.3)?

### 4.1 Experimental Setting

The NCF framework is evaluated on five seminal benchmarks. The Simple Pendulum (SP) models the periodic motion of a point mass suspended from a fixed point. The Lotka-Volterra (LV) system models the dynamics of an ecosystem in which two species interact. Additional ODEs include a simple model for yeast glycolysis: Glycolytic-Oscillator (GO) Daniels and Nemenman (), and the more advanced Sel’kov Model (SM) Strogatz (); Sel’Kov (). Like GO, SM is non-linearly parameterized, but additionally exhibits starkly different behaviors when its key parameter is varied (see LABEL:fig:attractors). Finally, we consider three PDEs, all with periodic boundary conditions, and cast as ODEs via the method of lines: the non-linear oscillatory Brusselator (BT) model for autocatalytic chemical reactions Prigogine and Lefever (), the Gray-Scott (GS) system also for reaction-diffusion in chemical settings Pearson (), and Navier-Stokes (NS) for incompressible fluid flow Stokes et al. ().

For all problems, the parameters and initial states are sampled from distributions representative of real-world problems observed in the scientific community, and the trajectories are generated using a time-adaptive 4th-order Runge-Kutta solver Virtanen et al. (). For LV, GO, GS, and NS, we reproduce the original guidelines set in Kirchmeyer et al. (), while exposing the data for ODE and PDE problems alike via a common interface. Such use of synthetic data is a common practice in this emerging field of generalizable dynamical systems, where the search for unifying benchmarks remains an open problem Massaroli et al. (). This need for shared datasets and APIs has motivated the GEN-DYNAMICS:[https://github.com/ddrous/gen-dynamics](https://github.com/ddrous/gen-dynamics) open-source initiative, our third and final contribution with this paper. Further details, along with the data generation process, are given in Appendix C.

We now highlight a few key practical considerations shared across experiments. We use the 3-networks architecture depicted in Figure 2b to suitably process the state and context variables. The dimension of the context vector, the context pool’s size and filling strategy, and the numerical integration scheme vary across problems. For instance, we set $d_{\xi}=1024$ for SP and LV, $d_{\xi}=202$ for NS, and $d_{\xi}=256$ for all other problems; while $p=2$ for LV and SM, $p=4$ for LV and GO, and $p=3$ for all PDE problems. Other hyperparameters are carefully discussed in Appendix B.

**Table 1: In-Domain (InD) and adaptation (OoD) test MSEs ($\downarrow$) for the LV, GO, SM, BT, GS and NS problems. The best is reported in bold. The best of the two NCF variants is shaded in grey.**
|  | LV ($\times 10^{-5}$) | GO ($\times 10^{-4}$) |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
|  | #Params | InD | OoD | #Params | InD | OoD |
| CAVIA | 305246 | 91.0$\pm$63.6 | 120.1$\pm$28.3 | 130711 | 64.0$\pm$14.1 | 463.4$\pm$84.9 |
| CoDA | 305793 | 1.40$\pm$0.13 | 2.19$\pm$0.78 | 135390 | 5.06$\pm$0.81 | 4.22$\pm$4.21 |
| NCF-$\mathit{t}_{1}$ | 308240 | 6.73$\pm$0.87 | 7.92$\pm$1.04 | 131149 | 40.3$\pm$9.1 | 19.4$\pm$1.24 |
| NCF-$\mathit{t}_{2}$ | 308240 | 1.68$\pm$0.32 | 1.99$\pm$0.31 | 131149 | 3.33$\pm$0.14 | 2.83$\pm$0.23 |

### 4.2 Interpolation Results

This experiment explores the SP problem discussed in Figure 1. During meta-training, we use 25 environments with the gravity $g$ regularly spaced in $[2,24]$. Each of these environments contains only 4 trajectories with the initial conditions $x^{e}_{i}(0)\sim\left(\mathcal{U}(-\frac{\pi}{3},\frac{\pi}{3}),\mathcal{U}(-1
,1)\right)^{T}$. During adaptation, we interpolate to 2 new environments with $g\in\{10.25,14.75\}$.

**Table 2: In-domain and adaptation testing MSEs ($\downarrow$) with OFA and NCF-$\mathit{t}_{1}$ on the SP problem.**
|  | In-Domain | OoD |
| --- | --- | --- |
| OFA ($\times 10^{0}$) | 1.14 $\pm$ 0.002 | 3.31 $\pm$ 0.026 |
| NCF ($\times 10^{-3}$) | 4.82 $\pm$ 0.08 | 0.25 $\pm$ 0.02 |

Results in Table 2 emphasize the merits of Meta-Learning approaches like NCFs in lieu of naive baseline like One-For-All (OFA), where one context-agnostic vector field is learned for all environments indiscriminately.(^3^33We refer the reader to Yin et al. () for a deeper discussion on One-For-All and One-Per-Env (one vector field for each environment).) Additional results examining the learned context vectors are reported in Section C.2.

### 4.3 Extrapolation Results

A Meta-Learning algorithm is only as good as its ability to extrapolated to unseen environments. In this section we tackle several one-shot generalization problems, i.e. $S^{\prime}=1$. We consider the LV, GO, SM, GS, BT, and NS problems, which all involve adaptation environments outside their meta-training distributions. We compare both variants of NCF to two baselines: CAVIA Zintgraf et al. () which is conceptually the closest GBML method to ours, and CoDA-$\ell_{1}$ Kirchmeyer et al. (). Their hyperparameters, laid out in Appendix B, are tuned for a balance of computational efficacy and performance, all the while respecting each baseline’s key assumptions.

Figure: Figure 3: Grid-wise adaptation on the LV problem showing low MAPEs ($\downarrow$).
Refer to caption: x3.png

At similar parameter counts(^4^44This accounts for the total number of learnable parameters in the framework excluding the context vectors, whose size might vary with the method as explained in Appendix B.), Table 1 shows that CAVIA is the least effective for learning all six physical systems, with a tendency to overfit on the few-shots it receives during both meta-training and meta-testing Mishra et al. (). NCF-$\mathit{t}_{2}$ achieves SoTA OoD results on 5 out of 6 problems. While CoDA retains its superiority on GS, we find that it struggles on non-linear problems, most notably on the SM problem whose trajectories go through 3 distinct attractors in a Hopf bifurcation. This shows that CoDA’s ability to automatically select low-rank adaptation subspaces via a linear hypernetwork decoder is limited, particularly for highly nonlinear systems. We note that NCF-$\mathit{t}_{1}$ equally fails to accurately resolve the SM problem, thus illustrating the value of second-order Taylor expansion and proximal alternating minimization lacking in its design but present in NCF-$\mathit{t}_{2}$. This suggests that proximal operators and Taylor-based regularization prevent NCF-$\mathit{t}_{2}$ from learning spurious associations, which is a problem commonly associated with poor OoD generalization Mouli et al. ().

#### Grid-wide adaptation on LV.

We report in Figure 3 how well NCF performs on 625 adaptation environments obtained by varying its parameters $\beta$ and $\delta$ on a 25$\times$25 uniform grid. It shows a consistently low MAPE below 2%, except in the bottom right corner where the MAPE rises to roughly 15%, still a remarkably low value for this problem. We highlight the 9 meta-training environments in yellow, and the 4 environments used for OoD adaptation for Table 1 in indigo. Remarkably, our training environment’s convex-hull where adaptation is particularly low is much larger compared to CoDA’s (Kirchmeyer et al., , Figure 3).

#### Nonlinear adaptation on SM.

We investigate how the models perform across the SM attractors depicted in LABEL:fig:attractors: training environments $e_{1}$ and $e_{2}$ fall into a limit cycle (L1), $e_{3}$ and $e_{4}$ collapse to a stable equilibrium (E), while $e_{5}$ and $e_{6}$ fall into another limit cycle (L2). LABEL:fig:erropergroup presents the adaptation MSE for each environment in the OoD testing dataset. While CoDA equally fails to capture all 3 attractors, we observe that CAVIA and NCF-$\mathit{t}_{2}$ tend to favor E. But unlike CAVIA, NCF-$\mathit{t}_{2}$ succeeds in capturing the other limit cycles as well, as evidenced by the low MSE on all environments.

## 5 Discussion

### 5.1 Benefits of NCFs

Neural Context Flows provide a robust but flexible framework for learning differential equations. They can handle irregularly-sampled sequences like time series, and can easily be extended to general regression tasks Finn et al. (). Some of their other desirable properties are highlighted below.

Figure: Figure 6: Interpretability with NCF.
Refer to caption: x4.png

#### Massively parallelizable.

Equation 10 indicates that NCFs are massively parallelizable along 3 directions. Indeed, evaluations of $\mathcal{L}$ can be vectorized across $m\times p$ environments, and across all $S$ trajectories. This leads to better use of computational resources for meta-training. We provide details on such vectorized NCF implementation in Appendix E, along with a thorough discussion on its scalability in Section A.6.

Figure: Figure 7: Uncertainty estimation with NCF.
Refer to caption: x5.png

#### Intrepretable.

Understanding how a model adapts to new physical settings is invaluable in many scientific scenarios. Moving information from one environment to another via context-affine transformations Blanke and Lelarge () provides a powerful framework for explaining Multi-Task Learning. With contextual self-modulation via Taylor expansion, we generalize this framework while maintaining computational efficiency. In Figure 6 for instance, we leverage Proposition 1 of Blanke and Lelarge () to showcase system identification with NCF-$\mathit{t}_{1}$, where the underlying physical parameters $\beta$ and $\delta$ of the LV problem are recovered up to a linear transform. The corresponding experiment is carefully detailled in Section A.2.

#### Extendable.

The relative simplicity of our formulation makes NCFs easily adaptable to other Neural ODE models. For instance, it is straightforward to augment the vector field as in APHYNITY Yin et al. (), augment the state space as in ANODE Dupont et al. (), or control the vector field as in Massaroli et al. ().

#### Provides uncertainty.

Uncertainty requirements are increasingly critical in Meta-learned dynamical systems . NCFs provide, by their very definition, a measure of the uncertainty in the model. During testing, all candidate trajectories stemming from available contexts (meta-trained, adapted, or both) can be used to ascertain when or where the model is most uncertain in its predictions. In Figure 7 for instance, we visualize the means and standard deviations (scaled 10 folds for visual exposition) of 9 and 21 candidate forecasts for the Lotka-Volterra and Sel’kov Model problems, respectively.

### 5.2 Limitations

#### Theory.

While massively parallelizable, extendable, interpretable, and applicable to countless areas well outside the learning of dynamics models, Neural Context Flows still need additional theoretical analysis to explain their effectiveness. While the provided uncertainty estimation partly compensates for this shortage, we believe that this creates an interesting avenue, along with other limitations outlined below, for future work.

#### Regularity assumptions.

Another limitation faced by Neural Context Flows lies within the assumptions they make. While differentiability of the vector field $f$ wrt $\xi$ is encountered with the majority of dynamical systems in science and engineering, some are bound to be fundamentally discontinuous. NCFs break down in such scenarios and would benefit from the vast body of research in numerical continuation Allgower and Georg ().

#### Hyperparameters.

Finally, NCF introduces several new hyperparameters such as the context pool size $p$, the context size $d_{\xi}$, the proximity parameter $\beta$ in Algorithm 1, and many more. Although we offer insights into their vital role in the Appendix D, we acknowledge that tuning them complicates the training process.

### 5.3 Broader Impact

While the benefits of NCFs are evidenced in Section 5.1, its negative impacts should not be neglected. For instance, malicious deployment of such adaptable models in scenarios they were not designed for could lead to serious adverse outcomes. This said, our data and models are openly available with minimal risk of misuse.

### 5.4 Conclusion and Future Work

This paper introduces Neural Context Flows (NCFs), an innovative framework that enhances model generalization across diverse environments by exploiting the differentiability of the predictor in its latent context vector. The novel application of Taylor expansion in NCFs facilitates vector field modulation for improved adaptation, enhances interpretability, and provides valuable uncertainty estimates for deeper model understanding. Our comprehensive experiments demonstrate the robustness and scalability of the NCF approach, particularly with respect to its most demanding hyperparameters. Future research will explore the limits of NCFs and their adaptation to even more complex scenarios. This work represents a promising step toward developing foundational models that generalize across scientific domains, offering a fresh and versatile approach to conditioning machine learning models.

## Acknowledgments and Disclosure of Funding

This work was supported by UK Research and Innovation grant EP/S022937/1: Interactive Artificial Intelligence, EPSRC program grant EP/R006768/1: Digital twins for improved dynamic design, EPSRC grant EP/X039137/1: The GW4 Isambard Tier-2 service for advanced computer architectures, and EPSRC grant EP/T022205/1: JADE: Joint Academic Data science Endeavour-2.

## Appendix A Algorithms & Proofs

### A.1 Proof of Proposition 3.1

For ease of demonstration, we propose an equivalent formulation of Proposition 3.1 that disregards the first input of $f$ not necessary for its proof. Although we will consistently write terms like $f(x)$, we emphasize that $x$ is meant to stand for the context, not the state variable.

###### Proposition A.1 (Second-order Taylor expansion with JVPs) .

Assume $f:\mathbb{R}^{d_{\xi}}\rightarrow\mathbb{R}^{d}$ is $\mathcal{C}^{2}$. Let $x\in\mathbb{R}^{\xi}$, and define $g:y\mapsto\nabla f(y)(x-y)$. The second-order Taylor expansion of $f$ around any $x_{0}\in\mathbb{R}^{d_{\xi}}$ is then expressed as

$$ $\displaystyle f(x)$ $\displaystyle=f(x_{0})+\frac{3}{2}g(x_{0})+\frac{1}{2}\nabla g(x_{0})(x-x_{0}) +o(\|x-x_{0}\|^{2}).$ $$

###### Proof.

Let $x,x_{0}\in\mathbb{R}^{d_{\xi}}$. The second-order Taylor-expansion of $f$ includes its Hessian that we view as a 3-dimensional tensor, and we contract along its last axis such that

$$ $\displaystyle f(x)=f(x_{0})+\nabla f(x_{0})(x-x_{0})+\frac{1}{2}\left[\nabla^{ 2}f(x_{0})(x-x_{0})\right](x-x_{0})+o(\|x-x_{0}\|^{2}).$ (12) $$

Next we define $g:y\mapsto\nabla f(y)(x-y)$, and we consider a small perturbation $h\in\mathbb{R}^{d_{\xi}}$ to write

$$ $\displaystyle g(y+h)$ $\displaystyle=\nabla f(y+h)(x-y-h)$ $\displaystyle=\left[\nabla f(y)+\nabla^{2}f(y)(h)+o(\|h\|)\right](x-y-h)\qquad \text{\char 43 }\text{Taylor expansion of $\nabla f$}$ $\displaystyle=\nabla f(y)(x-y)-\nabla f(y)h+\left[\nabla^{2}f(y)h\right](x-y)+ o(\|h\|)+\cancel{o(\|h\|^{2})}$ $$

The Hessian is by definition symmetric along its last two axes, which allows us to rewrite the third term as $\left[\nabla^{2}f(y)(x-y)\right]h$. We then have

$$ $\displaystyle g(y+h)=g(y)+\left[\nabla^{2}f(y)(x-y)-\nabla f(y)\right]h+o(\|h \|),$ $$

which indicates that $\nabla g(y):=\nabla^{2}f(y)(x-y)-\nabla f(y)$, from which we derive

$$ $\forall y\in\mathbb{R}^{d_{\xi}},\quad\nabla^{2}f(y)(x-y)=\nabla g(y)+\nabla f (y).$ $$

In particular,

$$ $\displaystyle\nabla^{2}f(x_{0})(x-x_{0})=\nabla g(x_{0})+\nabla f(x_{0}).$ (13) $$

Plugging this into Equation 12, we have

$$ $\displaystyle f(x)$ $\displaystyle=f(x_{0})+\nabla f(x_{0})(x-x_{0})+\frac{1}{2}\left[\nabla g(x_{0 })+\nabla f(x_{0})\right](x-x_{0})+o(\|x-x_{0}\|^{2})$ $\displaystyle=f(x_{0})+\frac{3}{2}\underbrace{\nabla f(x_{0})(x-x_{0})}_{g(x_{ 0})}+\frac{1}{2}\nabla g(x_{0})(x-x_{0})+o(\|x-x_{0}\|^{2}).$ $$

This concludes the proof.
∎

While this expression of the second-order Taylor expansion of a vector-valued function makes its implementation memory-efficient via Automatic Differentiation (AD), it still relies on nested derivatives, which scale exponentially with the order of the Taylor expansion due to avoidable recomputations. For even higher-order Taylor expansions, this scaling is frightfully inefficient. Taylor-Mode AD is a promising avenue to address this issue Bettencourt et al. [].

###### Proposition A.1 (Second-order Taylor expansion with JVPs) .

###### Proof.

### A.2 Identifiability of physical parameters

Our work critically builds on CAMEL Blanke and Lelarge [] which extensively studies a model similar to NCF-$\mathit{t}_{1}$ like a linearly parametrized system. Its Proposition 1 states, informally, that in the limit of vanishing training loss $\mathcal{L}(\cdot,\cdot,\cdot)=0$, the relationship between the learned context vectors and the system parameters is linear and can be estimated using ordinary least square. This forces the model to learn a meaningful representation of the system instead of overfitting the examples from the training tasks. We aim to validate this proposition in our setting by adapting the LV experiment.

For CoDA Kirchmeyer et al. [] and NCF-$\mathit{t}_{1}$, we use the exact same network architecture: a 4-layer MLP with 224 hidden units. This means no context nor state network is used in NCF-$\mathit{t}_{1}$: the context vector of size $d_{\xi}=2$ is directly concatenated to the state vector as done by Zintgraf et al. []. We note that the results obtained using this configuration further indicate the superiority of NCF, even when model comparison centers on their main/root networks (see also Park et al. [] for a similar model comparison based on parameter count)(^5^55As reported in Appendix B, the primary model comparison approach in this work counts all learnable parameters, including hypernetworks if involved. Only the context vectors are exempt from this count.).

Upon meta-training and meta-testing, we set out to recover the underlying parameters of the Lotka-Volterra systems via a linear transformation of the learned context vectors. We fit a linear regression model to the 9 meta-training context vectors, using the true physical parameters as supervision signal. We test on the 4 adaptation contexts. The results, displayed in Figure 6, adequately illustrate interpretability as stated in [Blanke and Lelarge, , Proposition 1]. They show that our trained meta-parameters recover the underlying system parameters up to a linear transform, and thus enable zero-shot (physical parameter-induced) adaptation via inverse regression.

### A.3 Convergence of proximal alternating minimization

For clarity of exposition, Theorem 3.2 expressing the convergence of Algorithm 1 to second-order stationary points is repeated below.

###### Theorem A.2 (Convergence to second-order stationary points) .

Assume that $\mathcal{L}(\cdot,\cdot,\operatorname{\mathcal{D}_{\text{tr}}})$ satisfies the Kurdyka-Lojasiewicz (KL) property, is $L$ bi-smooth, and $\nabla\mathcal{L}(\cdot,\cdot,\operatorname{\mathcal{D}_{\text{tr}}})$ is Lipschitz continuous on any bounded subset of domain $\mathbb{R}^{d_{\theta}}\times\mathbb{R}^{d_{\xi}\times m}$. Under those assumptions, let $(\theta_{0},\xi^{1:m}_{0})$ be a random initialization and $(\theta_{q},\xi^{1:m}_{q})$ be the sequence generated by Algorithm 1. If the sequence $(\theta_{q},\xi^{1:m}_{q})$ is bounded, then it converges to a second-order stationary point of $\mathcal{L}(\cdot,\cdot,\operatorname{\mathcal{D}_{\text{tr}}})$ almost surely.

We use this section to briefly emphasize that the assumptions of KL property Attouch et al. [] and Lipschitz continuity are mild and easily achievable with neural networks. Interestingly, boundedness of the sequence $(\theta_{k},\xi^{1:m}_{k})$ is guaranteed a priory if $\mathcal{L}$ is coercive Li et al. [], a property we encouraged by regularizing $\ell$ wrt weights and contexts as in Equation 9.

###### Theorem A.2 (Convergence to second-order stationary points) .

### A.4 Ordinary alternating minimization

Here we provide an additional procedure (Algorithm 3) for training Neural Context Flows which we reserved for the NCF-$\mathit{t}_{1}$ variant. While stronger assumptions Li et al. [] are needed to establish convergence guarantees like its proximal extension, is it relatively easier to implement, and exposes fewer hyperparameters to tune.

Figure: Algorithm 3 Ordinary Alternating Minimization

Figure: Algorithm 4 Bulk Adaptation of NCF

### A.5 Bulk adaptation algorithm

Neural Context Flows are designed to be fast at adaptation time. However, one might want to adapt to hundreds or thousands of environments, perhaps to identify where the performance degrades on downstream tasks. In such cases, Algorithm 2 could be slow due to its sequential nature. We provide Algorithm 4 that leverages the same parallelism exploited during training (while restricting information flow from one adaptation environment to the next).

Although highly parallelizable, we realize in practice that Algorithm 4 is susceptible to two pitfalls:

- 1.
Memory scarcity: the operating system needs to allocate enough resources to store the data $\operatorname{\mathcal{D}_{\text{ad}}}$, the combined context vectors $\xi^{a:b}\in\mathbb{R}^{d_{\xi}\times(b-a+1)}$ and backpropagate their gradients in order to adapt all environments at once, which might be impossible if the context pool size $p$ is set too high. To avoid this issue, we recommend using $p\ll m$, or ideally $p=1$ if contextual self-modulation is disabled.
- 2.
Slower convergence: the bulk algorithm could take longer to converge to poorer context vectors when the jointly adapted environments are all very far apart. This is because the contextual self-modulation process would be rendered useless by such task discrepancy, and sampling $j=e$ in the context pool $\mathrm{P}$ to return to standard a Neural ODE in Section 3 would be harder. We find in practice that manually setting $j=e$ by adjusting the vector field to *forego the Taylor expansion* works even better than the sequential counterpart (see Figure 3). A more direct way of achieving the same result is to disregard $\mathrm{P}$ and retain $f_{\theta}$ (rather than its Taylor expansion) in Equation 5.

All adaptation results in this paper use the sequential adaptation procedure in Algorithm 2, except for the grid-wise adaptation in Section 4.

### A.6 Scalability of the NCF algorithms

The Neural Context Flow (NCF) framework demonstrates excellent scalability with respect to various hyperparameters, owing to its innovative design and implementation. This scalability is evident in three key aspects: distributed training capabilities, efficient handling of large context vectors, and utilization of first-order optimization techniques.

Primarily, NCF’s training process can be distributed and parallelized across environments and trajectories (see Equation 10 and Algorithm 4), a feature that distinguishes it from baseline methods which are limited to parallelization across trajectories. Furthermore, the framework employs an efficient approach, as outlined in Proposition 3.1, to avoid materializing Jacobians or Hessians wrt potentially large context vectors, thereby significantly enhancing scalability. Additionally, NCF utilizes only first-order gradients wrt model weights $\theta$, in contrast to methods like CAVIA that require second-order information in their bi-level optimization loop.

Scalability can also be evaluated in terms of component size, particularly in meta-learning adaptation rules that incorporate additional contextual parameters beyond shared model weights. These components may include encoders, hypernetworks, and other mechanisms for generating or refining context vectors necessary for task adaptation. The complexity and memory requirements of contextual meta-learning rules, which are directly related to the size of these components, can be quantified by their parameter count among other metrics. In this regard, NCF maintains a constant memory cost $O(1)$, while baseline methods such as CAVIA and CoDA require additional memory to produce better contexts (see [Park et al., , Table 1]).

Despite these advantages, the NCF framework may face challenges related to memory and computational efficiency due to the requirement of solving $p$ Neural ODEs in (3), as opposed to a single one. In scenarios where all training environments are utilized (i.e., $p=m$), this results in a quadratic cost $O(m^{2})$ for Algorithms 1 and 3. However, our ablation studies in Section D.1 demonstrate that competitive performance can be achieved on most problems using as few as $p=2$ neighboring environments. Additional studies in Sections D.2 and D.3 establish the necessity of expressive context vectors Voynov and Babenko [] and validate the efficacy of the 3-networks architecture, respectively.

It is worth noting that limiting these quantities could directly contribute to improved parameter counts and more interpretable models. Moreover, restricting the total number of environments contributing to the loss (10) at each iteration may further enhance efficiency.

## Appendix B Experimental Details

In this section, we share crucial details that went into the experiments conducted in Section 4 and expanded upon in Appendix C. Unless otherwise specified, the hyperparameters in the next 3 paragraphs were applied to produce all figures and tables in this work.

Each method was tuned to provide an excellent balance of accuracy and efficiency. The main motivation behind our choices was to restraint model sizes both to allow efficient training of all baselines, and fair comparison at roughly equal parameter counts. Irreconcilable differences with the baseline adaptation rules make it difficult to perform a systematic comparison, which is why we focused on parameter count as done in ; all the while noting the all models involved had sufficient capacity to approximate the problems at hand. For NCF and CAVIA, we count the number of learnable parameters in the main network being modulated. For CoDA, we count the root network plus the hypernetwork’s learnable parameters as outlined in Equation 3.

Our main workstation for training and adaptation was fitted with an Nvidia GeForce RTX 4080 graphics card which was used for the SP, LV, SM, and NS problems. Additionally, we used an RTX 3090 GPU of the same generation for the GO problem, and an NVIDIA A100 Tensor Core GPU for the GS problem as its CNN-based architecture was more memory-intensive. In addition to the training hardware and deep learning frameworks (JAX for NCFs, and PyTorch for CAVIA and CoDA) that made apple-to-apple comparison challenging, our network architectures varied slightly across methods, as the subsections below will highlight. This is part of the reason we decided to focus on total parameter count as the great equalizer.

All in all, to ensure fair comparison, we made sure each model was sufficiently large to represent the task at hand (by comparing to commonly used hyperparameters in the literature, e.g., in Kirchmeyer et al. []). Since a full iteration means something different for the frameworks under comparison, we simply made sure the models were trained for sufficiently long, i.e. the lowest-performing method for each task was trained for at least as long as the second-best for that task, and its best checkpoint restored at the end.

### B.1 NCFs

#### Model architecture

The MLPs used as our context, state and main networks (as depicted in Figure 2b) typically have a uniform width. When that is not the case, we use the notation $[h_{\text{in}}\rightarrow h_{1}\rightarrow...\rightarrow h_{d-1}\rightarrow h_
{\text{out}}]$ to summarize the number of units in each layer of an MLP of depth $d$ for brevity. We use a similar notation to summarize the number of channels in a Convolutional Neural Network (CNN). For instance, on the (SP) problem, the context network was the MLP $[256\rightarrow 64\rightarrow 64\rightarrow 64]$ (which indicates a context size of $d_{\xi}=256$), the state network was $[2\rightarrow 64\rightarrow 64\rightarrow 64]$; the main network $[128\rightarrow 64\rightarrow 64\rightarrow 64\rightarrow 2]$ ; all with Swish activations Ramachandran et al. []. Bellow, we follow the same convention to detail the networks used for all problems in this paper’s main comparison (see Table 1).

• LV (context) $[1024\rightarrow 256\rightarrow 256\rightarrow 64]$; (state) $[2\rightarrow 64\rightarrow 64\rightarrow 64]$; (main) $[128\rightarrow 64\rightarrow 64\rightarrow 64\rightarrow 2]$ • GO (context) $[256\rightarrow 64\rightarrow 64\rightarrow 122]$; (state) $[7\rightarrow 122\rightarrow 122\rightarrow 122]$; (main) $[244\rightarrow 64\rightarrow 64\rightarrow 64\rightarrow 7]$ • SM (context, state, main) Exact same architectures as with SP above • BT (context) linear layer with 256 input and 128 outputs features followed by an appropriate reshaping then a 2D convolution with 8 channels and circular padding to maintain image size, a kernel of size $3\times 3$, and Swish activation(^6^66Except for the number of channels, these are default convolution settings we use throughout. The circular padding enforces the periodic boundary condition observed in the trajectory data.); (state) appropriate reshaping before 2D convolution with 8 channels; (main) CNN $[16\rightarrow 64\rightarrow 64\rightarrow 64\rightarrow 2]$ • GS Same network architectures as for the previous BT problem, with the exception that the linear layer in the context network outputs 2048 features • NS (context) A single linear layer with 202 input and 1024 output features reshaped and concatenated to the state and the grid coordinates before processing; (state) No state network is used for this problem; (main) 2-dimensional Fourier Neural Operator (FNO) Li et al. [] with 4 spectral convolution layers, each with 8 frequency modes and hidden layers of width 10. Its lifting operator is a $1\times 1$ convolution with 4 input and 10 output channels respectively – functionally similar to a fully connected layer with weight sharing. The projection operator has two such layers, with 10 inputs, 1 output, and 16 hidden channels in-between.

#### Training hyperparameters

As for the propagation of gradients from the loss function to the learnable parameters in the right-hand-side of the neural ODEs (5), we decided to differentiate through our solvers, rather than numerically computing adjoints Chen et al. []. Since our main solver were the less demanding RK4 and Dopri5 solvers, we didn’t incur the heavy memory cost Differentiable Programming methods are typically known for Nzoyem et al. [], Kidger et al. [].

We use the integration schemes: Dopri5 Wanner and Hairer [] for GO, SM, GS, and BT, with relative and absolute tolerances of $10^{-3}$ and $10^{-6}$ respectively; RK4 for LV with $\Delta t=0.1$; Explicit Euler for NS with $\Delta t=1$. Across all NCF variants, we kept the context pool size under 4 to reduce computational workload ($p=2$ for LV and SM, and $p=4$ for LV and GO, and $p=3$ for all PDE problems). We use the Adam optimizers Diederik [] for both model weights and contexts on ODE problems, and Adabelief Zhuang et al. [] for PDE problems. Their initial learning rates were as follows: $3\times 10^{-4}$ for LV and NS, $10^{-3}$ for GO, BT, and GS, $10^{-4}$ for SM. That learning rate is kept constant throughout the various trainings, except for BT, GS, and NS where it is multiplied by a factor (0.1, 0.5, and 0.1 respectively) after a third of the total number of training steps, and again by the same factor at two-thirds completion. The same initial learning rates are used during adaptation, with the number of iterations typically 1500.

As per Equation 10, the NCF implementation is batched across both environments and trajectories for distributed and faster training. This is in contrast to the two baselines, in which predictions are only batched across trajectories. We use full batches to accelerate training. This means that with LV for instance, all 4, 32, or 1 will be used at once depending on the data split we use (meta-trainingmeta-testing, supportquery). For regularization of the loss function Equation 10, we set $\lambda_{1}=10^{-3}$ for all problems; but $\lambda_{2}=0$ for ODE problems and $\lambda_{2}=10^{-3}$ for PDE problems. For NCF-$\mathit{t}_{2}$ we always use a proximal coefficient of $\beta=10$ except for the LV where $\beta=100$.

The simple pendulum SP problem used in this work leveraged NCF-$\mathit{t}_{2}$. The integrator was Dopri5 with default tolerances, a proximal coefficient $\beta=10^{2}$, a constant learning rate of $10^{-4}$, 500 outer and 20 inner iterations(^7^77The inner iterations used in NCF-$\mathit{t}_{2}$ are not to be confused with the inner gradient update steps in CAVIA, since NCF performs alternating rather than bi-level optimization.). For the other problems highlighted in Table 1, the remaining crucial hyperparameters are defined in Tables 3 and 4.

**Table 3: Hyperparameters for NCF-$\mathit{t}_{1}$**
| Hyperparameters | LV | GO | SM | BT | GS | NS |
| --- | --- | --- | --- | --- | --- | --- |
| Context size $d_{\xi}$ | 1024 | 256 | 256 | 256 | 256 | 202 |
| Pool-filling strategy | NF | RA | RA | NF | RA | NF |
| # Epochs | 10000 | 2400 | 24000 | 10000 | 10000 | 5000 |

**Table 4: Hyperparameters for NCF-$\mathit{t}_{2}$**
| Hyperparameters | LV | GO | SM | BT | GS | NS |
| --- | --- | --- | --- | --- | --- | --- |
| Context size $d_{\xi}$ | 1024 | 256 | 256 | 256 | 256 | 202 |
| Pool-filling strategy | NF | NF | RA | NF | RA | NF |
| # Inner iterations | 25 | 10 | 10 | 20 | 20 | 25 |
| # Outer iterations | 250 | 1000 | 1500 | 1000 | 700 | 250 |

### B.2 CoDA

The reference CoDA implementation from Kirchmeyer et al. [] was readily usable for most problems. We adapted its data generation process to incorporate two new benchmarks introduced in this paper (the SM and BT problems).

#### Model architecture

Neural networks without physics priors are employed throughout, all with the Swish activation function. • LV: 4-layers MLP with 224 neurons per hidden layers (width) • GO: 4-layers MLP with width 146 • SM: 4-layers MLP with width 90 • BT: 4-layers ConvNet with 46 hidden convolutional filters and $3\times 3$ kernels. Its outputs is rescaled by $10^{-4}$ to stabilize rollouts • GS: Same ConvNet as BT, but with 106 filters. • NS: 2-dimensional Fourier Neural Operator (FNO) Li et al. [] with 4 spectral convolution layers, each with 8 frequency modes and hidden layers of width 10. Its lifting operator is a single fully connected layer, while its projection operator has two such layers, with 16 hidden neurons.

#### Training hyperparameters

Using TorchDiffEq Chen [], CoDA backpropagates gradients through the numerical integrator. We use the same integration schemes as NCF above. It stabilizes its training by applying exponential Scheduled Sampling Bengio et al. [] with constant $k=0.99$ and initial gain $\epsilon_{0}=0.99$ updated every 10 epochs (except for the NS problem updated every 15 epochs). The Adam optimizer with a constant learning rate of $10^{-4}$ is used. The batch size and the number of epochs is given in the Table 5.

**Table 5: Hyperparameters for CoDA**
| Hyperparameters | LV | GO | SM | BT | GS | NS |
| --- | --- | --- | --- | --- | --- | --- |
| Context size $d_{\xi}$ | 2 | 2 | 2 | 2 | 2 | 2 |
| Minibatch size | 4 | 32 | 4 | 1 | 1 | 16 |
| # Epochs | 40000 | 40000 | 12000 | 10000 | 120000 | 30000 |

In line with CoDA’s fundamental low-rank assumption, we maintain a context size $d_{\xi}=2$ for all problems, since no more than 2 parameters vary for any given problem.

### B.3 CAVIA

The reference implementation of CAVIA-Concat Zintgraf et al. [] required substantial modifications to fit dynamical systems. Importantly, we incorporated the TorchDiffEq Chen [] open-source package and we adjusted other hyperparameters accordingly to match NCF and CoDA on parameter count and other relevant aspects for fair comparison. The performance of CAVIA-Concat is heavily-dependent on the context size, which we adjust depending on the problem.

#### Model architecture

Like CoDA, the Swish activation is used for all neural networks • LV: 4-layers MLP with width 278 • GO: 4-layers MLP with width 168 • SM: 4-layers MLP with width 84 • BT: 7-layers ConvNet with 46 hidden convolutional filters and $3\times 3$ kernels. Its outputs is rescaled by $0.1$ to stabilize rollouts • GS: 4-layers ConvNet with 184 hidden convolutional filters and $3\times 3$ kernels. Its outputs is rescaled by $0.1$ to stabilize rollouts. • NS: Same FNO as with NCF above, including the single-layer context network with 1024 output neurons.

#### Training hyperparameters

We use full batches to accelerate training. We use the Adam optimizer with a learning rate of 0.001 for the meta-update step. The single inner update step has a learning rate of 0.1. The number of iterations is given in Table 6.

**Table 6: Hyperparameters for CAVIA**
| Hyperparameters | LV | GO | SM | BT | GS | NS |
| --- | --- | --- | --- | --- | --- | --- |
| Context size $d_{\xi}$ | 1024 | 256 | 256 | 64 | 1024 | 202 |
| # Iterations | 200 | 100 | 1200 | 500 | 50 | 1200 |

## Appendix C Datasets & Additional Results

### C.1 Gen-Dynamics

Given the lack of benchmark consistency in Scientific Machine Learning Massaroli et al. [], we launched the Gen-Dynamics initiative(^8^88[https://github.com/ddrous/gen-dynamics](https://github.com/ddrous/gen-dynamics)). This is a call for fellow authors to upload their metrics and datasets, synthetic or otherwise, while following a consistent interface. In the context of OoD generalization for instance, we suggest the dataset be split in 4 parts: • (1) train: For In-Domain meta-training; • (2) test: For In-Domain evaluation; • (3) ood_train: For Out-of-Distribution adaptation to new environments (meta-testing); • (4) ood_test: For OoD evaluation

Each split should contain trajectories X and the time points t at which the states were recorded. The time t is a 1-dimensional array, and we recommend a 4-dimensional X tensor with dimensions described as follows: • (1) nb_envs: Number of distinct environments; • (2) nb_trajs_per_env: Number of trajectories per environment; • (3) nb_steps_per_traj: Number of time steps per trajectory (matching the size of t); • (4) state_size: Size of the state space

While the suggestions apply mostly to dynamical systems’ Meta-Learning, we believe they are generalizable to other problems, and are open to suggestions from the community. All problems described below are now represented in Gen-Dynamics.

### C.2 Simple Pendulum (SP)

The autonomous dynamical system at play here corresponds to a frictionless pendulum suspended from a stationary point by a string of fixed length. The state space $x=(\alpha,\omega)$, comprises the angle the pendulum makes with the vertical, and its angular velocity, respectively

$$ $\left\{\begin{aligned} \frac{{\text{d}\alpha}}{{\text{d}t}}&=\omega,\\ \frac{{\text{d}\omega}}{{\text{d}t}}&=-\frac{{\color[rgb]{1,0,0}g}}{L}\sin( \alpha).\end{aligned}\right.$ (14) $$

For this problem, each environment corresponds to a different gravity(^9^99To give an intuition behind the term “environment”, one might consider the surface of a celestial body in the solar system (see Section 1).), to which we don’t have access during the data collection stage. Therefore, with $L=1$ completely known, the goal is to learn a dynamical system that easily generalizes across the unobserved $g$. Only the NCF-$\mathit{t}_{1}$ variant is used in all experiments involving this problem.

The in-domain and adaptation mean squared error (MSE) metrics (see below) were reported in Table 2 during our favorable comparison to the One-For-All (OFA) baseline. Our loss values during training are reported in Figure 8. Unsurprisingly, learning one context-agnostic vector field for all environments (OFA) is suboptimal given the differences in gravity from one environment to the next. On the other hand, our method learns to discriminate between environments and produces accurate trajectories, one of which was presented in Figure 1b.

Here, we define the mean absolute percentage error (MAPE) and MSE criteria as they apply to trajectory data. These definitions are reused for subsequent problems.

$$ $\displaystyle\text{MAPE}(x,\hat{x})=\frac{1}{N\times d}\sum_{t=1}^{N}\left| \frac{x(t_{n})-\hat{x}(t_{n})}{x(t_{n})}\right|\times 100,$ (15) $$

and

$$ $\displaystyle\text{MSE}(x,\hat{x})=\frac{1}{N\times d}\sum_{t=1}^{N}\|x(t_{n}) -\hat{x}(t_{n})\|_{2}^{2}.$ (16) $$

Figure: Figure 8: MSE loss values when training a Neural Context Flow on the SP problem, compared with a baseline One-For-All formulation. OFA fails catastrophically since the diversity of environments in the training dataset prevents the approximation of any meaningful vector field.
Refer to caption: x6.png

Figure: Figure 9: Representation of the first and second dimensions of the learned contexts for the SP problem. The labels 1 to 25 identify the training environments, while $a_{1}$ and $a_{2}$ (in red), indicate the adaptation environments. We observe that environments close in indices (and thus in gravity) share similar contexts along these dimensions.
Refer to caption: x7.png

Beyond serving as a control signal for the vector fields, the contexts encode useful representations. In Figure 9, we visualize the first two dimensions of the various $\{\xi^{e}\}_{1\leq e\leq 25}$ after training. We observe that environments close in indices(^10^1010The indices of the training environments correspond to their ordering in increasing values of gravity.) are equally close in the two context dimensions. Similarly, distant environments are noticeably far apart in context space. The same observation is made during adaptation, where, for instance, $e^{\prime}=a_{2}$ (corresponding to $g=14.75$) gets a context close to $e=15$ (corresponding to $g=14.83$). This observation indicates that the latent context vector is encoding features related to gravity, which may be used for further downstream representation learning tasks.

### C.3 Lotka-Volterra (LV)

With dynamics that are of continued interest to many fields (epidemiology Venturino [], economy Wu et al. [], etc.), the Lotka-Volterra (LV) ODE models the evolution of the concentration of preys $x$ and predators $y$ in a closed ecosystem. The behavior of the system is controlled by the prey’s natural growth rate $\alpha$, the predation rate $\beta$, the predator’s increase rate upon consuming prey $\delta$, and the predator’s natural death rate $\gamma$

$$ $\displaystyle\begin{dcases}\frac{{\text{d}x}}{{\text{d}t}}&=\alpha x-\beta xy, \\ \frac{{\text{d}y}}{{\text{d}t}}&=\delta xy-\gamma y.\end{dcases}$ (17) $$

We repeat the experiment as designed in Kirchmeyer et al. []. All synthetic ground truth data is generated with the two initial states both following the $\mathcal{U}(1,3)$ distribution. Once sampled, we note that the same initial condition is used to generate trajectories for all environments. The parameters that vary across training environments are $\beta\in\{0.5,0.75,1\}$ and $\delta\in\{0.5,0.75,1\}$. In each training environment, we generate 4 trajectories with a Runge-Kutta time-adaptive 4th-order scheme, while we generate 32 for In-Domain evaluation. For one-shot adaptation, we extrapolate to $\beta\in\{0.625,1.125\}$ and $\delta\in\{0.625,1.125\}$, with only 1 trajectory per environment, and 32 for OoD testing. The observed parameters $\alpha$ and $\delta$ are always fixed at $0.5$.

The large grid-wise adaptation experiment conducted on the LV problem in Section 4.3 showed that NCFs were powerful at extrapolating. Here, through the MSE training losses in LABEL:fig:huge_adapt_2, we highlight the fact that the bottom and left edges of the grid display quite high errors. This illustrates the importance of observing several loss metrics when learning time series Hewamalage et al. [].

### C.4 Glycolytic-Oscillator (GO)

This equation defines the evolution of the concentration of 7 biochemical species $\{s_{i}\}_{i=1,2,\ldots,7}$ according to the ODE:

$$ $\displaystyle\frac{ds_{1}}{dt}$ $\displaystyle=J_{0}-k_{1}\frac{s_{1}s_{6}}{1+(s_{6}/K_{1})^{q}}$ $\displaystyle\frac{ds_{2}}{dt}$ $\displaystyle=2k_{1}\frac{s_{1}s_{6}}{1+(s_{6}/K_{1})^{q}}-k_{2}s_{2}(N-s_{5}) -k_{6}s_{2}s_{5}$ $\displaystyle\frac{ds_{3}}{dt}$ $\displaystyle=k_{2}s_{2}(N-s_{5})-k_{3}s_{3}(A-s_{6})$ $\displaystyle\frac{ds_{4}}{dt}$ $\displaystyle=k_{3}s_{3}(A-s_{6})-k_{4}s_{4}s_{5}-\kappa(s_{4}-s_{7})$ $\displaystyle\frac{ds_{5}}{dt}$ $\displaystyle=k_{2}s_{2}(N-s_{5})-k_{4}s_{4}s_{5}-k_{6}s_{2}s_{5}$ $\displaystyle\frac{ds_{6}}{dt}$ $\displaystyle=-2k_{1}\frac{s_{1}s_{6}}{1+(s_{6}/K_{1})^{q}}+2k_{3}s_{3}(A-s_{6 })-k_{5}s_{6}$ $\displaystyle\frac{ds_{7}}{dt}$ $\displaystyle=\psi\kappa(s_{4}-s_{7})-ks_{7}$ $$

where the parameters either vary or are set fixed as per Table 7.

**Table 7: Parameters and their values for the GO problem.**
| Parameter | $J_{0}$ | $k_{1}$ | $k_{2}$ | $k_{3}$ | $k_{4}$ | $k_{5}$ | $k_{6}$ | $K_{1}$ | $q$ | $N$ | $A$ | $\kappa$ | $\psi$ | $k$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Value | 2.5 | (varies) | 6 | 16 | 100 | 1.28 | 12 | (varies) | 4 | 1 | 4 | 13 | 0.1 | 1.8 |

Again, we follow the same procedure as in Kirchmeyer et al. [] to generate trajectories for each environment. Namely, we sample initial conditions from a specific distribution [Daniels and Nemenman, , Table 2]. We vary $k_{1}\in\{100,90,80\}$ and $K_{1}\in\{1,0.75,0.5\}$ to create 9 training environments with 32 trajectories each, both for InD training and testing. We use 1 trajectory for OoD adaptation to 4 environments defined by $k_{1}\in\{85,95\}$ and $K_{1}\in\{0.625,0.875\}$. We use 32 trajectories for OoD evaluation.

Intuitive post-processing of the contexts offers many insights into the meta-training process, particularly because it clusters the 9 environments into groups of 3, as illustrated in Figure 11.

Figure: Figure 11: Illustration of 256-dimensional context vectors consistently clustered with t-SNE embeddings with perplexity of 2 on the Glycolytic Oscillator (GO). The training environments are in green and labeled as their IDs, while the adaptation environments are in purple. This clustering mirrors the distribution of the system parameters.
Refer to caption: x8.png

### C.5 Sel’kov Model (SM)

We introduce the Sel’kov model in dimensionless form Strogatz [], a highly non-linear ODE mainly studied for its application to yeast glycolysis :

$$ $\displaystyle\frac{dx}{dt}$ $\displaystyle=-x+ay+x^{2}y$ $\displaystyle\frac{dy}{dt}$ $\displaystyle=b-ay-x^{2}y$ $$

where $a=0.1$ and $b$ is the parameter that changes. The trajectories are generated with a time horizon of 40, with 11 regularly spaced time steps. We sample each initial condition state from the distribution $\mathcal{U}\{0,3\}$, and we observe the appearance of a limit cycle (L1), then an equilibrium point (E), and another limit cycle as $b$ changes (see LABEL:fig:attractors).

Specifically, our 21 training environments are a union of 7 environments evenly distributed in $b\in[-1,-0.25]$ , then 7 evenly distributed in $b\in[-0.1,-0.1]$, and finally 7 others with $b\in[0.25,1]$ . We generate 4 trajectories for training and 4 for InD testing. As for adaptation, we choose 6 environments, with $b\in\{-1.25,-0.65,-0.05,0.02,0.6,1.2\}$ (see Figure 12). We set aside 1 trajectory for adaptation and 4 for OoD testing.

In Section 4, we showed that NCF-$\mathit{t}_{2}$ outperformed other adaptation rules. We note, however, that training such systems is very difficult, even with NCFs. We observed in practice that convergence is extremely sensitive to weight initialization, as the shaded regions on the loss curves of LABEL:fig:selkov_train_losses attest. The curves are complemented with Figure 13, emphasizing that the model performs best in environments near the equilibrium.

Figure: Figure 12: Visualization of a trajectory with the same initial condition in various environments of the Sel’kov Model’s dataset. In blue are the meta-training environments, in red the meta-testing ones.
Refer to caption: x9.png

Figure: Figure 13: Per-environment training and adaptatation losses. In blue are the meta-training environments, in red the meta-testing ones.
Refer to caption: x10.png

### C.6 Brusselator (BT)

The Brusselator model Prigogine and Lefever [] describes the reaction-diffusion dynamics of two chemical species, $U$ and $V$, and is given by the following system of partial differential equations (PDEs) defined on a $8\times 8$ grid:

$$ $\displaystyle\frac{\partial U}{\partial t}$ $\displaystyle=D_{u}\Delta U+A-(B+1)U+U^{2}V,$ $\displaystyle\frac{\partial V}{\partial t}$ $\displaystyle=D_{v}\Delta V+BU-U^{2}V,$ $$

where:

- •
$U$ and $V$ are the concentrations of the chemical reactants.
- •
$D_{u}=1$ and $D_{v}=0.1$ are the diffusion coefficients for $U$ and $V$, respectively.
- •
$A$ and $B$ are constants representing the rate parameters of the chemical reactions, both varying across environments.

Using an RK4 adaptive time-step solver, the system is simulated up to $T=10$ (excluded), with the step reported every time step $\Delta t=0.5$. The initial condition for each trajectory for all environments is sampled
and broadcasted as follows:

$$ $\displaystyle U_{0}$ $\displaystyle=\bar{A}$ $\displaystyle V_{0}$ $\displaystyle=\frac{\bar{B}}{\bar{A}}+0.1\eta_{ij}$ $$

where $\bar{A}\sim\mathcal{U}(0.5,2.0)$, $\bar{B}\sim\mathcal{U}(1.25,5.0)$, and $\eta_{ij}\sim\mathcal{N}(0,1)$ for each grid position $(i,j)$.

We aim to keep all environments involved in this problem outside the Brusselator’s oscillatory regime $B^{2}>1+A^{2}$. For Meta-training, $A$ and $B$ are selected from $\{0.75,1,1.25\}\times\{3.25,3.5,3.75\}$ yielding 12 environments, each with 4 trajectories for training, and 32 for InD testing. For adaptation, we select 12 environments from $\{0.875,1.125,1.375\}\times\{3.125,3.375,3.625,3.875\}$, with 1 trajectory for training, and 32 for OoD testing.

### C.7 Gray-Scott (GS)

We aim to validate the empirical potential of our method on spatiotemporal systems beyond ODEs. We thus identify the Gray-Scott (GS) model for reaction-diffusion, a partial differential equation (PDE) defined over a spatial $32\times 32$ grid:

$$ $\displaystyle\frac{\partial U}{\partial t}$ $\displaystyle=D_{u}\Delta U-UV^{2}+F(1-U)$ $\displaystyle\frac{\partial V}{\partial t}$ $\displaystyle=D_{v}\Delta V+UV^{2}-(F+k)V$ $$

Like in the Brusselator case above, $U$ and $V$ represent the concentrations of the two chemical components in the spatial domain with periodic boundary conditions. $D_{u},D_{v}$ denote their diffusion coefficients respectively, while $F$ and $k$ are the reaction parameters. We generate trajectories on a temporal grid with $\Delta t=40$ and temporal horizon $T=400$.

The parameters we use to generate our environments are the reaction parameters $F$ and $k$. We consider 4 environments for meta-training: $F\in\{0.30,0.39\},k\in\{0.058,0.062\}$; and 4 others for adaption: $F\in\{0.33,0.36\},k\in\{0.59,0.61\}$. Other simulation parameters, as well as the initial condition-generating distribution, are inherited from Kirchmeyer et al. [], where we direct readers for additional details.

### C.8 Navier-Stokes (NS)

Like the Gray-Scot case, the 2D incompressible Navier-Stokes case is inherited from Kirchmeyer et al. []. The PDE is defined on a $32\times 32$ spatial grid as:

$$ $\displaystyle\frac{\partial\omega}{\partial t}$ $\displaystyle=-v\nabla\omega+\nu\Delta\omega+f$ $\displaystyle\nabla v$ $\displaystyle=0$ $$

where:

- •
$\omega=\nabla\times v$ is the vorticity,
- •
$v$ is the velocity field,
- •
$\nu$ is the kinematic viscosity.

The trajectory data from $t=0$ to $T=10$ with $\Delta t=1$ is obtained through a custom Euler integration scheme. By varying the viscosity from $\nu\in\{8\cdot 10^{-4},9\cdot 10^{-4},1.0\cdot 10^{-3},1.1\cdot 10^{-3},1.2
\cdot 10^{-3}\}$ we gather 5 meta-training environments, each with 16 trajectories for training and 32 for testing. Similarly, we collect 4 adaptation environments with $\nu\in\{8.5\cdot 10^{-4},9.5\cdot 10^{-4},1.05\cdot 10^{-3},1.15\cdot 10^{-3}$ each with 1 trajectory for training, and 32 for testing. Other parameters of the simulation, as well as the initial condition generating distribution, are inherited from Kirchmeyer et al. [], where we encourage the readers to find more details.

## Appendix D Ablation Studies

These ablation studies are designed to investigate how the context pool size, the context size, and the 3-networks architecture affect the performance of NCFs.

### D.1 Limiting the context pool size

The context pool $\mathrm{P}$ from which environments $j$ are randomly sampled contributes significantly to the computational and memory complexity of the algorithm. As evidenced in Algorithms 3 and 1, the computation of the loss and hence its gradient can be parallelized across the contexts vectors $\xi^{j}$. With the goal of assessing the associated computational burden, we vary the size of the pool size for the LV problem from 1 to 9, reporting the In-Domain and OoD metrics, and computational time in LABEL:fig:contextpoolsize.

We dig further into the influence of the context pool size. While its ideal value might be difficult to conclude based on LABEL:fig:contextpoolsize, LABEL:fig:contextpoolsize_T1 investigates NCF-$\mathit{t}_{1}$ to highlight how important any value bigger than 1 is (with $p=6$ proving to be exceptionally adequate). Indeed, 0* indicates no context pool was used, and the vector field is evaluated without Taylor expansion. This ablation results in a catastrophic loss in MSE accuracy of about one order of magnitude. While computational time (see LABEL:fig:poolsizetimes) scales linearly with the pool size, the loss is not. This suggests that the context pool size $p$ is a hyperparameter that should be tuned for maximum accuracy.

All in all, the impact of the context size on the MSE loss values can be explained by the concept of “task relatedness”, or “context proximity”. The Taylor expansion is the key to information flowing from one environment to another. That said, some environments may be prohibitively far from others in context space. As such, the Taylor approximation incurs a non-reducible residual error term. Restricting the pool size limits the impact of those environments since they are less often sampled from the pool.

### D.2 Restriction of the context size

Rather than limiting the pool size, what if the context vectors themselves were limited? The latent context vectors are the building blocks of NCF, and having shown that they encode useful representations vital for downstream tasks in Section C.2, we now inquire as to how their size influences the overall learning performance. This further provides the opportunity to test the contextual self-modulation process. Indeed, over-parametrization should not degrade the performance of NCFs since context vectors must be automatically kept small and close to each other (in $L^{1}$ norm) for the Taylor approximations to be accurate.

Like in CoDA Kirchmeyer et al. [], the context size $d_{\xi}$ is directly related to the parameter count of the model; and limiting the parameter count bears practical importance for computational efficiency and interpretability. Thus, we perform the LV experiment as described in Section C.3 with $d_{\xi}\in\{2,4,8,16,32,64,128,256,512,1024\}$. The results observed in LABEL:subfig:spctsize align with our intuitive understanding of increased expressiveness with bigger latent vectors. Together, LABEL:subfig:spctsize and LABEL:subfig:loss_landscape shed light on the relationship between $\xi$ and the underlying physical parameters pair $(\beta,\delta)$, suggesting that NCFs would benefit from the vast body of research in representation learning.

We appreciate that while selecting $d_{\xi}$ small is more interpretable (as can be observed by the clustering of the training environment LABEL:subfig:loss_landscape), it is important to choose $d_{\xi}$ sufficiently big. Indeed, using the JVP-based implementation we provide in Appendix E, large context vectors come at a reduced extra cost in both speed and memory.

### D.3 Ablation of the 3-Networks Architecture

Another key element of the NCF framework is the 3-networks architecture described in Figure 2b. In it, the vector field consists of state and context-specific networks that help bring the inputs into the same latent representational space before passing to another network to approximate the local derivative. Like the context size, its ablation directly contributes to a reduction in parameter count.

For this experiment, we only consider NCF-$\mathit{t}_{1}$. We remove the data- and context-specific neural networks in the vector field, and we directly concatenate $x^{e}_{i}(t)$, and $\xi$; the result of which is passed to a (single) neural network. In this regard, NCF* (which designates NCF when its 3-networks architecture is removed) is analogous to CAVIA Zintgraf et al. []. That said, NCF still benefits from the flow of contextual information, which allows for self-modulation.

Figure: Figure 15: MSE loss curves when training the LV problem on a complete NCF, and on NCF*, the NCF variant deprived of the 3-networks architecture.
Refer to caption: extracted/5960501/figures/losses_ablation.png

**Table 8: MSE upon ablation of the 3-networks architecture (NCF*) on both SP and LV problems, highlighting difficulties adapting to new environments.**
|  | SP ($\times 10^{-3}$) | LV ($\times 10^{-5}$) |  |  |
| --- | --- | --- | --- | --- |
|  | In-Domain | Adaptation | In-Domain | Adaptation |
| NCF | 4.82 $\pm$ 0.08 | 0.25 $\pm$ 0.02 | 1.28 $\pm$ 0.25 | 0.8 $\pm$ 1.01 |
| NCF* | 8.39 $\pm$ 0.009 | 1239.64 $\pm$ 867.6 | 11.03 $\pm$ 1.1 | 69.98 $\pm$ 84.36 |

This study is of particular significance since, if in addition to foregoing the 3-networks architecture, we eliminated the Taylor expansion step, and NCF would turn into the data-controlled neural ODE Massaroli et al. []. So we run both dynamics forecasting problems without the 3-networks architecture, and we report the training MSE in Figure 15. In-domain and adaptation MSEs for both LV and SP problems are reported in Table 8.

While Figure 15 highlights a performance discrepancy of nearly 1 order of magnitude during training, the key insight is hidden in the adaptation columns of Table 8. Indeed, the removal of the data- and context-specific networks in the vector field considerably restricts the model’s ability to generalize to unseen environments, both for the SP and LV problems. This consolidates the 3-networks architecture as an essential piece of the NCF framework.

## Appendix E Example Implementation of NCFs

#### Vector field

The vector field takes center-stage when modeling using ODEs. It is critical to use JVPs, since we never want to materialize the Jacobian (the context size can be prohibitively large). For a highly performant implementation of neural ODEs, we leveraged JAX Bradbury et al. [] and its ever-growing ecosystem; in particular Optax DeepMind et al. [] for optimization, and Equinox Kidger and Garcia [] for neural network definition. The Jacobian-Vector Product primitive filter_jvp as illustrated below, coupled with jit-compilation, allowed commendable runtimes for all problems.

Figure: Listing 1: Second-order Taylor expansion in the NCF vector field.

#### Loss function

Equation 10 provides a structured summation layout particularly suited for a function transformation like JAX’s vmap. During implementation, that loss functions can be defined in two stages, the innermost summation term along with Equation 9 making up the first stage; while the two outermost summations of (10) make up the second stage. Once the first stage is complete, the second is relatively easy to implement. We provide a vectorized implementation of the first stage below.

Figure: Listing 2: Inner NCF loss function with vectorization support

## Appendix F Trajectories Visualization

Figure: Figure 16: Visualizing the first ground truth and predicted trajectories and phase spaces in all 4 adaptation environments of the LV problem. The initial condition is the same across all 4 adaptation environments.
Refer to caption: x11.png

Figure: Figure 17: Sample absolute reconstruction error of a trajectory during adaptation of the Gray-Scott system with NCF-$\mathit{t}_{2}$. This describes the sole trajectory in the fourth adaptation environment.
Refer to caption: x12.png

Figure: Figure 18: Sample absolute reconstruction error of a trajectory during adaptation of the Navier-Stokes system with NCF-$\mathit{t}_{2}$. The viscosity for this environment is $\nu=1.15\times 10^{-3}$.
Refer to caption: x13.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: Question:

- With respect to the 3-network model, can you remove the state-network but keep the context network? Does it make any difference compared to the one-network model which performs similar to CAVIA?


Suggestions:

- Include sample efficiency experiments for some example ODEs and PDEs. For example MSE for a model trained from scratch vs one finetunes via a new context vector.
- Include uncertainty as a function of forecast. One would assume that uncertainty increases as the forecast becomes longer. Could you provide such an estimate?
soundness: 3
strengths: - Writing is very clear and the methodology is well explained. This allows readers to understand the differences between this method and previous ones.
- Interesting use of context vectors through the 3-network model. Ablation studies in supplementary material show the need for such an architecture. 
- The combination via context through the Taylor expansion seems to be an interesting and novel application, which I can see being used in other fields outside of ODE and PDE simulations.
- The estimation of uncertainty via different context vectors is very simple yet very clear and useful.
summary: The paper introduces Neural Context Flows, a method for meta learning. The main contributions of the work are focusing on how to combine context vectors in a way that allows OoD generalization and interpretability.
weaknesses: Manuscript makes reference to sample efficiency of using such adaptive models for new context. However the manuscript does not include any experiments to support such a statement.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 2
questions: You introduce different context pool strategy (RA, NF, SF). What are the differences in terms of performance for your method? Some ablations should have been done to compare the different methods, e.g.:
- Comparing the performance (e.g. MSE, adaptation time) of RA, NF, and SF for different datasets.
- Analyzing how the choice of strategy impacts the learned context representations.
Then, propose guidelines for choosing the appropriate strategy.

You mentioned LEADS, a multi-task learning problem. It would have been nice to add this baseline to the different generalization experiments, especially for the new datasets that are not present in [3].

Can you provide training time details for the different meta-learning frameworks?
soundness: 3
strengths: The targeted problem is important. Building neural-ode like solvers able to generalize to changes in the PDE coefficients is important, often referred to solving parametric PDEs.

The method seems novel for learning dynamical systems with changes in pde coefficients. The use of taylor expansion is intuitive and natural. I particularly liked the intuition given line 202-205. Existing context-based methods do not try to leverage information from each context vectors, each describing the environment information. NCF fills this gap.

The method is evaluated on a wide range of PDE problems and is SOTA or competitive when considering a second order talyor expansion.
summary: This work proposes a new meta-learning strategy for learning dynamical systems governed by PDEs. It introduces a new multi-environment framework, where environments are defined by specific PDE coefficients, each describing a specific behavior. To do so, the paper proposes a Taylor expansion of a forecaster network at a context vector $\xi^e$ around other context vectors $\xi^j$. It thus collects information from other environments via a context pool of P indices, modulating the forecaster network and allowing also contexts themselves to be self-modulated. The method has other properties: parallelization, interpretability, uncertainty estimation, ... The method is evaluated both for in-domain and out-domain environments agianst two baselines (CAVIA & CODA) on multiple datasets and show competitive or sota performance.
weaknesses: Regarding the writing style of the paper:
- I think there is room for improvement. In the introduction, I think the problem of solving parametric PDEs / learning dynamical systems with varying PDE coefficients should be stated more clearly and explains what is an environment in your specific setting. The introduction should 1) clearly defines the problem of building generalizable neural PDE solvers, 2) what are the different directions taken to do so [1, 2, 3] and 3) explain how your work fits into these different directions and advances the field.
- There are especially 2 paragraphs, where neural ODEs and physics-based (hybrid) approaches are introduced, that are not necessary or too much detailed in my opinion.

The authors state that the method can be interpretable, provides uncertainty quantification and is parallelizable. These are important properties that lack for instance for CoDA, as you mentioned. I think that the paper should have exploit these properties in more depth, provide more ablation studies to show that NCF can exploit these properties such as:
- a detailed analysis showing how the learned context vectors relate to physical parameters, demonstrating interpretability
- benchmarks showing how parallelization impacts training time as the number of environments increases

[1] Subramanian et al.,  Towards Foundation Models for Scientific Machine Learning: Characterizing Scaling and Transfer Behavior, NeurIPS, 2023.

[2] Takamoto et al., Learning Neural PDE Solvers with Parameter-Guided Channel Attention, ICML, 2023.

[3] Kirchmeyer et al., Generalizing to New Physical Systems via Context-Informed Dynamics Model, ICML, 2022.

---

confidence: 3
contribution: 2
initial_rating: 1
presentation: 2
questions: 1) Did you used **different** trajectories to compute the metrics in the tables than the one used to perform the adaptation in OoD setup?
2) Could you please explain why your model performs better in OoD than InD ? This is an unexpected resuls : we would expect the model to perform always better on seen data.
3) Could you clarify what *robust* means in this context ? (robustness to OoD, to noise, to unseen initial condition ? ) and point to the corresponding experiment that support the claim ?
4) Could you elaborate on the motivation of using Taylor expansion for training the model ? If the main reason behind this is to force context to remain close, then a crucial ablation is missing where the Taylor expansion is replaced by a simple distance loss between context vector. 
5) The context vector is a learned embedding, so what is the point of learning a supplementary dedicated network to convert it into $\tilde \xi$, since you could directly learn $\tilde \xi$ ? In the paper, you explains that "this allows the framework to automatically balance the potentially non linear influence of the context with that of the state vector". What does this means ? Could you please elaborate on this ?
6) Could you please justify the claim line 248 that $L^1$ regularization is crucial ? Did I miss the ablation of this regularization in the paper ?
7) I am not sure to understand fig. 6. It seems that you collected the prediction of NCF for a same initial condition with different context vectors (including unrelated ones) and consider the max/min of these prediction as uncertainty measure. Hence my questions:
	- What is the probability of the true solution to lie within these bounds ?
	- It seems that all model gives fairly similar outputs, which is surprising since they are adapted to different dynamics. How do you explain that ?

**Minor**
 - I do appreciate code release. May I ask you to document and clean the code before publication ?
 - Could you please fix the url link to the associated github repo in the paper ?
 - Could clarify if Theorem 3.2 should be considered as a contribution ? It seems to be a straightforward application of Theorem 2 in Li et al, 2019. If not, it might be interesting to provide more details about what changed, and for which reason.

# Rebuttal update
The authors have provided a thorough reply to most of my concerns and have clearly demonstrated that my suspicion of overfitting was unfounded. This significantly impacts my evaluation, and I have raised my score from 1 to 5.

However, I still recommend rejecting this paper. While the proposed method is interesting and well-supported, it would benefit from further development, particularly in the analysis of its behavior. Although the authors made commendable efforts to improve their manuscript, the short timeframe allocated for rebuttals limited the depth of the supplementary results included. Specifically:

- The emergence of clusters of "similar" environments due to the context pool-filling method, coupled with the Taylor expansion constraint, is a compelling idea. However, it warrants deeper exploration, perhaps with more environments and a detailed analysis of the content of these clusters. These behavior is at the heart of the success of NCF and deserve a better, deeper exploration.
- Some design choices remain unclear, particularly the use of a projection network on the context embeddings.
- The authors have strengthened the "uncertainty"-related experiments, but this discussion should have been made clearer in the main paper. However, addressing this would require significant restructuring of the paper, which is impractical within the rebuttal timeline.

For these reasons, I have updated my evaluation to a borderline rejection. While the paper presents interesting results, I believe it is not yet ready for acceptance.
soundness: 2
strengths: The paper introduces a method that seems reasonable, with several interesting analysis of the behavior of the learned vector field. The related works section does present the most popular baselines and methods for this problem. The datasets used to evaluate the method are in par with what is currently used in the litterature.
summary: The paper introduces a new method for meta-learning of dynamical systems by enforcing context vector to be smooth and close to each other using a training method based on Taylor expansion of the vector field.
weaknesses: I have several major concerns on this paper that I will try to rank from higher to lower importance.
1) **Overfit**: I am very confused by the OoD adaption protocol used in the paper. From what I understood, Algorithm 2 is used on the validation set to tune the value of $\xi$ using gradient descent. Then, the prediction error is computed **on the same trajectories as the one used to pick $\xi$**. Given the size of $\xi$, it is most certainly overfitting to the (small) set of trajectories used to tune the context vector. I suspect that the value of $\xi$ is specifically tuned to match a small set of trajectories during validation, hence the results are not representative of the true performance of the model. A quick look at the code seems to confirm that, although it is hardly readable (many commented functions, no comments, empty README, residual files from project development). Moreover, the model performs almost systematically better in OoD setup than in In-domain, (table 1 and 2), which is unexpected, and very unusual.
2) **Unsupported claims**: 
	- "we introduce a *robust*  ..." ? The word "robust" is never used in the experimental section.
	- line 248 : "*The $L^1$ regularization term is particularly vital in promoting sparsity in possibly high dimensional context vector*". Sparsity is indeed used a lot for equation retrieval, such as SinDy-like methods, but I don't see why this is crucial here. Yet it probably is, but I don't see which experiment in the paper justify this.
	- "*straightforward method for uncertainty quantification*" : seems like an overstatement according to my understanding of the experiment. Quantifying the uncertainty means that the model can provides an interval to which the true trajectory must belong (up to some probability score). From what I understand from your experiment, you consider that the predictions from your model obtained with different contexts (including unrelated to the current environment) gives an uncertainty on the prediction.
1) **Motivation of the method**: the use of Taylor development is motivated (section 3, l.200 to 206) by the fact that it forces context vector to remain close. It is not clear to me why closeness is important for the task, and more importantly why simpler regularization to force $\xi_i$ to remain close from each other would not perform well.
2) **Complexity of the method**: training NCF requires to compute the second order derivative of a neural network, and then back-propagate through the entire computation graph, including the use of a Neural ODE. The method seems utterly complex and computationally demanding. Moreover, some design choices are unclear to me (see questions). Many hyper-parameters have to be set, and the paper provides no clear indication on how to set them (mostly, size of the context vector and context pool mode (RA, NF or SF)). However, I did appreciated the experiments on the size of the context pool. 

	Finally, Figure 7 and 16 shows training curves of the model exhibiting high spikes, including one (figure 7) from which the model never recover. This seems to indicate a highly unstable and difficult training, which is not a good sign for extrapolation to more complex dynamics.
	
1) **Experimental section**: It seems that the order 1 NCF model is trained with a different algorithm than the order 2 (l.264-268). I am not sure to see why, and more importantly, if the gap between the two models in table 1 is due to the supplementary order of the Taylor expansion, or the different training algorithm.

**Minor** (no effect on rating)
- in the introduction: "*Its dynamics are heavily depending on its parameter*" looks like an overstatement. The dependance of a dynamical system on its parameter can vary significantly from a system to another.
-  In section 2, you mentioned that CoDA involves two networks instead of one, hence requiring more computational resources to train. Your approach uses three networks. You might clarify this statement to explicitly mentioned hyper-networks as the bottleneck, and not the use of two networks ?
- in the beginning of section 3 (l.168), $\mathcal{D}_{ad}$ is mentioned before being introduced.
- L. 171, you refer to the smoothness assumption as an inductive bias. It's closer to a constraint than a real inductive bias (or at least, it is an inductive bias on the type of dataset you will test your model on, but not an inductive bias of general physical systems).

**Motivation for my grade**
My grade is mostly driven by my suspicion of overfitting. However, I am also concerned by the (in my opinion) poorly motivated design choices (use of Taylor expansion to promote closeness of context vector, supplementary encoder for context vector, back-propagating through the Hessian) and my difficulty to understand several claims and experiments.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: Refer to weakness.
soundness: 3
strengths: (1) Introducing a meta-learning framework to study complex dynamical systems driven by ODEs and PDEs is a novel and interesting direction.
(2) This article includes both experimental validation and some theoretical analysis.
summary: This paper proposes NCF to solve parametric ordinary differential equations. It is designed to improve the adaptability and generalization of learning dynamical systems across various environments. NCF introduces a meta-learning approach that employs a context-modulation mechanism, incorporating uncertainty estimation. Specificaly,  NCF uses k-th order Taylor expansion to enable contextual self-modulation. The authors demonstrate the performance of NCF on a variety of ODEs and PDEs problems and illustrate the effectiveness of the proposed methods.
weaknesses: (1) Currently, there are several methods for parameterizing equations, such as those for ordinary differential equations [1] and partial differential equations [2]. In the paper, the authors mainly compared their approach with meta-learning methods like CAVIA and CODA. So, does it have a competitive advantage over other non-meta-learning methods?
[1] Parameterized Neural Ordinary Differential Equations: Applications to Computational Physics Problems
[2] Identification of the flux function of nonlinear conservation laws with variable parameters

(2) The authors mentioned that their method is robust. One important aspect of proving robustness is how the model performs when the observed data contains noise. However, this was not demonstrated in the experiments.

(3) The authors converted three PDE problems into ODE problems for their study. They should provide a detailed description of how this conversion was done and how the errors were controlled during the process. PDE problems are quite sensitive to the choice of numerical schemes, and different discretization methods can significantly affect the accuracy of the solution. Converting them into ODEs and solving them using an ODE solver will inevitably introduce errors.

(4) The link to the code provided in the paper is inaccessible.

Meta-review: This paper proposes a novel meta-learning strategy for modeling dynamical systems across varying environment parameters. The core idea is to introduce environment-specific latent context vectors and expand the vector field using a Taylor series about these context vectors. This approach facilitates information sharing across environments, enhancing data efficiency and improving adaptation to new environments. The method is evaluated for both in-domain and out-of-domain environments against several state-of-the-art meta-learning baselines on multiple datasets, demonstrating competitive performance.

All reviewers agree that the targeted problem is important, and the use of Taylor expansion is intuitive and well-supported. Given the promising applications in meta-learning and the compelling results, I recommend this paper for publication.

Decision: Accept (Poster)