System-Aware Neural ODE Processes for Few-Shot Bayesian Optimization
Authors: Jixiang Qing
Rebecca D. Langdon
Robert Matthew Lee
Behrang Shafei
Mark van der Wilk
Calvin Tsay
Ruth Misener
Venue: Accepted by TMLR ()

Abstract: We consider the problem of optimizing initial conditions and termination time in dynamical systems governed by unknown ordinary differential equations (ODEs), where evaluating different initial conditions is costly and the state's value can not be measured in real-time but only with a delay while the measuring device processes the sample. To identify the optimal conditions in limited trials, we introduce a few-shot Bayesian Optimization (BO) framework based on the system's prior information. At the core of our approach is the System-Aware Neural ODE Processes (SANODEP), an extension of Neural ODE Processes (NODEP) designed to meta-learn ODE systems from multiple trajectories using a novel context embedding block. We further develop a two-stage BO framework to effectively incorporate search space constraints, enabling efficient optimization of both initial conditions and observation timings. We conduct extensive experiments showcasing SANODEP's potential for few-shot BO within dynamical systems. We also explore SANODEP's adaptability to varying levels of prior information, highlighting the trade-off between prior flexibility and model fitting accuracy.

Published in Transactions on Machine Learning Research (5/2025)
System-Aware Neural ODE Processes for Few-Shot Bayesian Optimization
Jixiang Qing'; Becky Langdon', Mark van der Wilk?, Calvin Tsay',
Robert M Lee², Behrang Shafei? Ruth Misener
"mmeerial College London
BASFF SE Univvesiity of Oxford
Reviewed on OpenReview: htips: // openreview. net/ forum? id= FFARLuWefK
Abstract
We consider the problem of optimizing initial conditions and termination time in dynamical systems governed by unknown ordinary differential equations (ODEs), where evaluating different initial conditions is costly and the state's value can not be measured in real-time but only with a delay while the measuring device processes the sample. To identify the optimal conditions in limited trials, we introduce a few-shot Bayesian Optimization (BO) framework based on the system's prior information. At the core of our approach is the System-Award Neural ODE Processes (SANODEP), an extension of Neural ODE Processes (NODEP) designed to meta-learn ODE systems from multiple trajectories using novel context embedding block. We further develop two-stage BO framework to effectively incorporate search space constraints, enabling efficient optimization of both initial conditions and observation timings. We conduct extensive experiments showcasing SANODEP's potential for few-shot BO within dynamical systems. We also explore SANODEP's adaptability to varying levels of prior information, highlighting the trade-off between prior flexibility and model fitting accuracy.
1 Introduction
This paper studies widely encountered, yet less investigated, problem: optimizing the initial conditions and termination time in unknown dynamical systems where evaluations are computationally expensive. We assume that the primary evaluation cost comes from switching initial conditions and wish to use as few trajectories as possible.
This issue is prevalent in multiple fields, including biotechnology, epidemiology, ecology, and chemistry. For instance, consider the optimization of (bio)chemical reactor. Here, the objective is to determine the optimal "recipe," " i.e,, the set of initial reactant concent selectivity, and/or minimize waste [Taylor et al. 2023: Schoepfer et al. 2024 Schilter et al., 2024 Thebelt
and reaction times that. maximiie yield, enhance
et al. 2022). Specifically, the high costs associated with changing reactants for multiple experimental runs highlight the need for developing efficient optimization algorithms (Paulson Tsay 2024).
Bayesian Optimization (BO) (Frazier 2018 Garnett 2023) is well-established method for optimizing expensive-to--valutee black-box objective functions. It relies on probabilistic surrogate models built from limited function evaluations to guide the optimization efficiently. However, standard Gaussian Processes (GPs) (Williams & Rasmussen 2006) (the de facto probabilistic surrogate model in BO) with conventional kernels do not capture the dynamics of these systems effectively. Recent efforts to introduce Bayesian ODE models aim to incorporate suitable assumptions; however, their inference often involves lime-consuming computations (Heinonen et al.
2018 Dandekar et al.) 2020), making them impractical for a time- sensisev
optimization scenario. Some Bayesian ODE models attempt to alleviate computational burdens by leveraging
crude approximation inference (Out al. 2023 or approximating numerical integrators (Hegde al. 2022
"Corresponding author: j. j.ging@imperial..cco

---

Published in Transactions on Machine Learning Research (5/2025)
Objective Function
2.0
Initial Condition
1.5
SANODA NP
1.0
-2
0.5
0.0
Time 1.5
Number of Evaluated Trajectories
Figure 1: Illustration of meta-learning based few-shot Bayesian Optimization (BO) with time-delay constraints (detailed Section 4) applicd to the Kotka-Volteryy (LV) system using three different models: Gaussian Process (GP), Neural Process (NP), and System Aware Neural ODE Process (SANODEP) (Section. Starting with one randomly sampled trajectory (marked with a o), The BO progressively selects promising initial conditions for experiments (Section. 4.1 and, based these, recommends subsequent mcasurement times while enforcing minimum time delay constraint (At) (Section. 4.2 The figure demonstrates that after evaluating just. one trajectory (the first, left-most column of the three rows), the meta-learned SANODEP model more closely approximates the original LV system compared to the non-meta-learnet GP model. This enhanced approximation facilitates the BO search in identifying more promising initial conditions.
Ensinger et al., 2024). While these approaches reduce computation time, the necessary approximations problem of
may degrade model performance. Consequently, the widely encountered yet performing BO in dynamical systems remains open, primarily due to the lack of appropriate model.
This work attempts to tackle this optimization problem by conceptualizing it as few-shot optimization, through leveraging prior information about dynamical systems to formulate the meta tasks to train a learning model, which is then able to adapt to new problems with very few data (a.k.a,, shots) and served as our BO surrogate model. To achieve this, we focus on the Neural ODE Process (NODEP) (Norcliffe et al. 2021) as our learning model due lo several useful model properties. First, NODEP combines Neural ODEs [Chen al. 2018 emphasizing the dynamical system perspective that may provide more informative inductive bias with Neural Processes (NP) |Garnclo et al. 2018b)). enhancing the meta-learning aspect for few-shot optimization. Second, meta-learning endows NODEP with fast adaptation to new observations, mitigating the potential concern of model training time with incoming measurement results. Third, NODEP enables continuous gradient-barn based optimization over time, compared with discrete-time meta-learn models (c.g., Foong et al. (2020)), integrating scamlessly within the BO framework. Finally, compared to other continuous-tinn time models involving computationally intensive operations (e.g., time domain attention (Chen)
al. 2024 NODEP is computationally efficient, making it suitable for BO.
However, building meta-learning based few-shot BO framework for unknown dynamical systems is not a straightforward downstream application. An obstacle is that NODEP's functionality is limited by its original learning objective: predicting a single trajectory state value given content data for that trajectory. For initial condition optimization of unknown ODE systems, it essential to learn to aware of the underlying governing systems using context data consisting of one or several different trajectory observations, as this benefits for forecasting trajectories starting from arbitrary initial conditions. Moreover, no BO frameworks are proposed for such problem settings.
In line with our optimization requirements, we introduce the System-Aware Neural ODE Process (SANODEP), generalization of NODEP in terms of meta-learning of ODE systems to plug-in our newly developed BО framework for few-shot BO in dynamical systems. Figure shows, through example, SANODEP outcomes in comparison to other approaches. Our contributions are as follows:
1. We propose novel context embedding mechanism enabling SANODEP to meta-learn from batch of trajectories with minimal model structure adjustments based on NODEP.
2

---

Published in Transactions on Machine Learning Research (5/2025)
2. We developed model-agnostic, time-delay constraint process BO framework for optimizing initial conditions and termination time in dynamical systems.
3. We compare SANODEP and other meta-learned models under our BO framework, validating the benefits of the ODE-aware model structure and demonstrating the effectiveness of the few-shot BO framework's model agnosticity property.
4. We conduct initial investigation how the different levels of prior information can be utilized for SANODEP. Strong prior information can enable physically-informed model structure with extended loss considering parameter inference, enabling novel few-shol parameler estimation functionalily, while weak prior information may still be useful through a properly designed task distribution, albeit. at a detriment to model fitting capability.
The rest of the paper is organized as follows: In Section 2 we describe the preliminaries. Scction B introduces the Swstem-Aware Neural ODE Processes (SANODEP). In Section 4 we develop the optimization framework based on SANODEP specific to our optimization problem. Section 5 presents the numerical experiments we conducted. In Section
discusses related work. Section 6 we investigate the impact of different levels
of prior information.
2
Problem Statement and Background
2.1 State Optimization
Consider dynamical system whose evolution is described by the following ordinary differential equation (ODE) system in the time domain, denoted t:
dæ f(x h), x<h = to) Zo.
(1)
dh
where stt) E Xe C Rd represents the state value of the system at time t, f()) : Rd x R + Rd is an unknown function (vector field) governing the system dynamics, and Dy E Xy S Xx C Rdd the initial system state. Define fevolve (xF t, f) : Xo x (R x R - - Rd= where T E [to, tmazl, to obtain the system state at time from an initial condition Ty as: fevolve (Io, f) = xo + S.e f(=(h), g(-) 1e R be practitioner-specinn consider the following multi-objective problem:
Let (known) function that aggregates the state values into scalar,
(1.00)erxx) ((fevolue(xo) t, f)), -f))),t
max
(2
The goal to identify the initial conditions wo along with the corresponding evolution termination time that provide the optimal trade-off between the objective 9 applied to the state values æ(t) and the amount of time to reach this state. Simply put, we wish to maximize the objective function early.
2.2 Neural ODE Processes (NODEP)
NODEP (Norcliffe et al. 2021 is latent variable-based Bayesian heta-learning model with generative processcs that can be summarized as follows: given set of N observations of state values at different times from single trajectory, represented by context set C := {(ff, > NODEP assumess the conditional prediction has been generated from a latent controlled ODE state dimensionality de- The stochasticity of the model is induced by stochastic latent initial condition Lo and the stochastic dynamics representation, which in practice is implemented via a time-invariant control LD. These two terms are obtained through first encoding the context elements d ((Cf, f)) and then applying average pooling (Zaheer et al. 2017) to produce a single representation vector r = N E ((tf,, which then mapped to corresponding distributions Lo and LDp Once we sample the latent initial condition and the trajectory dynamics representation 3 the
With slight abuse of notation, follow the common notation in control to use (Instead of d) to represent the realization of time-invariant control Lp.
3

---

Published in Transactions on Machine Learning Research (5/2025)
evolution of the latent ODE is:
/"' Jan (1(h), u, h Oode) ) dh. Lg
[()= = Sevatue (la, 1, (un(u, Onde)) = ((to) +
(3)
Here fon represents the vector field as neural network parameterized by Oode- We compute fevotve using mumerical ODE solvers. To predict the state values at time t, after solving the latent ODE, the latent state ((t) E IR4 can be decoded back to real state predictions through a decoder p(a((t)((t), t). Thus, the overall joint probability of the model can be represented as:
p (u, to x t) -p ( u, p(u((())(())))))
(4)
The latent state's dimensionality de typically chosen be larger than dæ as additional state variables benefit the model flexibility (Dupont et al. 2019 and enable learning of higher-order dynamics implicitly (Noreliffe et al. 2020).
Given a context set C and a target set T== 1(4 x DD to calculate the log- log-likelihood on T for inference, the intractable posterior of (4o[C) p(u has been approximated through the mean field approximation g(La g(l//T) (T))( which practice has been implemented through encoder qr (lo|C), qn (u (C) an amortized fashion, eventually leading the following evidence lower bound (ELBO):
log (((=1)) 11 1.....
(5)
uTT)
Σ log p (20x to, (x)
KL (an. (lo/T) llrr. (l,(C) KL [gn (ulT) ||yo (u(C)),
where T := ((f)) represents the (irregularly sampled) target times NODEP has competitive perfor-
mance in predicting single trajectory, and we therefore consider extending its functionality to meta-learn ODE system distributions for our subsequent optimization purpose.
3 System Aware Neural ODE Process
TC
x
TC new
x
Few-Shot Optimization with Prior Information Given that evaluating fevotve(Wa, L, f) is limited by cost to small number of different initial conditions. To enable fast adaptation with few evaluations, we take a meta- lcarning approach by assuming that f is a rcalization a stochastic function F, and we can access data from its different realizations to formulate the meta-tast distribu- tions.
Loaew
D Days
TT Thee
XT
More precisely, consider a random function F: Xo x dd representing the distribution of dynamical systems. From specific realization of F, suppose we have observed M distinct trajectories, the context set then encompasses observations from all M trajectories, denoted as T (TF ..., T Each trajectory, labeled T,C, includes its own set of N₁ context observations: T.F = ((tf x1 N Furthermore, consider new trajectory Taew- After ob- serving additional context observations on this trajectory, TC (e.g., the initial condition of this new trajectory), can define the extended context set as C= 7 U T new
Figure 2: Graphical Model of SANODEP, the model predicts any time point in the new tra- jectory by knowing both observations from new trajectories and from past trajectories. Depend- ing on whether the X and T consists of more than the initial condition, the model focuses on forecasting or interpolating tasks. The solid and dashed lines represent the generative and in- ference processes, respectively. convention by assuming is superset of C.
We follow the Garnelo al. 2018b and Norcliffe et all (0021 We denote context trajectory data either as a set of pairs {TC. a()) 1 or as vectors (IC,X) x²), where 2 = [ef
C of T. We similar notation (I", x² for the target set.

---

Published in Transactions on Machine Learning Research (5/2025)
Additionally, the target set is defined as T = new a visual claboration is provided in Figure 5 in Appendix A
It is straightforward to see why NODEP is suboptimal in such scenario as it makes predictions based on single time series (C- = Tnew T = =new), unable leverage M trajectorys' information. Below, we demonstrate how SANODEP efficiently enabled through the set-based representation (Zaheer et al. 2017
3.1 Set-based Dynamical System Representation
Similar to the latent variable u (Eq. 3) capturing trajectory dynamics in NODEP, we adapt a latent variable Days g(Ways/CC (C) in SANODEP, ellectively replacing u in the model structure but with enhanced conditioning on context observations from M + trajectories, to capture the dynamical systems properties. Follow this mechanism, will have a feature extraction for multi-start multivariate irregular time series, in addition with an efficiency requirement for our optimization purposes.
Such questions are less considered in contemporary time series models, as existing approaches (Shukla & Marlin 2021; Schirmer et al., 2022) are mainly developed for single initial condition start trajectories. Specifically for multi-start scenario, Jiang al. (2023) proposes first extracting trajectory- wise aggregated feature vector, and then averaging feature vectors final representation. However, both extraction are implemented through convolution operation, which poses challenges with irregular time series without additional modifications.
A straightforward thinking would be still using average pooling across all context elements as a set-based representation. While the Picard-Lindclöf theorem [Lindelff) 1894) guarantecs the uniqueness of state values in an initial value problem when fax 1) is Lipschitz continuous, in case when only part of the underline system states can be measured as xt), different trajectories might still share identical æ(t) values, leading to identifiability issues (duplicate context elements come from different initial conditions). Consequently, augment each observation with its corresponding initial condition: (4, To, ,; to enhance the model's s ability to differentiate between trajectories that might otherwise appear identical. Then we perform the average pooling on the flattened context set to obtain the context representation Tays 1-1 N; En Em ( 21 rf,, which is then mapped to a distribution
propose an augmented sets-based approach:
representing possible system realizations from F that has generated C. Aside from the theoretical intuition, also empirically compare the average pooling without initial condition augmentation in the Section. 6.1 and find that the augmentation leads to more robust performance even F only represents first-order systems with fully observable states.
For the rest of the model, SANODEP follows the NODEP structure, hence requiring minimal adjustment, only with additional care on activation function choices to enable differentiability (Appendix A.3). Appendix A.2 provides detailed description of the model structure.
3.2 Bi-scenario Loss Function
In line with the principles of episode learning (Vinyals et al. 2016), SANODEP's training is structured through multiple episodes. Motivated by our optimization problem that will be detailod in Section. the design of each episode's problem ensures good model performance under two primary scenarios:
1. Forecasting: Using M context trajectories TC, the model predicts future state values XT new for a new trajectory initiated from T {(Inew) C xnews)) at. designated target times TT
new
2. Interpolating: From the same M trajectories, the model interpolates and extrapolates state values for a new trajectory that already includes K > observations: T = {tnnew) Inew; K predicting the states at times T new-
Note that again Th. T
5

---

Published in Transactions on Machine Learning Research (5/2025)
Both scenarios can be considered under bi-scenario loss function, designed to enhance model's accuracy for predicting new trajectory states X new at Tew
La FM.T.,, farecast: P fareona) log pe (x new! I U T forecast)) Tnew)
(6)
where po(.) represents the SANODEP prediction when parameterized by
Algorithm Learning and Inference in System Aware Neural ODE Processes (SANODEP)
Require: ODE system inducing distributions P, known trajectory range [Miin: Marl) batch size of Monte Carlo approximation of initial condition sample Na,, batch size of dynamical systems sample Naysr prespecified time grid grrid linspace(to: Emaxs grrid), minimum and maximum context points within trajectory Mmiin Mmax. minimum and maximum extra target points mmmn Amaz, initial condition space Xy.
1: Initialize SANODEP model parameters with random seeds.
2: for step in training_steps do
3: # Training data generation
for j= to Nays do
Sample ODE system f F and Nm, different initial conditions from Xo; solve Na, odes at time grids Tarid to obtain the dataset (Ti,..,, TN Sample known trajectory number M from Uniform Uniform(Mmin M.nan.
6:
7:
for | = to M do
8:
# Sample context and target elements within each trajectory Sample my from Uniform(Mmnn Mmaa), ) sample ni from Uniform(Ammnn Aman). Randomly subsample from Ti to extract the context dataset T. and target dataset Ti",, where
10:
|T my and [T: my + ny and Ti C T.T
11:
end for
12:
Concatenate M trajectories context set T = (TT, ... TW) and larget set T (TF ... Til end for
13:
14:
# Model Prediction
15:
for j =
Neys do for 1= No. do
16:
17:
sample forecast from Bernoulli(A) # determine forecasting or interpolating if forecast then
18:
19:
T The = (1los: Ror)) where (log: vog) E Te else if forecast = O and & >> M then Randomly subsample from Ta to extract the context dataset Th and target dataset Tx. where [T IT. = m, [T m + n and To C C T
20:
21:
22:
else
23:
T = Th -
24:
end if
25:
Augmented the context trajectorics and target trajectorics with the new trajectory T UTT
UT
TT=TFUTE Thew
26:
Compute variational posterior q(usgel TT), q( g(Lo (TT), g( g(Weys (T5,, q( q(Lo (T) through the encoder block.
27:
Sample 1(4)), Mays from ((((1)) q q(tsyss777 (T)).
28:
Solve latent odes as in Eq. (3) for all times and decode to obtain the model prediction.
29:
Calculate the trajectory wise loss LELBOR based on Eq.
30:
end for
31:
end for
Average the trajectory wise loss LELBO No, E Neo LELBOO
32:
33:
Update the model parameter through optimizer +- - 17 V end for
34:
During training, dynamical system realization f from F has been drawn, randomly sample the mumber of context trajectories M that we have observed, where the number of context elements, the initial

---

Published in Transactions on Machine Learning Research (5/2025)
condition of the trajectory and the observation time are all sampled. For the new trajectory to be predicted, besides sampling the target set (T X a Bernoulli indicator 1 forecost ~ Bernoulli(A) parameterize by À is sampled to determine whether the episode will address a forocasting or interpolating sccnario. This indicator directly influences the part of the context set sampled from the new trajectory
- [(tnew) C Kneww)) {)
if 1 forecast = 1,
The foreeast) = [[(fnes) Fnew: Fnnew M forecast 0.
For each system realization, we train in a mini-batch way by making prediction on a batch of different new trajectories.
The intractable log- ligliikeldoood in Eq. 6 can be approximated via the following evidence lower bound (ELBO):
log po (X new
T U The (1 forecast) Thew)
E q(u.ys ITEUT (1 (heeee
((65 log P (Xnew T Thew Uays D
(8)
IT UTU Them (1 forecast) U Tnew) ll 9 (usys - KL (q (LJ... : x xboww) p( (nnLn )
- KL [q] (Usys
( U Thew (Iforcastt)
where q (usys T U Tee new (1 forecast) U T and q (I 1 ($6 (n Tf....) (..... has been obtained through the encoder in amortized fashion, the prior p( (_..... is isotropic Gaussian. Appendix B derives the ELBO and Algorithm. provides implementation details of the model inference procedure. Compared to Eq. 5 Eq. (S) is designed to generalize across any new trajectory within the same dynamical system. This aligns with our Bayesian Optimization (BO) setting, where we further introduce the BO algorithm in the next section.
4 Time Delay Constraint Process Bayesian Optimization
Through maximizing the ELBO for 0, we can obtain SANODEP's predictive distribution po( T .o) for batch of time T ..., (N) at specified initial condition To, which is sufficient for few-shot learning tasks. In the section, specifically for our practically motivated optimization problem that incorporates one additional time delay constraint, we propose an optimization framework, leveraging SANODEP for few-shot BO in ODE and benchmark in Section. 6.2
Minimum Observation Delay Constraint When optimizing Eq. (2 we assume that, while one can observe state values any chosen time on specific initial condition wa as illustrated in Figure D the next observation can only commence after fixed known time duration Al. Tn practice, these delays stem from the need to sequentially conduct separate, smaller experiments for state value measurements, each requiring known period to complete before the next can begin.
With the above context, the proposed optimization framework consists of the following two steps:
4.1 Initial Condition Identification
The first optimization stage identifies the optimal initial conditions necessary for starting experiments. The optimality of the initial conditions is defined as achieving the maximum expected reward after the completion of observations starting at this location. Inspired by the batch strategy achieve a similar non-myopic objective (González el. al.; 2016 Jiang et al. 2020), we propose an adaptive batch size based optimization strategy for the initial condition identification:
maximizen, (aaxtniiicnnei....tt/aa--veeeee (ttrr a (Zo. to..t.,.......... tN.,,NN (C,,T..2))
s.t. Vi [1...., NJ - At
where a is the batch acquisition function to be maximized, Nopt = (1, N is the search space for the total number of observation queries, with Nmax = [(mmax - to)/ Atj. T (4,,..,,yy ,,,] represents the set of observation times. We omit the dependence of (=[[ C, L, mo) in a(.) for brevity thereafter.
Forr notation simplicity, always subscript represent the initial condition of trajectory.

---

Published in Transactions on Machine Learning Research (5/2025)
4.2
Choose the next query time
Once the imitial condition has been chosen, and our last observation query made is at time tas when still have query opportunity (i.e., fmar > in + At), redefine the maximum remaining trajectory observation Nmar J, the batch size search space Nope and time search space [tn + At, (marl, we choose the next query time recurrently via the following optimization problem:
88 Eime--ay A
(areniieee..........eree...e..e s.t. Vi E (l,... NJ : - 1--1 At.
(Ion, 11, 1..,,,,.1,1 (N)
(10)
Search Space Reduction The integer variable N's search space Nopt in Eq. 9.. 10 though one dimensional, can be cumbersome to optimize in practice. However, for the batch expected hypervolume improvement
acquisition function (QEHVI) (Daulton et al. 2020) that will use as acquisition function the search space can be reduced without losing optimality:
Proposition 1. For acquisition functions defined as a() = Exx [HVI ( F C)]
maximnie
#c alzo, (I...,t = maximize
Bo (2.... EeT, Nef NNaee !m..)
a(xo, ,,,.,,,,,, tN),
(2. Ne
where HVT stands for hyperuolume improvement based on Pareto frontier F², see Definition 2 of Daulton et al. (020). Appendix C.I provides the proof and shows that the search space reduction property holds for generic acquisition functions that are monotonic w.r.l. set inclusions. Consequently, we define Nopt as mmn N, maz] thereafter.
2
A Algorithm Model Assisted Ordinary Differential Equation Optimization Framework
Input: maximum number of experiment trajectories Nexp; minimumn spacing At, maximum observations per trajectory Nmaz initial evaluated context datasets C, acquisition function a, evaluated trajectory number 1, initial time lower bound to tor time range [fib maa]) observer fevolue for measurement.
2: while S Neap do
3: ==
4: Construct the probabilistic model for pt C, xn)
5: # Initial Condition Optimization
27. N* arg max no m.
Nagur 1.N a (a.,t
t. Vi E (1, ..., NJ 4 - 1--1 Al
4
tu = to + At, = Lb, maz)) T = T Utu 35) C = CUT
7: Neraj = Nmaa -1
8: # Within trajectory Next Measurement Time Scheduling
while Ntraj >= do
10:
N (i arg max
(b, INE.) Niras N.Nerasl a (li,..,,NN) 4-11 At
Vi e (1...., N)
11:
i Sevotve (2 (1)
12: to + At, =
bh, twbl, C= CC (1)
13:
Ntraj
At
-in
14: end while
15: n = 1 #1
16: end while
Optimization Framework outlined Algorithm. We refer to the Appendix C.1 for details of how SANODEP utilized within QEHVI, together with acquisition function optimizer described in Appendix C.I Finally, highlight that this Bayesian Optimization framework is agnostic to the model choice of (X (X" will compare in Section. both NP and GP (as non-meta-leat model) can be used.
C, T To),
8

---

Published in Transactions on Machine Learning Research (5/2025)
5
Related Work
Meta-learning of Dynamical Systems: Our work is built upon a meta-learned continuous time-basee model for state prediction in dynamical systems. Singh et al. (019) extended Neural Processes (Garnelo et al,, 2018b) by incorporating temporal information from the perspective of the state-space model. However, their approach is limited to discrete-time systems. Following the NODEP framework, [Dag et al. (2023) extended the model for forecasting purposes in multi-modal scenarios. Recently, [Jiang et al. (2023) used meta-learning to tackle high-dimensional time series forecasting problems for generic sequential latent variable models, they leverage apatial-tempone convolution to extract the Usys, which only works on regularly spaced time steps (e.g,, image frames in their use cases). However, SANODEP extracts system dynamic's feature from multiple irregularly sampled multivariate time series using the newly proposed encoder structure. Sonn & Jeong 2023 explored Hamiltonian representations, which are flexible enough for cross cross-domain generalizations. Beyond the Bayesian meta-lcarning paradigm, Li et al. (2023) used classical gradient-basee
leaanien Finn
et al. 2017) lo meta-learn dynamics with shared properties. Ausina et al. (2024) investigated the separation and modeling of both dynamic variables that influence state evolution and static variables that correspond to invariant properties, enhancing model performance. Motivated by recent advances in extending Graph Neural Networks (GNNs) to dynamical systems, several studies have explored leveraging GNN-ODEs for modeling multi-agent interactions (e.g-. [uo et al. (2025); Yuan et al.] (2024)). However, such advantages is not applicable to single-agent problem setting. Additionally, several works (Luo et al.) 2025 Fotiadis et al. 2023) investigate robustness to distribution shifts in system-level dynamics by separating system-leval and trajectory-level dynamics encoding. While this separation improves trajectory-level prediction robustness, it is less beneficial to guard against trajectory-level predictive performance if the initial condition optimization is already misleading due to distribution shift.
Bayesian Optimization: Vellanki et al. (017 address process-constraingg batch optimization where each batch shares identical constrained variables within a subspacc. Our problem similarly treats initial conditions as constrained variables in batch opt but differs by using a joint expected utility method instead of greedy strategy. Additionally, our approach involves on-the-fly optimization with system evolution, leading to novel timo-constrained problem addressed by adaptive batch BO algorithm with search space reduction. Folch et al. (2022 2024 also consider optimization under movement (initial condition) constraints.
Few-Shot Bayesian Optimization: A significant amount of recent work has focused on meta-learning stochastic processes, to be used as alternatives to GP for BO. Within the (Conditional) Neural Process framework ((C)NP) [Garnelo al. 2018a b). Transformer Neural Processes (TNP) (Nguyen & Grover 2022
have demonstrated robust uncertainty quantification capabilities for BO. B0, Dutordoir et al. (2023) introduced a
diffusion model-based Neural Process, providing novel framework for modeling stochastic processes which enables joint sampling and inverse design functionality. Similar to TNP, Müller et al. (2021) approached the meta-learning problem as Bayesian inference, subsequently integrating a BO framework into this model Müller et al. 2023). We note that while a majority of the aforementioned approaches lead to promising empirical performance on BO for regular black- box functions, their backend models (if there are any) are mainly motivated by a drop-in replacement of GP or non-tempore stochastic process. meta-learned BO method specifically for dynamical systems is, to our best of knowledge, not revealed yet.
6 Experiments
This section conducts experiments on meta-learning and few-shot BO for dynamical systems. All models are implemented using Flax (Heek et al. 2023) and are open source, available in: https://github.ccm TsingGAg/SANODEP The optimization framework is based on Trieste (Picheny al. 2023).
6.1 Modeling Comparison
We evaluate model performance comparisons on meta-learning diverse ODE system distributions with a twofold purpose. First, investigate whether leveraging the system information can help modeling. We also investigate whether the inductive bias could bring performance benefits.
9

---

Published in Transactions on Machine Learning Research (5/2025)
Baseline: To see if adding system information helps with the modeling, we compare SANODEP against NODEP under its original testing scenario (interpolation). For the second purpose, we compare SANODEP with Neural Processos (NP) (Garnclo et al.
2018bp. As mentioned in Section D we were not able to compare with the meta ODE models of Jiang el al.
(2023) as it not been able to be applied on irregular time series.
Although GP-ODE (Heinonen al. 2018 appears to be a promising non- non-meta-lee meta-learrning baseline for modeling dynamical systems, it is not directly applicable to our meta-testing scenario. This is due to its prohibitively expensive training process and the fact that its performance on modeling multiple irregularly sampled trajectories from the same dynamical system has not been investigated in the original work. Therefore, we instead provide results using a standard Gaussian Process with independent outputs for reference, where we use Matérn 3/2 kornel as a common choice (c.g., Xie et al. (2024); Back et al. (2024)) and trained with maximum likelihood estimation [ For models that do not explicitly incorporate the ODE in their structure, we use the trajectory's initial condition augmented with time as its imput. We also note that due to the context set consisting of batches of multiple trajectories sampled irregularly, prevents efficient use of scaled dot product bttention-based models (e.g.. Nguyen & Grover (2022) (see Appendix
B.2 for details) hence we
omit for comparing these none ODE based NP models.
Training: Following Norcliffe et al. 2021), we treat F as a parametric function of a specific kinetic model with stochasticity induced by model parameter distributions P 7 We refer to Appendix [D1 for the model formulations, its parameter distribution P's form, and the training epochs. Interpolation-only experiments use SANODEP trained using Eq. ( setting À 0, and NODEP trained using Eq. (5), We refer to Appendix. [.2 for the specific training parameter setting for Algorithm. D
[
For the second experimental purpose, we train all models using Eq. (6) setting À 0 for interpolation only training and À 0.5 for bi-scenario training, with 5 different random sceds for model initialization.
We evaluate the model performance on target data generated from systems sampled from the same F, which is known as in-distribution generalization evaluation. Excluding GP, each model was evaluated on 10 random systems, cach consisting of 100 trajectorics to prodict in a minibatch fashion. As it computationally infeasible to evaluate GPs on the same scale, we used a random subset of the test set, 5,000 trajectories. Appendix B.2 provides full data generation settings. Following Norcliffe et al. (2021), we assess model performance in terms of Mean-Squared Error by conducting sequential evaluation with varying numbers of context trajectories, as illustrated in Figure 3 (and Figure g in Appendix D.2 for Negative-log- Likelioood....... We also provide the model prediction illustration in Appendix D.3 and the averagdd.onttet--trajctttyyyyi performance results in Tables 2 and 3 in Appendix D.2
Evaluation:
Table 1: Model comparison on interpolation tasks (meean-squared-error x10-2) for range of dynamical systems.
Models
Lotks- Volterra (2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
Lotka- Volterra (3d)
SIRD (4d)
NODEP-A- - 0 38.9 i 0.52
15.5 1 0.379
7.12 1 0.56
367.5 1 58.0
26.8 1.53
188.1 21.8
SANODEP-A== =0 32.1 + 0.55
8.7 + 1.11
0.52 + 0.03
254.3 + 60.3
25.1 + 1.85
137.2 + 12.8
Does System Awareness Improve Predictive Performance? Table. D demonstrates that SANODEP is better than NODEP in all problems, indicating the benefit of capturing system information.
Does Incorporating Temporal Information Enhance Model Performance? Illustrated in Figure comparing with NP, except. for the Lotka-Voterra problem where the SANODEP trained with mixing probability À = 0.5 shows larger predictive variance, and Lotka- Voterra (3d) where NP demonstrates better performance when context trajectory numbers are small, for the rest of the problems, SANODEP cither on par or noticeably better than NP irrespective with the mixing probability. Indicating that the incorporation of right. inductive bias has the potential to help modeling.
Does Augmenting Initial Conditions Improve Modeling? SANODEP by comparing models with and without initial condition augmentation, as depicted in Figure g
We further evaluate the performance of
We use the GP implementation within Trieste and additionally follow the common numerical trick implemented in linear_ iinear_peerator" (used by CPyTorch), where the jitter term is progressively increased if Cholesky decomposition fails.
We discuss Appendix [ on how the flexibility of F may change and affect the SANODEP.
10

---

Published in Transactions on Machine Learning Research (5/2025)
Lotka-Voterra (2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
Lotka-Voterra (3d)
SIRD (4d)
10a
10"
10
10"
10a
10a
10
MSE
10
10
0.4
10
0.5
10a
0.4
0.3
10
10
10"
10
0.2
5
10
10
0
10
10
10
5
10
Number of Context Trajectories
SANODEP-KKr - 0.5
NP-X
NP-A = 0.5
NODEP
SANODEP-A D.0
SANODEP-A 0.5
Figure 3: Model evaluation performance comparison different number of context trajectories. The first row corresponds to forecasting performance (prediction with only the initial condition known), and the second row represents an interpolating setting. It can be seen that SANODEP is either on-par or marginally better than NP for most problems, and the initial condition augmented encoder-based model variant provides more robust performance across problems.
(labeled as SANODEP- for without augmentation). The results intriguingly show that even for first order systems, models augmented with initial conditions exhibit greater robustness across various systems. This demonstrates the empirical benefits of our context aggregation approach.
How Do Meta-Learned Models Compare with Non-Meta-Learnee Models? We provide the sequential model evaluation plot in Figure including GPs in Appendix D.2 As expected, the meta-learned models show superior accuracy when the number of trajectories is small (as we have also illustrated in Figure aligning with the problem setting of this study. We also remark that GP can quickly overtake when the context trajetcory number abundant, since the Neural Process typed model is known to suffer from under fit when the data is abundant.
6.2 Few-Shot Bayesian Optimization
Lotka Vitte-Voteeraaaao (2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
(3d)
SIRD (4d)
Hypervolume Indicator
30
1.5
2.6
0.10
30
2.4
0.05
20
1.0
0.1
10
10
0.00
0.0
0.5
Scaled Time
Scaled Time
Scaled Time
Scaled Time
Scaled Time
Scaled Time
NP-A 0.0
NP-1 0.5
NODEP
SANODEP-A
SANODEP-A 0.5
Figure 4: Comparison of few-shot BO performances in terms of mcan + standard deviation. Excopt for the Lotka-Voterra (3d) problem, BANODEP-basedd few-shot BO demonstrates competitive performance on all the rest of the problems.
Baselines: We conduct few-shot BO using the meta-learned models (SANODEP, NODEP, NP) trained in the previous section using our optimization framework (ark((Allu (Algorithm 2 For reference, we provide GP-based BO (GP-BO) using the same framework. Since the optimization problem we consider is novel and not addressed by existing BO algorithms, we conduct our comparisons within the proposed BO framework, evaluating different models under a consistent setting.
Optimization:
We optimize dynamical system realizations sampled from F. The parameter settings, the objective function definitions and performance indicator calculation are provided in Table For simplicity, in each optimization problem, assume the maximum number of observations per trajectory is At fuog -
11

---

Published in Transactions on Machine Learning Research (5/2025)
We start with one randomly sampled trajectory with uniformly spaced observations and have an additional budget to query 10 trajectories. Each optimization is repeated 10 times with different random seeds. Results: are reported in terms of the scaled experimental time versus the hypervolume indicator, with the reference point also provided in Table (we also report the results of PI-SANODEP that will be introduced in next section.).
Docs Meta-lcarning Help Few-Shot Optimization in Dynamical Systems? We highlight the merits of few-shot BO, as it consistenlly demonstrates signilicant convergence speed improvements during the early stages of BO.
Does Dynamical System Informed Model Behaves Better than Standard Model for Optimization? Exept for Lotka-Voterra (3d), SANODEP demonsttates better performance compared with Neural Processes. Again validate the potential of a deliberate consideration of temperal information.
Does Weight Training Objective help Optimization? For SANODEP, except in the case of the Brusselator, combining training with bi-scenario losses shows either comparable slightly better performance than training solely in the interpolation setting. This suggests that careful design of the loss function for optimization purposes may offer modest benefits.
Finally, we provide the time profile of different models in the Appendix D.5] in general, as the involvement of simulation in inference process, SANODEP and NODEP takes more time than NP and GP when perform optimization.
7 Discussion
Bayesian Optimization (BO) of unknown dynamical systems is an underexplored area, with existing approaches often relying on unsuitable surrogale models. Addressing this gap, while minimising the number of costly optimization steps, we have developed a novel few-shot BO framework. Our approach extends the Neural ODE Processes (NODEP) into the System-Aware Neural ODE Processes (SANODEP), which meta-learns the system's prior information to enhance the optimization process. Through extensive benchmark experiments, demonstrrated SANODEP's potential over NODEP. Our results show that SANODEP, equipped with ortinization-driven loss function, offers competitive performance compared to other non-dynamic-aware meta-learning models, and can offer additional possibilties including parameter estimation which can't be said with other non-dynamic-aware meta-learning models.
Limitations: A key limitation of SANODEP, inherited from Neural Processes, is the issue of underfitting. This necessitates a trade-off between few-shot functionality and model fitting capability. While strong prior information enhances performance in scenarios with well-defined priors, the model performance is limited by lack of prior information about the dynamical system. Although We have discussed extensively variety of prior information in Appendix L and our proposed task distribution (Appendix E.2 shows promise for cross-domain generalization, SANODEP still faces trade-off between model fitting and prior flexibility, indicating that it cannot be used out of the box as a generic probabilistic model to optimize arbitrary dynamical system without knowing its kinetic form. Furthermore, unlike GP-BO which has established sublinear cumulative regret bounds typically due to the connection between maximum information gain and model predictive variance from the kernel (Srinivas et al. 2009), amortized latent variable models like SANODEP face several upstream open issues, including the convergence properties of amortized variational inference. Consequently, regret bounds for SANODEP based BO remain an open.
the
Future investigations: Future work will aim to develop more adaptable conditional distributions to better capture a wider range of dynamical systems. This will address the trade-off between model performance and the flexibility needed to handle diverse and complex dynamical systems without relying heavily on strong priors.
Acknowledgments
The authors gratefully acknowledge support from BASF SE, Ludwigshafen am Rhein (BL), Engineering and Physical Sciences Research Council [grants EP EP/W003317/1, EP EP/X025292/1, EP BP/Y028775/1, and
12

---

Published in Transactions on Machine Learning Research (5/2025)
EP EP/S023151/1] (RM, CT, JQ, BL), a BASF /RAEng Research Chair in Pata-Driven Optimisation (RM), and a BASF/RAEng Senior Research Fellowship (CT).
References
Tize Amanda Auzina, Çagatay Yildız, Sara Magliacane, Matthias Bethge, and Efstratios Gavves. Modulated neural ODEs. NewrIPS, 36, 2024.
Jaewan Back, Jinwoo Kim, Hyunho Lee, Minki Lee, and Mingi Choi. Bayesin-optimizztoon--sss efficient operation for direct ammonia solid oxide fuel cells. Journal of Power Sources, 619:235194 94, 2024.
Ricky TQ Chen, Yulia, Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In NewrIPS, volume 31, 2018.
Yuqi Chen, Kan Ren, Yansen Wang, Yuchen Fang, Weiwei Sun, and Dongshong Li. ContiFormer: Contimuous- time transformer for irregular time series modeling. In NeurIPS, volume 36, 2024.
Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation based inference. Proceedings of the National Academy of Sciences, (17(48):30055-3000) 2020.
Raj Dandekar, Karen Chung, Vaibhav Dixit, Mohamed Tarek, Aslan Carcia-Valadie Krishna Vishal Vemula, and Chris Rackauckas. Bayesian neural ordinary differential equations. arXiv preprint orXic:2012.07244, 2020.
Ting Dang, Jing Han, Tong Xia, Erika Bondareva, Chloë Siegele-Brown, Jagmohan Chauhan, Andreas Grammenos, Dimitris Spathis, Pietro Cicuta, and Cecilia Mascolo. Conditional neural ODE processes for individual discase progression forocasting: case study on COVID-19. In 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 3914-3925, 2023.
Samuel Daulton, Maximilian Balandat, and Eytan Bakshy. Differentiable expected hypervolume improvement for parallel multi-objective Bayesian optimization. In NeurIPS, volume 33, pp. 9851-9864, 2020.
Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural ODEs. In NewrlPS, volume 32, 2019.
Vincent Dutordoir, Alan Saul, Zoubin Ghahramani, and Fergus Simpson. Neural diffusion processes. ICML, pp. 10000-902,200 2023.
Katharina Ensinger, Nicholas Tagliapietra, Sebastian Ziesche, and Sebastian Trimpe. Exact inference for contimuous-time Gaussian process dynamics. In AAAI, volume 38, pp. 11883-11891, 2024.
Lco Feng, Hossein Hajimirsadeghi, Yoshua. Bengio, and Mohamed Osama Ahmed. Latent bottlenecked attentive neural processes. In ICLR, 2023.
Chelsea Finn, Picter Abbeel, and Sergey Levine. Model-agnostic meta-learning for fast adaptation of deep networks. In ICML, pp. 1126-1135, 2017.
Jose Pablo Folch, Shiqiang Zhang, Robert Lee, Behrang Shafei, David Walz, Calvin Tsay, Mark van der Wilk, and Ruth Misener. SnAKe: Bayesian optimization with pathwise exploration. In NewrIPS, volume 35, pp. 35226 35239, 2022.
Jose Pablo Folch, Calvin Tsay, Robert M Lee, Behrang Shafei, Weronika Ormaniec, Andreas Krause, Mark van der Wilk, Ruth Misener, and Mojmir Mutaý. Transition constrained Bayesian optimization via Markov decision processes. In NearIPS, volume 37, 2024.
Andrew Foong, Wessel Bruinsma, Jonathan Gordon, Yann Dubois, James Requeima, and Richard Turner. Meta-learning stationary stochastic process prediction with convolutional neural processes. In NewrIPS, volume 33, pp. 8284-8295, 2020.
13

---

Published in Transactions on Machine Learning Research (5/2025)
Stathi Fotiadis, Mario Lino Valencia, Shunlong Hu, Stef Garasto, Chris D Cantwell, and Anil Anthony Bharath. Disentangled generative models for robust prediction of system dynamics. In International Conference on Machine Learning, pp. 10222-10248. PMLR, 2023.
Peter I Frazier. A tutorial on bayesian optimization. arviv preprint ar arXi:: 2018.
Marta Garnelo, Dan Rosenbaum, Christopher Maddison, Tiago Ramalho, David Saxton, Murray Shanahan, Yee Whye Tch, Danilo Rezende, and SM Ali Eslami. Conditional neural processes. In ICML, pp- 1704-1713, 2018a.
Marta Garnelo, Jonathan Schwarz, Dan Rosenbaum, Fabio Viola, Danilo Rezende, SM Eslami, and Yee Whye Teh. Neural processes. arXiv preprint arX:::1807.01:1 0122,,, 2018b.
Roman Garnett. Bayesian optimization. Cambridge University Press, 2023.
Javier González, Michael Osborne, and Neil Lawrence. GLASSES: Relieving the myopia of Bayesian oplimisation. In AISTATS, pp. 790-799, 2016.
Jonathan Heek, Anselm Levskaya, Avital Oliver, Marvin Ritter, Bertrand Rondepierre, Andreas Steiner, and Marc van Zee. Flax: A neural network library and ecosystem for JAX, 2023. URL htttp///itthub.com
google/flax
Pashupati Hegde, Çagatay Yildiz, Harri Lähdesmäki, Samuel Kaski, and Markus Heinonen. Variational multiple shooting for Bayesian ODEs with Gaussian processes. In UAI, pp. 790 799, 2022.
Markus Heinonen, Cagatay Yildiz, Henrik Mannerström, Jukka Intosalmi, and Harri Lähdesmäki. Learning unknown ODE models with Gaussian processes. In JCML, pp. 1959 1968, 2018
Shali Jiang, Henry Chai, Javier Gonzalez, and Roman Garnett. Binoculars for efficient, nonmyopic sequential experimental design. In ICML, pp. 4794 4803, 2020.
Xiajun Jiang, Ryan Missol, Zhiyuan Li, and Linwei Wang. Sequential latent variable models for few-shot high-dimensional time-series forecasting. In ICLR, volume 11, 2023.
Patrick Kidger. On neural diffcrential equations. arXiv preprint Xxxx:2202.02435 2022.
Hyunjik Kim, Andriy Mnih, Jonathan Schwarz, Marta Garnelo, Ali Eslami, Dan Rosenbaum, Oriol Vinyals, and Yee Whye Tch. Attentive neural processes. In ICLR, 2019.
Qiaofeng Li, Tianyi Wang, Vwani Roychowdhury, and Mohammad Khalid Jawed. Metalearning generalizable dynamics from trajectories. Physical Review Letters, 131(6))066301, 2023.
Ernest Lindelöf. Sur l'application de la méthode des approximations successives aux équations différentielles ordinaires du premier ordre. Comptes rendus hebdomadaires des séances de l'Académée des sciences, 116 (3):454-457, 1894.
Xiao Luo, Yiyang Gu, Huiyu Jiang, Hang Zhou, Jinsheng Huang, Wei Ju, Zhiping Xiao, Ming Zhang, and Yizhou Sun. Pgode: towards high-quality system dynamics modeling. In Proceedings of the 41st International Conference on Machine Learning, ICML'24. JMLR.org, 2025.
Samuel Müller, Noah Hollmann, Scbastian Pineda Arango, Josif Grabocka, and Frank Hutter. Transformers do Bayesian inference. In ICLR, 2021.
Samuel Müllor, Matthias Feurer, Noah Hollmann, and Frank Hutter. PFNs4BO: In-context learning for Bayesian optimization. In ICML, pp. 25444 25470, 2023.
Tung Nguyen and Aditya Grover. Transformer neural processes: Uncertainty-awer meta learning via sequence modeling. In JCML, pp. 16569-16500 2022.
Alexander Norcliffe, Cristian Bodnar, Ben Day, Nikola Simidjievski, and Pietro Liò. On second order behaviour in augmented neural ODEs. In NewrIPS, volume 33, pp. 2020.
14

---

Published in Transactions on Machine Learning Research (5/2025)
Alexander Norclife, Cristian Bodnar, Ben Day, Jacob Moss, and Pietro Liò. Neural ODE processes. In ICLR, 2021.
Derek Onken and Lars Ruthotto. Discretize-optimize vs. optimize-diseretize for time-series regression and continuous normalizing flows. arXie preprint orXiv:2005.13..0 2020.
Katharina Ott, Michael Tiemann, and Philipp Hennig. Uncertainty and structure in neural ordinary differential equations. arXie prepriul arXiv:2305.1322 2023.
Joel A Paulson and Calvin Tsay. Bayesian optimization as a flexible and efficient design framework for sustainable process systems. arXiv preprint orXiv:2401.16375 2024.
Victor Picheny, Joel Berkeley, Henry B Moss, Hrvoje Stojic, Uri Granta, Sebastian W Ober, Artem Artemev, Khurram Ghani, Alexander Goodall, Andrei Paleyes, et al. Triesto: Efficiently exploring the depths of black-box functions with TensorFlow. Xix preprinl orXie:2302.08/p. 2023.
Thomas Pinder and Daniel Dodd. GPJax: A Gaussian process framework in JAX. Journal of Open Source Software, 7(755:4455, 2022. doi: 10.21105/joss.04455 URL https:://bii.rgg// ottp:/06...r//111.221/////////// 003...4..55555
Ilya Prigogine and René Lefever. Symmetry breaking instabilities in dissipative systems. ii. The Journal of Chemical Physics, 48 44):1166050,000, 1968.
Jixiang Qing, Tom Dhaene, and Ivo Couckuyt. Spectral representation of robustness mcasures for optimization under input uncertainly. In ICML, volume 38, 2022.
Jixiang Qing, Henry Moss, Tom Dhaene, and Ivo Couckuyt. {PF) ees: Parallel feasible pareto frontier entropy search for multi-objective Bayesian optimization. In AISTATS, 2023.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. NIPS, 20, 2007.
Oliver Schilter, Daniel Pacheco Gutierrez, Linnea M Folkmann, Alessandro Castrogiovanni, Alberto García-
Durán, Federico Zipoli, Loic M Roch, and Toodoro Laino. Combining Bayesian optimization and automation to simultaneously optimize reaction conditions and routes. Chemical Science, 2024.
Mona Schirmer, Mazin Eltayeb, Stefan Lessmann, and Maja Rudolph. Modeling irregular time series with continuous recurrent units. In ICML, pp. 19388-19405, 2022.
Alexandre Schoepfer, Jan Weinreich, Ruben Laplaza, Jerome Waser, and Clemence Corminboeuf. Cost- informed Bayesian reaction optimization. 2024. doi: /he2444//hemmxxi-220000
Satya Narayan Shukla and Benjamin M Marlin. Multi-time attention networks for irregularly sampled time series. In ICLR, 2021.
Gautam Singh, Jaesik Yoon, Youngsung Son, and Sungjin Ahn. Sequential neural processes. In NewrIPS, volume 32, 2019.
Yeongwoo Song and Hawoong Jeong. Towards cross domain generalization of Hamiltonian representation via meta learning. In ICLR, volume 12, 2023.
Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3993 2009.
Connor J Taylor, Alexander Pomberger, Kobi C Felton, Rachel Grainger, Magda Barecka, Thomas W Chamberlain, Richard A Bourne, Christopher Johnson, and Alexei A Lapkin. brief introduction to chemical reaction optimization. Chemical Reviews, 12336)83 3126, 2023.
Alexander Thebelt, Johannes Wiebe, Jan Kronqvist, Calvin Tsay, and Ruth Misener. Maximizing information from chemical engineering data sets: Applications to machine learning. Chemical Engineering Science, 252: 117469, 2022.
15

---

Published in Transactions on Machine Learning Research (5/2025)
Pratibha Vellanki, Santu Rana, Sunil Gupta, David Rubin, Alessandra Sutti, Thomas Dorin, Murray Height, Paul Sanders, and Svetha Venkatesh. Process-constrained batch Bayesian optimisation. NearIPS, 30, 2017.
Oriol Vinyals, Charles Blundell, Timothy Lillierap, Daan Wierstra, et al. Matching networks for shot. learning. In NewrIPS, volume 29, 2016.
Christopher KT Williams and Carl Edward Rasmussen, Gaussian processes for machine learning, volume MIT press Cambridge, MA, 2006.
Mohan Wu and Martin Lysy. Data-adaptive probabilistic likelihood approximation for ordinary differential equations. In AISTATS, pp. 1018-1026, 2024.
Yilin Xie, Shiqiang Zhang, Joel Paulson, and Calvin Tsay, Global optimization of gaussian process acquisition functions using a piecewise-lineal kernel approximation. ar preprinl 1rXiv:2410.16880 2024.
Jingyang Yuan, Gongbo Sun, Zhiping Xiao, Hang Zhou, Xiao Luo, Junyu Luo, Yusheng Zhao, Wei Ju, and Ming Zhang. Egode: An event-attended graph ode framework for modeling rigid dynamics. In The Tiirth-eighth Annual Conference on Neural Information Processing Systems, 2024.
Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola. Deep sets In NeurIPS, volume 30, 2017.
16

---

Published in Transactions on Machine Learning Research (5/2025)
Appendix
Table of Contents
A Problem Setting and Model Structures
17
A.1 Problem Setting Illustrations
17
A.2 Model Structure
17
A.3 On the Differentiability of the Encoder and Decoder
19
B
ELBO Derivation and Training Details
20
SANODEP EL derivation
20
B.2 Learning and Prediction Process of SANOEP
20
Time Delay Constraint Optimization
21
C.1
Acquisition Functions
21
Search Space Reduction
D
Experimental Details and Additional Results
22
D.1
Mela Training Data Definition
22
D.2
Additional Experimental Results
23
D.3
Model Prediction Visual Comparison
................................
23
D4
Optimization Problem Definition
23
D.5
Inference Time
24
E
Invsstggtionn of the prior strength, }
26
E.1
Strong Prior information: Plysicc-Informed Informed SANODEP
26
E2 Weaker Prior of F: Cross-domain genera lization possi bi
29
Notations
31
A Problem Setting and Model Structures
A.1
Problem Setting Illustrations
A.2
Model Structure
The model structure and hyperparameters of SANODEP
summarized as follows:
Initial State Encoder:
Tinit demit(((oo æol), Dinit IR++1 RT consists of three dense layers and activation functions in between.
Augmented State Encoder: r; (1i,,, Ca mil), or = -R-22+++ 2d+++1 Rr consists of three dense layers and activation functions in between.
System Context Aggregation: Tays N Σ
Context to Hidden Representation heys daye (Teys), hanas (oinse(Tsniil Deys()) and dina(.) one dense layer followed by context encoder activation function, haye 12th Minit 12
17

---

Published in Transactions on Machine Learning Research (5/2025)
3D
T
TC new
T22
T
fnewi]
E
Epes N
Figure 5: Illustration of the initial condition optimization problem within one state variable (d 1) ODE system. The optimization involves forecasting scenario: having observed the context set C == (U T,9) U. new (including the new trajectory's initial condition, illustrated +), one tries to predict the target state values X [newl : Tnew at target times TT
nea [tnew ..... knew Nvl (illustrated *).
Hidden to Variational posterior of Usys: ((usyal+) N (Dnage (heys). diag (go sys (hsye) )
12))
PHays Rh > Rdva is a dense layer, Do == R Rdvva is defined as ow + softplus( softplus(Dense()) (Dense(.)) Tih hyperparameter.
Hidden to Variational posterior of Lo: q( (e) N (Ou (hmmt), ) P imit (m))))?). Puanat Rh R' and doinit R R' has the same structure as above, but with independent weights.
dODELl, Usys. t| + 1,, ΦODE == Bo++ - Rd,, ΦODE is consisted of three dense layers with nonlinear activation functions in between.
ODE block:
Decoder: Qdee (t), Usys. f] xx),, L= /' (Ørdec [(t), Usys, t), d [t),, Usya f)), where Ordec and dee has the same structure variational posterior, but with independent weights.
Model Hyperparameters
Encoder o, output dimension r: 50
Context encoder (Ør, Pinit) activation function: SiLU
Encoder hidden dimension h: 50
ODE layer (ΦODE) activation function: Tanh
ODE layer hidden dimension: 50
Decoder hidden dimension: 50
Latent ODE state dimension I: 10
Context to latent dynamics (Ouager sys))) activation function: SiLU
Context to latent initial (Ørense Peinitt condition activation function: SiLU
Latent dynamics dimension days: 45
Variational posterior q ((Lozewl.) variance lower bound: dib 0.1
Variational posterior q(Usy=l.) variance lower bound: du = 0.1
Decoder (Øéecs &pece dosee) activation function: SiLU
18

---

Published in Transactions on Machine Learning Research (5/2025)
ODE solver: Dopri5 with rtol = le --- 5 and atol - 1e-5
We use SiLU instead of ReLU as in the original NODEP is to enforce the differentiability, which we claborate on in the following subsection.
A.3 On the Differentiability of the Encoder and Decoder
Since the optimization takes the initial condition together with the time as decision variables, we need the model output to be differentiable w.r.t. these quantities. NODEP overlooked this part and utilized ReLu activation functions originally, we elaborate that the differentiability requirement practically affects our choice of activation functions.
Time derivative Without loss of generality, we assume we optimize a function g that takes the output of the decoder &pdec Usys: t] (defined in Appendix A.2 as its input The derivative of g w.r.t is:
dg (Oyden (t), ueys: t))tt) V
Oupdee (oonee
d dec
0g
(ODE>!((), Usyas 8] + ai(!)
+
(11)
dl
01
01
where ΦODE. represents the neural network structure modeling the latent ODEs. is clear that the differentiability of model output w.r.t time f needs to be enforeed through 8t
Oonaes as long as the time t is
taken into account in the decoder part.
Initial Condition Derivative The derivative of the model output w.r.t the initial condition æo can be reprosented as:
dg (ousen [((t), Usys: t)
dxx(ta)
b ude at)
al(I)
dopdec
Ous
Vopacod
(x((t) +
Ouays
@x((o)
a (Oopayo (haye) + diag (0o (go (heys)s))
(12)
Doudec Ol(t)
@xo(t) +.)
@l((a)
dopdec Orssys
V Opaco
+
Oxr(t)
Thays
orsys
dr sye
where Tays and hays are defined in Appendix A.2 0n sys (hays) + diag sys (haya))) represent the reparame- terization based sampling of the latent variable Uays. where is a random sample from Gaussian prior. Note that the derivative of the state value w.r.t to initial condition at(1)
Baitoo) can be obtained through the adjoint method, derivation has been provided in Appendix C of Chen et al. (2018), without full specification of such cumbersome terms, is sufficient to indicate the sensibility of having differentiable activation
well the decoder functions in context encoding (Br.ys) context latent representations Dr.y.
(theys
(B (0ensen) For practical calculation of the gradient of the model w.r.t decision variables both training al(1) and optimization, we simply leverage the standard backpropagation through the ODE solver (also known the 'diseretiza-and-pptinni (e.g. Onken & Ruthotto (2020); Kidger 2022 to obtain the gradient due to its speed advantage and minimum implementation effort.
In Figure i we also empirically show the necessity of these differentiable choices by inspecting the time and initial condition derivative with respect to g(.).
we note that this is sufficient show that the requirements of time differentiabllity irrespective the consideration of do dea ((t), Mays.t)
19

---

Published in Transactions on Machine Learning Research (5/2025)
04
30
62
-10
-44
-35
-15
-10
(a) (a(t)
using ReLU.
(6) using ReLU.
(c)
Oe(O) 81
using SiLU.
Pn((
using SiLU,
Figure 6: Initial condition gradient and time derivative w.r.t different activation function choices, the comparison is conducted with the same model structure but different activation functions.
B ELBO Derivation and Training Details
B.1 SANODEP ELBO derivation
Note that since the forecasting scenario is a special case of interpolating scenario, below we derive through a unified representation:
log (X IT UT T Thew)
= llog ( (X hoodd Teew Usys: Dj p (Weys T U U p (Lnnew) d'utsys Donew
I
(1og) I plusyalTe Trew) p (Ldmew) I (uays T U Th UT UTT new)
new
(Maya T U Thew Treww
(D (t) )now!)
p (XT new Teow Usyer Linee A'uays Lonew 9 (inew (6
(13)
(uayal IT T E Busss TTuT UTew Thewlo (D (f log @ (Usys TC Thew U UTT TTow)
p (L
X
T new news Usys: LJ Idncw))
(L (1g a
E
I "(uuaall (etuuyllEETAn (T UTT (eewl(( 1e( Lynow (f log (XT new nnw Usysr D.an))
-KT (Msyss IT U The UTT 1|q (usys / U Tre
KI, [g (Ig (16 I (new)) 1|p (L)
We note the variational posterior q(u T encoder in an amortized approach similar in Garnelo et. al. (2018b).
[LT
has been obtained through the (+C xf ()..))
B.2 Learning and Prediction Process of SANOEP
The practical learning process of SANODEP is summarized in Algorithm. [ with parameters listed in section. F aas well. In all of our subsequent experiments, use Mm 0, M = 10, Nz, = 100, Noys 20, N Ngrid 100, Mmin = 1, Mmaz 10, Umin = 0, Pmaz = 45. We note that all the loops have been practically implemented in batch, and all the subsampling approach has been implemented practically through masking in Jax.
Challenges to incorporate attention Finally, we remark that meta-learning of ODE systems is computa- tionally challenging with models having an attention block. With lines 15-16 of Algorithm. D the training model takes a batch dataset of Nays Dimn [Tj] points simultaneously, while the data are subsampled context and target data separately. This is handled with mask. Contemporary deep learning libraries (e.g., Jax) cannot leverage sparse masks to reduce attention computation complexity inherently. Hence, it is
20

---

Published in Transactions on Machine Learning Research (5/2025)
infeasible to calculate self or cross-attention over keys with the number Noo [Tyl. which in our cases is tens of thousands. Generally, it is well recognized (e.g.. Feng et al. (2023)) that utilizing attention blocks in neural processes with large data is computationally intensive. Additionally, lines 6 of Algorithm. proferably require a dynamic update of C with the initial condition, which is also computationally intensive for self' attention.
C Time Delay Constraint Optimization
C.1 Acquisition Functions
We advocate that the framework is agnostic the acquisition function forms long as it can be written expected utility form and satisfy the property mentioned in Lemma. C.1 however mention the acquisition function form we specifically utilized for performing few-shot BO in Scction 6.2 as:
NMc
a (Io. f1...., iN. Po (x C, T" mo)) ~
Σ
HVI
(Oudeol
(TT)., Uaysy T TITT""
(14)
NMC
In practice, we infer the current Pareto frontier F² from SANODEP as well from the M context trajectories TC, We use NMc = 32 as the MC sample size to approximate the expected hypervolume improvement in all our experiments.
Acquisition Function Optimization We note the acquisition optimizer both initial condition (Eq. 9 as well as optimal time scheduling (Eq. (10) involves constraint optimization. We utilize trust region- based constraint optimization available in Scipy to optimize acquisition function in both optimization problems starting with uniform time scheduling (i.e., (inspace(l,, t1, N)) Note that for the initial stage optimization, we additionally leverage (10 instances) multi-start procedure on the initial condition decision variable only to boost the optimization performance.
C.2 Search Space Reduction
Proof of the search space reduction We prove that the search space can be reduced, without eliminating the global maximum in acquisition function optimization processes (Eq. (9))
Lemma C.1. For batch acquisition function defined as expected utility En( where the utility function u() is monotonic w.r.t set inclusion (S
u(T) 2 u(S)), the aquisition function is also
monotonic w.r.l set inclusion: ((T) 2 a(S).
Proof. Let T S U T/S, then a(T) aSS) = Ensss) ERCSS) EP()) BP(//S99) (n(T) ((S) 2 0, hence prove complete.
Corollary C.1.1. Assume that the ubility function u(-) is monotonic with respect to the inclusion of the set in each set seperately: Let X = X2..........
and Y = YY,,, be two collections of sets, if Y Vi E {l...,n] n} then x < w(Y). Then, the batch acquisition function defined as = Ex(() (u())) satisfies a(Y) 2 ((X).
Lemma C.2. QEHVI is monotonic acquisition function with respect set inclusion
Proof. Using Corollary, C.1.1 Lemma C.2 and the well-acknowledged fact that the hypervolume indicator is monotonic with respect to set inclusion on each output dimensionality, the proof is straightforward. We note that the above property is asymptotically hold when qEHVI is monte carlo approximated.
Finally, we start the proof of Theorem.
Proof. when N E point in between, forming an augmented set T as input, with Lemma C.2 and Corollary. C.1.I the proof hence complete.
1], 3i [1, qu 1] s.t. 4 f:-1 > 2At meaning that can insert new batch
httpp:///occs.ee http:////.ccccceooo.oo//////e/tn
21

---

Published in Transactions on Machine Learning Research (5/2025)
We finally remark that, with Lemma. C11 a large series of batch improvement- based acquisition functions (e.g,, parallel Expected Improvement, parallel Probability of Improvement, and parallel information-theeoreiee acquisition functions that can be written as expected utility form and also preserve the set inclusion monotonicity property (e.g., Qing et al. (2023)) )) also embrace the same search space reduction benefit.
D Experimental Details and Additional Results
D.1 Meta Training Data Definition
Lotka- Voterra
2D Cases
da:
Ta - B
(15)
dt
6z122 - y22
The initial condition is sample from a=lo) U(0.1,3) the dynamics F sampled from U( 1), ß ~ U (1,,)) 2), d U(0.5,1.5), 7 (0.5, 1.5), tmax = 15. We run the meta-learn this system distribution with 300 epochs.
3D cases
031 - -3 43
dx
802211 - -132 - Cx2T3
(16)
dt
1/8332 0r3
The initial condition sampled from æ(to) (0.1, The dynamics F are sampled from a U (3. 1).
8 2(1 5 ~ U(0.5, 1.5), y (00.,, U(0.5, 1.5), Ç (0..5, " U 1.5), and
0 ~~ U(0,5,155) with 4mar 15. We run the meta-learning algorithm on this system distribution for 300 epochs.
Brussclator Prigogine & Lefever 1968
da:
A + 25 (B + 1)m1 Br
(17)
dl
22
we leverage A ~ U), B ~ 3),,,3)) æ ~ ((0.1,200) Uman 15. We run the meta-learning algorithm on this system distribution for 300 epochs.
SIR model
dS = BSI,
dI
= 3SI - II, dR
= =1I
(18)
dt
dt
dt
The initial condition sampled from S(to) U(1.0,3.00 [(to) 0.01, R(to) 0, and 8 UO0.1, y U(0.1, 10). kmax 1. We run the meta-learning algorithm this system distribution for 300 epochs.
SIRD model
dS --SS,,
dl = BST 7 nI, dR dD dt
-1,
=
(19)
dt
dt
dt
where S(to) U(10,30) T(to) = 0.01, R(to) (tt) 0. tmax 1, 8 10.5,2.0] 7 [0.1, 10), and / 5001, We run the meta-learn on this system distribution with 300 epochs.
GP Vector Field
dx f(x)
(20)
dt
We assume vector-valued GP as the non-parametric prior for the vector field f GP x(,, with correlations between each output as the most generic cases. To sample the vector field, we leverage the
22

---

Published in Transactions on Machine Learning Research (5/2025)
parametric approximation of GPs through the random Fourier feature approximation of the RBF kernel Rahimi & Recht (007) as a Bayesian linear model, which common approach (e.g. Qing et al. (2022) to obtain diffcrentiable GP samples, the vector fields practice are generated from kernel lengthscale 0.8 and signal variance 1. The implementation also utilizes the paramelrie sampling approach of the GPJax (Pinder
& Dodd 2022) library.
D.2 Additional Experimental Results
We report all model's performance interpolating tasks and forecasting tasks, averaged over known trajectory ranges M, in Table. and Table. 3
Table 2: Model comparison (mean-squarcd-crror x10-2) for range of dynamical systems.
Models
Lotka- Volterra (2d)
Brusselttor (2d)
Selkov (2d)
(3d)
Lotku- Volterra (3d)
SIRD (4d)
Forecasting
GP
87.2 31.9
39.3 = 61.0
31.8 +: 6.72
1769.33
990.6
41.4 :: 15.8
1322.845 489.4
NODEP-A
102.9 6.85
105.2 10.8
16.1 1.14
6367. 6367.784 72..88444 161.8
64.8 2.56
3368.685 189.5
NP-A
3.1
21.6
1.6 0.0916
952.3
32.8 0.85
721.2 108.8
NP-1
47.3
1.7 0.0882
1037.981 207.2
34.8
88.3
SANODEP-A=
48.4
20.1 0.816
0.0376
+ 89.0 873.1
36.2 1.17
426.4 44.0
SANODEP-A 0.5
50.6
18.8 0.467
1.28 1.28==00m
35.5 0.824
410.5 27.3
Interpolating
GP
52.2 + 20.2
9.62 19.0
11.7
173.7 296.8
13.7 26.8
280.7 137.3
NODEP-İ
38.9
7.12 0.559
367.5
=P-
33.2 2.59
0..4
0.742
286.2
24.7
302.0 110.9
NP-À 0.5
33.2
1.407
0.807 + 0.0638
25.9 1.39
2446
SANODEP =
32.1
0.556
8.74 1.09
0.53 0.03
257.5 60.3
25.2 1.84
139.5 + 13.5
SANODEP-A 0.5
36.6 + 7.73
8.06 + 0.265
0.573 =: 0.0247
307.8 33.5
25.2 1 0.464
143.8 = 14.0
Table 3: Model comparison (negative log-likelihood x10²) for a range of dynamical systems.
Modela
Volka-Volterra (2d)
Brusselator (2d)
Seikov (2d)
SIR (3d)
Locke- Locka-Voltest Volterra (3d)
SIRD (ld)
Forecasting
18.2::77 7.71
30.2 : 25.7
17.1
11143.352 31169.575
19.3 5.75
27321.033
57173.003
NODEP-A=
119.3 21.7
79.1 13.9
136.4 32.1
1066 331.3
625.5 174.2
447.7 86.0
NP-À - 0.0
98.9 52.2
17.7 1..7::40 4.41
4.28
4: 1.34
162.6 1: 78.6
491.2 122.1
197.5 151.7
NP-A = 0.5
200.9 :: 58.0
37.7
18.5
5.93 0.836
279.7 198.2
961.1 454.1
75.7 12.8
SANODEP-A=D
22.6 2.99
206.9 79.4
277.5
17.4
SANODEP-A 0.5
27.0 + 13.0
23.7 + 4.55
8.55 + 1.42
219.4 + 139.7
329.0
103.8
48.3 11.0
Interpolating
GP
3.27 1 1.24
1.58 2.0
5.27 1: 1.87
718.0 1: 2426.002
2.69 1.61
1045. 054..453-4 WW17 646
NODEP-1=0
32.6 6.73
2.88
66.1
NP-A = 0.0
18.7 7.46
0.254 0.61
-2.99
0.215
22.4 10.1
03.8 11.2
30.1 23.1
NP- - 0.5
34.8 14.0
2.88 2.01
-1.97 0.147
40.1 28.4
122.1 49.6
12.5 6.12
SANODEP-A=
-1.74
41.321.7
SANODEP-A 0.5
6.08 + 3.63
1.13 0.596
-1.5+00284
31.3+ 19.3
47.3 14.6
7.4 1.97
D.3 Model Prediction Visual Comparison
We provide the prediction comparison of different models on a specific realization of ODE systems in Fig.
13
D.4
Optimization Problem Definition
The definitions of the optimization problem provided in Table.
23

---

Published in Transactions on Machine Learning Research (5/2025)
-Vaterra (2d)
Brusselator (2d)
(2d)
SIR (3d)
Voterra
SIRD (4d)
10a
10a
10
10*
10a
10
10a
10
10a
10-2
10-1
10
10a
10
100
ISEE
10
10
10
104
104
10
10-3
10-1
10
10"
10
Number of Context Trajectories
MANODDE-VK--A-0X-
NP-A- D
NP-A-02 0.5
NODEP
SANODEP.A D.D D.D
SANODEP A 0.5
GP
Figure 7: Sequential model evaluation performance including GPs. The first row corresponds to forecasting performance (prediction with only the initial condition known), and the second row represents an interpolating
setting.
Lotka-Votera (2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
Lotka-Voterra (3d)
SIRD (4d)
10
10a
10a
10'
107
10
10a
10a
10a
10
10
10
10a
10
10'
10a
10
-10
10
-10-
-10-
-10
-10
10
10
-102
-10
10
10
10
10
10
10
10a
10a
10a
107
10a
10
10a
10a
10a
10a
102
10
10a
10a
10a
10a
10
10
-10
10
-10*
-10
-10
-10
18
-10a
-10*
10a
Number of Context Trajectories
SANODEP-Ix--A 0.5
NP-A
NP-A 0.5
NODEP
SANODEP-A 0.0
SANODEP-A 0.5
GP
Figure 8: Sequential model evaluation performance in terms of Soggttive-oo--Likeihooe (NLL) including GPs. The first row corresponds to forecasting performance (prediction with only the initial condition known), and the second row represents an interpolating setting.
Table 4: Optimization problem formulations
Problem
Lotks- Volserta (2d)
Brusselator (2x)
Selkov (2:)
SIR moddel (3d)
Lotka- Volterra (3d)
SIRD (4d)
formulation
#2
>>. 0.05
m// >> 0.05 >
r = 10.
r 10.
7 [0, 10/,
[0, 1].
7 e J0, 15/, F [0.0) 20
T [0.
Design space
2.01
0..1
[10]
[10
ato andes
Zde
ape (itoo, 0.1 Zdees 0j
20 aiice
ato [andes
D.1
utdee 0.0
1.5
8 0.5
0.5.
Optimization
A
ODE definition
8 1.2,
A = 0.8,
B=
0.25, 6=
10= 1.5
7 5
1.5
0.5 P
5 1.5
( 1.2
1/
8 1.5
Reference point
|-777, 12.686
--1.66,3.8877
[-.117..5.110)
[0.51151, 0.79,.677616
[-177557, 13.1667]
(0.52198, 1.01]
D.5 Inference Time
We provide the parameter counts, training time and prediction time comparison of SANODEP vs. baselines in Table. 5 note that all the models have the same latent dimensionality for fair performance comparison, which is regarded as the key bottleneck of the latent variable typed meta-learning model (Kim et al. 2019
24

---

Published in Transactions on Machine Learning Research (5/2025)
Lotka- Makka-Voeera (2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
Lotka-Voterra (3d)
SIRD (4d)
7.5
20
2
2.5
10
1
0-
0
0.0
0.5
1.0
15
0.0
0.5
1.0
1.5 0.0
0.5
1.0 0.0
0.5
1.0
0.0
0.5
1.0
1.5
0.0
0.5
1.0
3
30
4
20
2
L0
2
0.5
0
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5 0.0
1.0 0.0
0.5
0.0
1.0
0.5
State 1 Prediction
State Prediction
Tie
Real State 2 Value
Real State 4 Value
State 2 Prediction
State Prediction
Real State
Value
Real State 3 Value
Figure Neural Processos (A 0.5) podel merformance on test system in different meta-lcarning ODE problems.
Vicka-Voterra (2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
Lotka-Voterra (3d)
SIRD (4d)
30
10
20
20
5
10
0
0
0
1.5
1.5
2.0
30
55
2
20
4
20
10
2
0.5
0
0.0
0
2
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
1.5
0.0
0.5
1.0
State Prediction
State Prediction
Then
Real State 2 Value
Real State Value
State 2 Prediction
State Prediction
Real State Value
Real State 3 Value
Figure 10: SANODEP-(A = 0.5) model performance test system in different meta-learning ODE problems.
The parameter increase of SANODEP compared to NODEP the additional encoder o, to handle the augmented state variable.
Table 5: Comparison of the number parameters and prediction time of different models. We measure run time on the Lotka- Voterra (d 2) problem using an NVIDIA A40 GPU, We note training time and prediction time per mini-batch, the optimization query time is per querying one observations.
Model
Number of Parameters
Training Time (a)
Prediction Time (a)
Optimization Query Time (s)
GP NP
NA
NA 0.3016 + 0.0791
NA
23.81 43.94
21814
0.2551 + 0.2337
18.01 + 18.36
NODEP
30286
0.3792 + 0.1892
0,3421 0.2412
34.71 = 20.62
SANODEP
35514
0.3813 + 0.2423
0.3482 + 0.2731
39.47 = 37.30
25

---

Published in Transactions on Machine Learning Research (5/2025)
Notka-Voterra (2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
Lotka-Voterra (3d)
SIRD (4d)
130
15
30
40
5
10
20
0
10
0
0
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0 0.0
0.5
1.0 0.0
0.5
1.0
1.5 0.0
0.5
1.0
3
30
10
10
20
2
5
10
10
0.5
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
1.0
0.0
0.5
0.5
State 1 Prediction
State Prediction
Tien
Real State 2 Value
Real State 4 Value
State 2 Prediction
State Prediction
Real State Value
Real State 3 Value
Figure 11: PISANNDEP-(A 0.5) model performance on test system in different mota-learning ODE problems.
Lotka-Voterra (2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
Lotka-Voterra (3d)
SIRD (4d)
30
2.
2
10
0.5
1
0
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5 0.0
0.5
1.0 0.0
0.5
1.0
0.0
0.5
1.0
1.5
0.0
0.5
1.0
3
30
4
20
2
20
2
10
0.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
1.0 0.0
0.5
State 1 Prediction
State Prediction
Tex
Real State 2 Value
Real State 4 Value
State 2 Prediction
State Prediction
Real State Value
Real State 3 Value
Figure 12: NODEP (A = 0) model performance on the test system different meta-learning ODE problems.
E Investigation of the prior strength, F
E.1 Strong Prior information: Physics-Informed SANODEP
When f is realization of F and we know the parametric form of the kinetic model, easily integrate this kinetic form in SANODEP's model structure. We call this model variant Physical Informed (PI)-SANODEP.
As a concrete example, consider f as realization of the Lotka- Volterra (LV) problem with unknown
parameters. Under the SANODEP framework, we can replace fan in Eq. 3 with the following physics- informed form:
x1
(21)
dt x2
8g 322 22
26

---

Published in Transactions on Machine Learning Research (5/2025)
(2d)
Brusselator (2d)
Selkov (2d)
SIR (3d)
Lotka-Voterra (3d)
SIRD (4d)
30
20
2
10
0
-2
-2
0
-4
4
-2
10
0.0
0.5
1.0
15
0.0
0.5
1.0
1.5
0.0
0.5
1.0 0.0
0.5
10 0.0
0.5
1.0
1.5 0.0
0.5
1.0
30
20
2
20
10
0
-2
0.0
0.5
1.0
0.0
0.5
0.5
Preducted State Value Predicted State Value
Predicted State Value
Real State Value
Real State Value
Predicted State Value
Real State Viatue
Real State Value
Figure 13: Performance of the GP model on the test system in different meta-lcarning ODE problems.
In this form, Usys = [av, S9, 80, yol. and I(ta) = 2o and =(t)= x(t are strictly positive. We log- log-normal distributions for the conditional prior and variational posterior for Uays in the encoder. For PI-SANODEP's training loss function L' since these parameters are known during meta-training, we augment our loss (Eq. by incorporating likelihood term for these parameters, leading to the following expression:
Since the parameters of the LV system
SANODEP
PI-SANODEP
α
6
5
3
Density
0
0.0
0.5
1.0
1.5
0.0
1.0
1.5
State
Prediction
State 2 Prediction
Ti
Real State
Value
Real State 2 Value
Figure 14: Model Comparison of SANODEP Physics Informed (PI) SANODEP, PI-SANODEP provides remarkable better prediction accuracy. In addition, PI-SANODEP provides a reasonable estimation (curve) on system parameters (dashed line).
log p (XT nows Usys 1 UTC
Teem: T TIT log p (XT new Usys: T TTE
Temm: T
log p (usya MTsyyITC. T T T
(22)
Since Usya is known during meta-training, network to calculate its likelihood p (Usya TTUTE new: new) = = (usysi Sp (hsys) diag (oo aye (hayss)
the encoder (described in Appendix A.2 inference
Regarding ( Mays; T U TE Tif T we however omit the dependence of usys and still make of Eq. (8) to calculate. Thus, the final loss function for PI-SANODEP extends the original SANODEP loss with an added likelihood term log (Msys|77 |T UTC estimation.
new: Th capturing both the dynamics and parameter
10tn PI-SANODEP, since the ODE is explicitly defined the kinetic form (e.g., Eq. 21 the decoder is only used to predict variance.
27

---

Published in Transactions on Machine Learning Research (5/2025)
Brusselator (2d) SANODEP
Brusselator (2d) PI-SANODEP
A
2
Value
1
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
Selkov (2d)
Selkov (2d)
SANODEP
PI-SANODEP
a
Value
2
SIR (3d)
SIR (3d)
SANODEP
PI-SANODEP
Value
20
200
0.0
0.5
1.0
0.0
0.5
1.0
10
Lotka-Voterra (3d) Lotka-Voteraa (3d)
SANODEP
PI-SANODEP
α
E
5
Value
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.5
1.0
1.0
1.5
2.0
0.5
1.0
0.5
1.0
1.5
C
n
0.5
1.5
0.5
1.0
1.5
0.5
1.0
1.5
0.5
1.0
1.5
SIRD (4d)
SIRD (4d)
SANODEP
PI-SANODEP
B
H
Value
20
0
0.0
0.5
1.0
0.0
0.5
1.0
1
2
10
2
Prediction
Real State 1 Value
State Prediction
State Prediction
State 2 Prediction
Real State 2 Value
Posterior
Real State 3 Value
Real State 4 Value
Fiew
Figure 15: Additional experiments for meta-parameter estimation based on PI-SANODEP, the shaded area represents the meta-training parameter support S supp(P). It be that PI-SANODEP not only improves upon SANODEP over prediction (e.g., in the SIR example) but also provides a reasonable parameter estimation of the test system. Parameter estimation tends to become difficult when the number of parameters is high as can been seen from the 3-dimensional Lotka-Votee Voterra system.
We demonstrate PI-SANODEPss performance with SANODEP on Potka-Votersy problem in Figure 14 and the rest of the problem in Figure 15 we note PI-SANODEP also demonstrates the potential of concurrent few-shot prediction together with parameter estimation of dynamical systems. Compared with Bayesian inference- based parameter estimation strategies (e.g. Wu & Lysy (0244), ), which may take minutes to perform
28

---

Published in Transactions on Machine Learning Research (5/2025)
estimation, PI-SANODEP shares the same advantage as simulation-based inference Cranmer et al. 2020 and can estimate almost instantly since it only requires a forward pass of the encoder.
E.2 Weaker Prior of F: Cross-domain generalization possibility
What if we know nothing about F? From meta-learning perspective, to still frame the problem as in- distribution evaluation, this demands an extremely flexible prior for F to be utilizable in the target dynamical system. To our best of knowledge, meta-learning dynamical systems to be leveraged cross-domain scenarios is still an open direction. Certain approaches have resorted to incorporating additional dynamical system properties (e.g., energy conserving system Song & Jeong (2023)) which limit the scope of system types. We propose and leverage vector valued GP prior for F flexible task distribution.
Sample (X)
Sample f(x)
Sample fixh
Sample fix);
Sample f(x)
Sample fixh
Sample f(x22
Sample
Sample fix12
Sample fixl2
Sample f(x)2
Sample f(x);
Vector Field
Vector Field
Vector Field
Vector Field
Vector Field
Vector Field
Sample Trajectory
Sample Trajectory
Sample Trajectory
Sample Trajectory
Sample Trajectory
Sample Trajectory
State Valuue
Sp0e vilua
State Value
State Valuue
Spate Value
S00 VOL
Time NP prediction on
Time NP prediction
Time
Time NP predi ctlon
Time NP prediction
Time
NP prediction
NP prodiction
7.75
2.5
0-
0.0
0.0
5
2.5
SANODEP prediction
SANODEP prediction
SANNODEPrrediction
SANODEP prediction
SANODEP prediction
SANODEP prediction
SANDDEP prediction
1.0
0.0
-5
-2
2.5
0.5 Time
1.0
0.0
0.5 Time
1.0
0.0
0.5 1.0
0.0
0.5
1.0
0.0
0.5
1.0 0.0
0.5
1.0
Time
Tima
Time
Time
Figure 16: Task distribution generated by vector-valued Gaussian Process (GP) priors. The first two display samples from parametric GP priors. The third and fourth rows show the sampled vector fields and the resulting trajectories starting from sampled initial conditions. The last two rows depict the meta-model predictions on these trajectories. It can be seen that while the GP-based vector field has introduced a very flexible task distribution, both meta-learn models tend to underfit on these trajectories.
Inspired by the recent approach of using Gaussian Processes (GPs) to model vector fields (e.g., Heinonen et al. (2018),, we propose use vector-valued GP the non-parametric prior for the vector field, denoted by (O, K(x, (()x)))), as task distributions for meta-learning for dynamical systems. Further details on the data used for meta-training are provided in Appendix D.I The dynamical systems sampled from this
29

---

Published in Transactions on Machine Learning Research (5/2025)
SANODEP Brusselator (2d)
Lotka-Voterra (2d)
Brusselator (2d)
Selkov (2d)
Lotka-Voteria (2d)
Sel (2d)
5.0
5.0
2.5
2.5
2.5
-2
2-
2.5
-2
0.0
0.5
1.0
1.3
0.0
0.
1.0
1.3
0.0
0.5
1.0
1.3
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
7.5
7.5
5.0
5.0
1
0.0
0.0
-2.5
0
0
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
State Prediction
State Prediction
Real State
Value
Real State Value
Figure 17: Cross domain generalization exploration: SANODEP and NP meta trained on GP-ODE seemss to able to capture the trajectorics from Bruscclator system.
vector-valued GP prior, as illustrated Figure T6 demonstrate ability represent a broad distribution of dynamical systems. Additionally, we offer visual comparison of SANODEP and NP on the sampled trajectories. This comparison reveals tradc-off between fitting accuracy and trajectory ficxibility for both models. Specifically, both models tend to significantly underfit trajectories that involve oscillations, indicating the noticable complicity of this meta-learning problem.
We also explore whether this flexible task distribution facilitates cross-domain generalization. To this end, we employ the trained model on different dynamical systems that share the same state variable dimensionality, as shown Figure 17 Both meta-learned models demonstrato some potential for cross-domain generalization particularly from the Brusselator system. However, generalization systems involving oscillatory behaviors proves to be especially challenging.
30

---

Published in Transactions on Machine Learning Research (5/2025)
F
Notations
Table 6: Nomenclature Table
Notation
Meaning
State variable vector Dimension of the state space Initial time for the dynamical system Trajectory time samples, sampled irregularly Termination time for a trajectory Context Sets
dæ
lo
t
t
f
Vector field governing the dynamics Aggregation function used in optimization Latent initial condition
g(-)
Lo
Ln
Dynamics of the system Context encoding
or
r
Context representation vector
lo
Realizations of latent initial conditions
I(t)
Latent states time t Realizations of the control term representing latent dynamics Neural network parameterized vector field Parameterisation of the ODE system Numerical solver implementation of the vector field Instantenous time variable in continuous dynamics Dimension of the latent state space Target Sets
ul uysss))
fnn
Oode
f Fevalue
d
TT
Target times for predictions Target state values
XT
F
Distribution of dynamical systems
P
Parametric Task distribution, often stochastic State variable space Initial condition space Time domain for the system evolution
X
X
M
Number context. trajectories
N
Number of context elements in trajectory Context observations from M trajectories Context observations from the Ith trajectory A new trajectory which augments the context set C Context observations from the new trajectory Target observations from the new trajectory Latent variable conditioned on + 1 trajectories Realization of system-aware latent dynamics System context represenlation vector Context encoding augmented with the initial state Target state values for the new trajectory Targel times for the new trajectory State vector samples from the new trajectory Time samples from the new trajectory Number observations in the new trajectory
T
T
Trew
new
T
D Days
ssys
Tsys
d sys
x
new
Th
Knews
knew;
K
Lo
SANODEP's multi scenario loss function parametrised by 0 SANODEP prediction parametrised by
pa
forecast
Bernoulli Indicator
Parameter for Bernoulli distribution in decision-making process
31

---

Published in Transactions on Machine Learning Research (5/2025)
Table 6: Nomenclature Table (Contimued)
Notation
Mcaning
Nopt
Search space of the number of observations in the trajectory
At
Minimum delay between observations
tmaz
Maximum time observation
Nmaz
Maximum number of observations per trajectory
a(-)
Batch acquisition function
HVI
Hypervolume improvement based on the Pareto frontier Pareto frontier
F*
Mmin
Minimum number of context trajectories observed in system Maximum number of context trajectories observed in system Sampled number of initial conditions within system
Maa
Nao
N Nays
Number of dynamical system samples per training iteralion Number points in the time grid for trajectory evaluation
Ngrid
Mmann Mona
Min and max context points within trajectory
Bmiin Amaa
Min and max targel points beyond context in a trajectory
S
The support set of the stochastic parameters for meta-task distribution
Torid
Time grid used for trajectory evaluations
32


Reviews:
[{"review_id": "CJ7fdBDChD", "reviewer": "Reviewer_2ir2", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n\n- It tackles an interesting application of Bayesian optimization to ODE systems.\n\nWeaknesses\n\n- Technical contributions are somewhat limited.\n- Baseline methods are outdated.\n- Intuition behind the components included in the proposed algorithm is missing.\n\nQuestions and Suggestions\n\n- What is the meaning of a curly parenthesis in Equation (1)?\n- Why is it a multi-objective optimization problem?  Why do you want to minimize evolution termination time?\n- In Equation (2), maximize can be just $\\max$.\n- For the sentence \"Simply put, we wish to maximize objective g early,\" please check if it is grammatically correct.\n- Every equation is considered as a sentence, so that it should end with comma (,) or period (.). For example, Equations (1), (2), and (3) don't end with either , or ..\n- Could you describe the difference between Equations (5) and (8) in plain text?  I believe that this should be added to improve the presentation.\n- I think \"section. 6.1\" and \"Appendix. A.3\" should be \"Section 6.1\" and \"Appendix A.3.\"  Specifically, there should be no period.  Also, Fig. is odd to me.  It can be just Figure.  There are too many similar issues; please fix all.\n- In Section 3.1, why do the authors use DeepSets?  There exist more sophisticated set-taking methods like methods with different pooling methods.\n- In Page 6, \"can obtained\" should be \"can obtain.\"\n- Basically, I don't understand Section 4.1.  Why do we need to find an initial condition using Equation (9)?  Isn't it given by practitioners?\n- In Page 7, \"Eq. 9,10\" should be Equations (9) and (10).\n- In Section 4.2, do the authors utilize the vanilla qEHVI?  Which implementation did you use?\n- The title of Section 5, \"Related Works\" should be \"Related Work.\"  Work is uncountable in this usage.\n- Why do the authors use GP with the Matern 3/2 kernel?  Which implementation did you use for GP?\n- In Page 9, \"Table 2 and Table 3\" should be \"Tables 2 and 3.\"\n- Font size of the caption of Table 1 is too small.\n- I think NP and vanilla GP with the Matern 3/2 kernel are somewhat outdated.  The authors should try to compare more recent methods, for example, TNP and more sophisticated BO algorithms, to the proposed method.\n- In Page 10, some paragraph heads have colon.  I believe they are not consistent.\n- In Page 10, question marks for some paragraph heads are missing.\n- In Page 10, \"demonstrate\" should be \"demonstrates.\"\n- I think that Algorithms 1 and 2 can be located in the main article, not in the appendix.  They are important to understand the proposed algorithm.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "- This study addresses optimizing initial conditions and termination time in dynamical systems governed by unknown ODEs, with costly evaluations and delayed state measurements.\n- Few-shot Bayesian optimization using prior information is introduced to efficiently identify optimal conditions in limited trials.\n- System-aware neural ODE processes (SANODEP) is developed, expanding neural ODE processes to meta-learn ODE systems via a novel context embedding block.\n- A two-stage Bayesian optimization framework is proposed to incorporate search space constraints, optimizing initial conditions and observation timings effectively.\n- Experiments demonstrate SANODEP's potential in few-shot Bayesian optimization, while exploring its adaptability to varying levels of prior information.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please see the textboxes above.", "broader_impact_concerns": "I don't think that this work has broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Ue83x2HPS8", "reviewer": "Reviewer_XihC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- The paper effectively addresses the practical challenge of optimizing dynamical systems with limited evaluative resources.\n- The proposed SANODEP method is innovative in its integration of system-awareness into the Neural ODE Processes.\n- The experimental section is thorough, with a variety of dynamical systems tested, demonstrating the robustness of SANODEP across different scenarios.\n\nWeaknesses\n- How does the SANODEP framework handle dynamical systems with highly nonlinear or complex interactions that may not be well-captured by the assumed ODE forms?\n- How sensitive is SANODEP to the choice of hyperparameters, particularly those related to the context embedding block and the bi-scenario loss function?\n- What are the theoretical guarantees for the convergence of the proposed two-stage Bayesian Optimization framework, especially when dealing with non-convex objective functions common in dynamical systems?\n- What is the robustness of SANODEP against adversarial attacks or perturbations in the input data？\n- There are several relevant methods missing in the manuscript that need to be discussed and analyzed: 'PGODE: Towards High-quality System Dynamics Modeling' and 'EGODE: An Event-attended Graph ODE Framework for Modeling Rigid Dynamics'.\n- Can the authors discuss the theoretical implications of the proposed method on the exploration-exploitation trade-off in the context of Bayesian Optimization, and how this might influence the optimization of initial conditions and observation timings?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents a meta-learning approach for few-shot bayesian optimization in dynamical systems, which addresses the challenges of optimizing initial conditions and termination times by leveraging prior system information.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See Weaknesses.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "hHTRmVRhwg", "reviewer": "Reviewer_jc5T", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The problem of finding the optimal initial conditions of unkown dynamical system is relevant.\n- The paper has a detailed related work section. \n\nWeakness:\n- The contribution is not clearly defined, as the change to NODEP framework seems minor. See Requested Changes in order to make the contribution clearer.\n- Fig 1. is hard to grasp and does not show the benefit of the model: difference filled vs blank red dots? GP seems to do as good?\n- The notation is partially confusing (See \"Minor\" und requested changes)\n- While the paper is overall well written, Sec. 4 is hard to follow. \n- (PI)-SANODEP seems orthogonal to the rest of the papers story.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes an extension of the NODEP framework towards dealing with multiple trajectories. The multiple trajectories are constructed by realization of different initial conditions. \nThe novel system is then used as proxy model in a Bayesian Optimization problem.  The authors frame the optimization problem as a two step problem, where first the initial condition is optimized and second the termination time is optimized. \nSimilar optimization problems may be found in real world production systems.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Clarifying difference to NODEP: \n-  Multiple Trajectories:  Is it possible deal with different initial conditions by increasing the feature dimension of $x$, i.e, replacing $\\phi=(t_i, x_i)$ with $\\phi=(t_i, z_i)$ with $z_i = [x_0, x_i]$ and treating it as one dataset? Similar as augmented neural odes are constructed. \n-  Loss: As the loss function is different in comparison to the original NODEP paper: Can you discuss the influence and emperically demonstrate the influence of each term in Eq. 6?\n\nGP Baseline:\n- Fig 4. suggests that the proposed framework is able to outperform the GP baseline. However, the authors provide little information about how the GP was constructed. Beyond that, the authors omit comparing against GP-ODEs, which was though cited.\nSee for example Heinonen 2018 and subsequent work for this. I think a comparison against this model class is crucial.\n\nMinors (Notation):\n- Section 2.2 --> notation switching $d$ vs $u$.\n- what is the benefit of having $f_{evolve}$ if you can just write $l(t)$?", "broader_impact_concerns": "Not applicable.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This submission considers the problem of optimizing initial conditions and termination time in dynamical systems governed by unknown ordinary differential equations (ODEs). It takes the ODE approach, extends the NODEP method, proposes a few-shot Bayesian optimization  framework to incorporate system's prior information.

Some reviewers consider it to address a relevant and practical challenge in optimizing dynamical systems with limited evaluative resources. The extension of NODEP to the meta-learning setting with system awareness improves learning efficiency and shows some level of novelty. However, reviewers consider the technical contribution to be limited and the scope of the problem is relatively narrow. The empirical comparison to the limited set of GP and NP-based meta-learning methods remains a concern to reviewers even though the authors explain why other more recent models and standalone approaches do not apply or require a significant effort to be adapted to the problem under study.

The authors provide a revision to address formatting issues and provide more clarification following the reviewer discussion. While some reviewers are still concerned of its limited scope and baselines, they all agree that the proposed method has demonstrated its value in addressing the specific optimization problem of dynamic systems.
Decision: Accept (poster)