Learning vector fields of differential equations on manifolds with geometrically constrained operator-valued kernels
Authors: Daning Huang
Hanyang He
John Harlim
Yan Li
Venue: ICLR 2025 Spotlight (2025)
Keywords: Dynamics on manifolds
Operator-valued kernel
Geometry-preserving time integration
Ordinary differential equations
Primary area: learning on time series and dynamical systems
TL;DR: We present a geometrically constrained kernel and a geometry-preserving ODE solver for learning ODEs on manifolds, improving computational efficiency and long-term prediction accuracy in high-dimensional systems.

Abstract: We address the problem of learning ordinary differential equations (ODEs) on manifolds. Existing machine learning methods, particularly those using neural networks, often struggle with high computational demands. To overcome this issue, we introduce a geometrically constrained operator-valued kernel that allows us to represent vector fields on tangent bundles of smooth manifolds. The construction of the kernel imposes the geometric constraints that are estimated from the data and ensures the computational feasibility for learning high dimensional systems of ODEs. Once the vector fields are estimated, e.g., by the kernel ridge regression, we need an ODE solver that guarantees the solution to stay on (or close to) the manifold. To overcome this issue, we propose a geometry-preserving ODE solver that approximates the exponential maps corresponding to the ODE solutions.  We deduce a theoretical error bound for the proposed solver that guarantees the approximate solutions to lie on the manifold in the limit of large data. We verify the effectiveness of the proposed approach on high-dimensional dynamical systems, including the cavity flow problem, the beating and travelling waves in Kuramoto-Sivashinsky equations, and the reaction-diffusion dynamics.

Published as a conference paper at ICLR 2025
LEARNING VECTOR FIELDS OF DIFFERENTIAL EQUATIONS ON MANIFOLDS WITH GEOMETRICALLY CONSTRAINED
OPERATOR- SERRATOR-VALUED KERNELS
Daning Huang
Hanyang He
Department of Aerospace Engineering The Pennsylvania State University University Park, PA 16802, USA dan i innepsg.uuee .ddu
Department of Electrical Engineering
The Pennsylvania State University University Park, PA 16802, USA
hfh531 www.ppppee...e
John Harlim
Yan Li
Department of Mathematics,
Department of Electrical Engineering The Pennsylvania State University University Park, PA 16802, USA
Department of Meteorology and Atmospheric Science Institute for Computational and Data Sciences The Pennsylvania State University
eg115225psu.ee
University Park, PA 16802, USA jharl marlimmpsuuueu
ABSTRACT
We address the problem of learning ordinary differential equations (ODEs) manifolds. Existing machine learning methods, particularly those using neural networks, often struggle with high computational demands. To overcome this issue, we introduce geometrically constrained operator-valued kernel that allows us to represent vector fields on tangent bundles of smooth manifolds. The construction of the kernel imposes the geometric constraints that are es- timated from the data and ensures the computational feasibility for learning high dimensional systems of ODEs. Once the vector fields are estimated, e.g., by the kernel ridge regression, we need an ODE solver that guarantees the solu- tion to stay (or close to) the manifold. To overcome this issue, we propose a geometry- seometry-preserring ODE solver that approximates the exponential maps corresponding to the ODE solutions. We deduce theoretical error bound for the proposed solver that guarantees the approximate solutions to lie on the manifold in the limit of large data. We verify the effectiveness of the proposed approach on nigh--immeniiona dynamical systems, including the cavity flow problem, the beating and travelling waves in Kuramoto-Sivashiee equations, and the reaction- diffusionn dynamics.
INTRODUCTION
In this paper, consider the problem of learning ODE whose solutions lie on manifold, which arises from wide range of applications, from mechanical multibody systems to electrical circuit simulation and power systems (see the references in [Ascher & Petzold (1998) Kunkel (2006)), where the system of ODEs manifolds is formulated based on Differential -aanntiil--lgeee Algebraic Equations [Rheinboldi 1 1984). One of the main challenges in this problem that the underlying manifold (geometric) constraints
not explicitly known and need be uncovered from the data.
A popular solution to this problem to employ a nonlinear dimensionality reduction approach, such as autoencoders, to represent the geometrical constraint of the dynamics, i.e., the manifold. A typical strategy is to learn a -ow-dimensional latent space to represent the original data and learn the dynamics in the latent space; the dynamics represented using, e.g., neural-network (NN) based discrete- time mapping (Linot Graham 2020), recurrent neural network (Maulik et al.
2021 Vlachas et al. 2022), pouuence-to-eqqeeee mapping (Wu al. 2024), Latent Dynamics

---

Published as a conference paper at ICLR 2025
Network |Regazzoni al. 2024), and SINDy Fukami al. 2021 Lin et al. 2024 The main pitfall of this class of methods is that the latent space only provides global parametrization of the manifold, whose dimension is typically larger than the intrinsic dimension of the manifold, and the dynamics are not guaranteed to learn the exponential maps (for discrete-time model) or vector fields (for continuous- titti model) of the manifold. As result, the predicted trajectories may deviate from the manifold, limiting the long-term prediction accuracy, as we shall see in several numerical examples with NN baselines. Furthermore, the computational cost train such NN-based nonlinear estimators is known to be expensive. In our numerical comparison with the proposed linear method, we found that while their testing times are comparable, the training time for NN models is over 800 times slower a simple example with two ambient dimensions.
Contribution. This motivates us to develop a linear estimator of the vector fields on the tan- gent bundles of smooth manifolds. Motivated by SINDy algorithm (Brunton et al. 2016 we construct a geometrically constrained dictionary to represent the unknown vector fields where these constraints will be approximated from the available point cloud data induced by the ob- served time series of the dynamical systems. Since such a dictionary is subjected to the curse of dimensionality, we employ the standard "kernel trick" mitigate this issue. Unlike previous works in non-manifold setting with scalar-valued kernels [Baddoo et al. 2022 Yang et al. 2024) and krig- ing/ nug/Gaussian process (Glaz et al. 2010), the geometrical constraints give rise an operator-values kernel that leverages the intrinsic dimension of the manifold to enable practical implementation. To numerically integrate the system ODEs with the estimated vector fields, we devise an ODE solver that guarantees the solutions to be on the manifolds in the limit of large data. We will demonstrate the effectiveness this approach numerically on several high-dimensional test prob- lems, including the cavity flow problem, the beating and travelling waves in Seramoto--Svashing equations, and reaction -diftiondiiffssion dynamics.
Paper organization. The remainder of this paper is organized as follows: In Section 2 discuss the geometrically constrained dictionary, extending the SINDy approach and generalizing it to an operator-value kernel to mitigate high dimensional problems. In Section 3 we introduce shometry-preseriing time integration scheme, provide an illuminating example that motivate this integrator, and discuss the convergence properties. In Section 4 we discuss closely related approaches that will be used as a baseline to quantify the performance of the proposed approaches documented in Section 5 In Section we give a brief summary. We supplement the paper with five appendices that report the detailed technical numerical tools needed in the algorithm, computational complexity, the theoretical proofs, and additional numerical results.
GEOMETRICALLY CONSTRAINED DICTIONARY
Consider dynamical systems governed by a system of ODEs,
x = f(x),
RER"
(1)
where the vector field f: R 1 - R" The Sparse Identification of Nonlinear Dynamics (SINDy)
approach (Wang et al. 2011 Brunton et al. 2016 to approximate components of (f',.... f")
by sparse regression on set of appropriate basis functions. Typical choices of basis functions can be polynomials and/or trigonometric functions. An example of polynomial dictionary is,
( (1
x
(x2)T ... )ERR
where x := x : j
jnj. Given a set of labeled training data K//E,,,E,....N1 where the subscript index denotes temporal information, x(t) and := x(t;), the SINDy approach is to approximate f*(x)= = fc (x:(e)) == 0(x),, with coefficients & e R ob- tained from solving the following sparse regression problem,
N
E = arr argmin Σ
|| - (xxx A||
(2)
where A > is a sparsity parameter, f, ...,, f(", and == ,..,, e eR"m by default || means l2 norm. There two key issues with this approach as it stands. First, the method is
2

---

Published as a conference paper at ICLR 2025
sensitive to the choice of dictionary. If the space spanned by the dictionary does not encompass the underlying function, the estimated vector field will not be accurate when evaluated on new sample data from the same distribution. Second, this method is computationally impractical as increases since the size of the dictionary increases (exponentially) as function of n. Particularly, if the dictionary consists of monomials of degree up p, then the size of the dictionary is p"
Let us now focus on the first issue for a class of dynamics where the solutions lie on a d- dimensional Riemannian sub-manifold M RR", In this context, the vector field f E (M) is a map f: M - TM that identifies the state x M to a vector in the tangent space f(x) E T.M. Denote the bases of the tangent space Xx, M = æd and the normal space as columns of E gn and N; E (x=x(=-d), respectively. We note that these basis vectors can be identified from point cloud data using the local SVD technique (Donoho & Grimes 2003 Zhang & Zha 2004 higher-order methods (Jiang et al. 2024). See Appendix(Ajov A for the details. For the remainder of this paper, we denote 1 and N, to be the point cloud approximation T, and N,, respectively.
Let the matrix P(x) : R - T&M = R" be an orthogonal projection to the local tangent space at x. One can show that P(x;) T where columns Tf E RRadd forms set of orthonormal vectors that span M. With this background, since f(x) T.M, it is clear that P(x)(x)= f(x) under the π n--iiensional Euclidean inner product. Practically, when the manifold is unknown, one can approximate P(x,) by p(x,) T, Based on this information, we propose the following modification on the SINDy dictionary for modeling the vector field,
(3) where 0(x) e nx is a block diagonal matrix with 0(x) as the diagonal block component, and the coefficients È e Ru are obtained by fitting the model to the observed vector field, x,. In
f(x) = (x; P(x)O(x)
practice, when the available data are only the time series X (x one needs to approximate the derivatives. We consider approximating &; with y; r. [x:-1 Xt)
As we mentioned before, another challenge with this dictionary is that it is subjected to the curse of dimensionality, that the number of candidate functions grows exponentially and eventually becomes computationally intractable as the dimension of states increases. To mitigate this problem, we propose an operator-valued kernel deduced from the dictionary in (3] that allows the vector field to lie on the tangent bundle in the limit of large data with a compact model having rank equal to the intrinsic dimension of the manifold.
In the remainder of this section, we first use the kernel trick to motivate a Geometrically con- strained Multivariate Kernel Ridge Regression (GMKRR) model in the ambient space. Then we formalize the GMKRR model rigorously using the Reproducing Kernel Hilbert Space (RKHS) with a family of operator-valued kernels. Lastly, we manipulate the GMKRR model into the intrinsic space to enable practical computational implementation.
2.1
KENNELIZATION OF THE GEOMETRICALLY CONSTRAINED DICTIONARY
While kernel regression with lg regularization (0 < q $ 1) has been studied extensively (see Shi et al. 2019 and the references therein), it is computationally much simpler to employ & regularization, which will be the focus in this paper. Specifically, we focus on the following as the primal form,
problem modifying [
argmin ||y VE| + +
(4)
=
where y = [y .y. ....y y E n with y; = TT (x;+1 x,)AAt, and V [w]. v...'.
WWT.V..... , TTT
m with vi (x;) 1(x;) R
Next, introduce the dual variable a RNN so that E = T a, and the dual form of 4 is
à = argmin || - a 1|
2 VT al) = argmin ||y Ka || + d Ka
(5)
where the gram matrix K T N and its jith block is K, (xx)xx( = xxx,,x)
Ran The solution the dual form a (AI+ K) - and the predictive model for new input x is given by,
fx (x)= = w k(x) K) K)
(6)
3

---

Published as a conference paper at ICLR 2025
where k(x) [k(x,)()),)))))p)))))))))))))))) k(x, k(x,x)))lll...........l..k.)l))) Since wxx) ((x) R
we have
x(x,x') = 1 (x) (x)) = PXX) P(x') p(x))
(7)
where p(x,x)) 0x)0((x) x(X, ERIX..
and the last equality is because 0xx)(())
Up this point the regression problem with is converted GMKRR problem. The aainetrrcally-costrainee function & in 7 is used as a matrix-valued kernel and constructed from finite of candidate functions defined in the ambient space.
2.2 THE RKHS OF THE INTRINSIC GMKRR MODEL.
In the following, we generalize the GMKRR model a family of matrix kernel functions that may include infinitely many candidate functions via the construction of a X(M)-valued RKHS H.
Definition 2.1. Let X be a non-empty set, W separable Hilbert space with inner product &,, (,)), and LW) Banach space of bounded linear operators on A function k: x X (WW SPD if ffr any pair (x,x)) X X, x(x,x)" = k(x', x)) and (2) for any finite set of points (x: in X and if I Wn f -i,==1 (())) 0. The function k operator-valuee kernel on X and W.
Definition 2.2. Following the notation from the previous definition, for each x E X and f,g W, define fx'')) k(x,x)) for allx' X. For Σ fx, and g' -i RxrBii, define the inner
product, (f, (f.g), v = -ijj (xx,,xx)))) Then H = span fix X, W) forms an RKHS with reproducing kernel k. The RKHS has the reproducing property that ( ((xx).); (f(.),(,, Je
Definition 2.3. The inner product ):,-),e also induces the RKHS norm, llf = / (F for all
IN Kx, G. When W is a n-dimensional Euclidean space, [f vf Kf, where f= [8] .2 ....))
and E RnNanN with the ji block Kij k(x(,x)).
Lemma 2.1. Consider a function k: M - where R x R - R is scalar-valued kernel. Then k is an operator-valued kernel.
- I (MMM), defined as k(x, x') = p(x, (xx) P(x'),
See Appendix D for the proof of Lemma 2.1
The operator-valuee kernel k forms the desired v(Muev RKHS, denoted as H... practice, we can use any SPD kernels such as the squared exponential (SE) kernel x(x,x') exp
[ || x|| /r) or the Matérn kernels (see Appendix E.1).
The function in is special case of the operator-valued kernel on M and X(M). Subsequently, the GMKRR model is reformulated via HM. Given a dataset y(); i= the unknown vector field fE Mys is parametrized f(x) Σ x(x; k(x) where [ax] a ........a.. TJ" is determined by minimizing the following objective function,
j(f- Σ
- f(x)) + =ly-Kall Axx Kax Ta,,
which is the same optimization problem solved in the previous section, and the solution is given by G. However, now the new formulation admits a family of operator-valuee kernels, that may involve an infinite set of candidate functions.
2.3 CONVERSION TO INTRINSIC SPACE
Subsequently, formulate the GMKRR model that the predictive model is effectively defined the intrinsic space and becomes computationally tractable to train and evaluate. Note that P(x)= = f the operator-valued kernel & in ambient space can be rewritten
k(x,x') = p(x
T fxx(xx))
where Ox = = Tx E dddd and
xxxxx)=((xxx))ff) = p(xx'') ( Emaxd
(8)
4

---

Published as a conference paper at ICLR 2025
Using (8) the gram matrix K in ambient space is decomposed K RT where I E RRnNsdN is a block diagonal matrix with diagonal block entries Tú, ,2....,TN, and R E RdNxdN is NxN block matrix with the (i. tt block Rij- = r(X;,))) p(x; xx) Similarly,
xx) = Ty (x, (xx,x2), ".. ... (X,xxx) Txr(x)
(9)
Using the above decompositions, the GMKRR formulation
6 in ambient space is converted to,
(x) (AI K)- (=x(())5'''
(10)
where the Woodbury identity is used in the last equality.
In the intrinsic GMKRR formulation (TP. the matrix (AI+ R) is of dimension dN x dN and com- putationally tractable invert, especially d is small; this regardless of the ambient dimension n. Furthermore, the term xx guarantees that the vector fields lie on the local tangent space of the underlying manifold at the limit of large data.
Lastly, we briefly discuss the intrinsic GMKRR model from an RKHS point of view. First, it can be proved that the function r in 8 is an operator- kernel on M and L(()) (the proof is the intrinsic operator-valued kernel. Then, induces an RKHS and the corresponding GMKRR model in the intrinsic space. The GMKRR is effectively applied
similar to that of Lemma 2.1 r is referred
modified dataset ((x, y; T / where the modified label ý; is the vector field expressed in the local tangent space xj. the kernel r(x,x', (1) p is chosen to be the Diffusion Map kernel and (2) the pairs of data points (x.x)) are sufficiently close so that Oxx = fx is always orthogonal, then the RKHS induced by r subset of 1.2 (M()) that spanned by smooth figenvector-fields of the Connection Laplacians (Singer & Wu 2012).
3
TIME INTEGRATOR
D-0.1
D=0.3
While standard ODE solvers such as Runge- Kutta methods often empirically produce solutions that close enough to the manifold (i.e., a manifold invariant scheme) for suf- ficiently small time step, it well known that the invariant manifold property is only valid when the solvers are em- ployed on special class manifolds (Calvo al. 1996).
D-0.5
0.5
0.0
-0.5
Various ODE solvers on manifolds have been proposed in
-1.0
literature, see Hairer 2011]; Crouch Grossman 1993
for general vector fields and Leimkuhler & Patrick 1996
1.0
1.5
for Hamiltonian systems. In this section, we will illustrate this issue simple example, propose normal correc- tion (NC) to the classical explicit Euler scheme which we call Euler+NC in the remainder of this paper, and provide a convergence study. This approach can be viewed as a realization of the local coordinate approach (see Section I11.2 in Hairer 2011) with local parameterization being estimated by GMLS. The proposed normal correction ap- proximates all of the higher-order terms in an exponential map, expx, (f(x;)A)) including the second fundamental form, and is computationally more attractive than the clas- sical Taylor's method that requires the derivatives of the estimated vector fields.
(a) Parameterized manifolds.
D-0.1
D-0.5
NODE
UpNet
fruth
110
(b) Predictions at a long time period.
Figure 1: Dynamics on a series of 1D manifolds.
3.1 A MOTIVATING EXAMPLE
Consider a scalar ODE: Ò = (0), 0(0) = 0, whose solution has period 2f, and embed the solution in 2D ambient space by (x), X2) (0) cos(0), r r00) sin(O)), where r(0) Coos(OO), and series of values, where the neighboring
The 2D embedding is illustrated in Fig. Ia for points separate apart by one step size Af 0.04. When D 0, the 1D manifold is a "simple" unit circle. As increases, the manifold becomes more distorted, which may pose a challenge in solving the dynamics.
5

---

Published as a conference paper at ICLR 2025
For each value of D, GMKRR model is trained using the trajectory data obtained in the ambient space (details in Appendix(E.1] Subsequently, the GMKRR model is solved using standard RK2 and the proposed Euler+NC method. Meanwhile, neural ODE (NODE) model (Chen et al. 2018 and a Latent Dynamics Networks (L.DNet) model (Regazzoni et al. 2024 are trained as baselines. we show the prediction quality for x1 time interval [96, 128] for D = 0.1 and 0.5. From
In Fig. 1b
Fig. lb we note that the proposed Euler+NC consistently preserves the geometrical constraint; yet, RK2 solutions deviate from the truth. The trajectories of both NODE and LDNet are off the
manifold for D = 0.1. For D 0.5, the NODE achieves a similar accuracy as GMKRR with Buler while LDNet is very accurate (see Table 3 for the RMSE comparison). Both NODE and Nuler+NC develop phase shifts where the former is less accurate relative to the latter. This result suggests the robustness of the proposed solver over the nonlinear methods, that are much more expensive to train. Specifically, for I.Netts, the results above are obtained using twice the size of training data compared to others (we fail to train competitive LDNet using the same data size) and we also found that LDNet with ID latent variable (which is the intrinsic dimension for this problem) does not work. See the complete results in Appendix(E.] where we also reported simulations for comparisons of training and prediction times,
other values of D, their noot-mean-soe ceaenn-suureeeeeeeee comparisons of results with Gaussian and Matérn kernels for p. and the numerical verification of the theoretical convergence rate reported in Theorem 3.1
3.2 EULER METHOD WITH NORMAL CORRECTION
N
T' E M
Tangent increment f(x) )åt Tr Jadt
To solve ODEs on manifold, we consider an integra- tor based the following identity:
X/+1
expx, (f(
Exact exponential map X1+1 exp, (x)))))
Normal Increment (T ((a) (r )a) (At)
=X, + T.T. T (X111
(xi+1 N
X;), (1I)
X1+1
fxyl&r
O((A)))
where the exponential map expx, Txx (M) M is
defined as a map of the vector field f alongg the lo- cal tangent space Tx, .M plus a second-ored -order correc- tion in the normal direction, which can be approxi- mated as a function of f (see Fi
Figure 2: Dynamics decomposed on the manifold.
Based on the identity I once f, is trained (e-g, by the GMKRR model in 100), the solution xe(()1)) 1) to the following approximate differential system,
X (1) = (x)) 1> fn. (x() R
is approximated by &c,++ through the following steps:
1. Use the algorithm in Appendix Aoto approximate bases of the local tangent and normal spaces of & We denote the estimates as Te./ and Ñe.
2. Let (x,,I...., xj,x] X be the kth nearest neighbors of &f. We model the normal compo- nent of by solving the following local regression problem,
8: mig
Σ g J=1
IN (x.) - (A0) geggg g (i. -/,1/11
(12)
Notice that &: Rdd n- takes input from the local coordinates basis Te.i. To solve this re- gression problem, we will employ the Generalized Moving Least Squares (GMLS) method that is documented in Appendix(A] In Appendix D.2 can see that g approximates the second fundamental form of M up to order-At.
3. We predict the state at time f 1; At with the following model that approximates
& = & A1f, / 2 & (T (&
(13)
Here, we simply evaluate g as illustrated Figure 2
We should point out that the proposed Euler+NC method has approximately the same complexity as RK2, as RK2 involves two function evaluations and each requires one local to approximate 1. In Euler+NC, we also employ two peast-squares fittings, one in Step 1 and the other in Step 2. See Appendix for the complexity of the three steps above.
fitting
6

---

Published as a conference paper at ICLR 2025
3.3 CONVERGENCE THEORY
In the following result, we state the error bounds predicting the state Xn, =xxx(())=x(])) for some fixed T> 0 with a general approximate vector field f under the numerical integrator (13). this derivation, assume that the local tangent and normal spaces are approximated with the GMLS method discussed in Appendix(A with a polynomial of degree- e, provided that the manifold M sufficiently smooth, c with l 2 Simultaneously, g, attained by solving [2] is also a polynomial of order-e.
Theorem 3.1. Let the solution of [ lie on a smooth (or at least C manifold M, with 2. Let be Lipschitz function that approximatest. such that If- 1. O(e) ase O. Here, μ denotes an invariant measure of (22). Let .x.] be independent and identically distributed (Li.d.) samples of μ that is absolutely continuous with respect to the volume measure, du = qol with
O < Qmin q S qmax. Then with probability higher than 1- N
N) 4
v Eu [Ixxn - 11 1 C E Ar
N
where T nåt. Here the constant C > 0 is independent of and At. See Proposition D.2 for the detailed expression of this constant.
See Appendix D.2 for the proof. We should point out that E is the bias, which can be specified further in terms of N, l, d, depending on the class of models being used to approximate f. For example, error bounds for KRR with i.i.d. data are reported in (Caponnetto & De Vito 2007] Steinwart al. 2009 The last term in the error bound highlights that the error induced by the proposed solver in (33 is an order-one scheme with an additional factor contributed by errors in the estimation of the tangent space. More work needs be done for non-i.i.d. data, especially when the training data consists of time series of solutions.
RELATED WORK
Learning dynamics with geometrical constraints has been an active research field for decades. Early studies considered constraining the high-dimensional dynamics in linear subspace to obtain oww-dimessional model [Benner et al. 2015). However, such a linear subspace approach is not adequate to represent even simple manifolds such as the sphere and torus. Recent approaches in approximating vector fields on the tangent bundles can be categorized into two classes.
The first class either assumes the knowledge of the manifolds (Elamvazhuthi et al. 2023 Jeong et al. 2023 or employs an empirical procedure that relies on the domain knowledge to globally parameterize the manifolds [Loiseau et al. 2021 Callaham et al. 2022). In the latter approach, the dynamics are identified using SINDy with a special set of polynomials that are only appropriate for describing simple manifolds, such as the torus. Since our method shares a similar vein to this approach, which they called Manifold-SIN SINNYy we provide numerical comparison of their example, the cavity flow problem, in Section 5 We shall see that our method yields significantly more accurate short and long-teriy perm- predictions of the dynamics.
The second class of methods addresses arbitrary manifold structure. In Crosskey & Maggioni ( & Graham 2002 Zeng et al. (2024, atlas of charts of manifold is learned from the data, and the dynamics are learned within each chart so that the dynamics are guaranteed to evolve on the manifold. Specifically, CANDyMan uses k-menss algorithm cluster the data, and employs an autoencoder to represent data clusters as an atlas of the charts and feed-forward neural networks model the dynamics the latent variable induced by the local chart. [oou al. (2020) employed a similar strategy but learned the atlas and dynamics simultaneously through neural ODE framework. The main drawback of this class of methods is the complicated model form, which involves not only the charts and local dynamics but also the transition maps among the charts. In Section we will compare GMKRR with one of the methods in this class, the CANDyMan proposed in Floryan Graham 20222), , several examples in their paper, including the beating and travelling waves in Kuramoto Sivashinsky equations, and reaction reaction-diffusiin diffusion dynamics. We shall see how GMKRR outperforms CANDyMan in training efficiency and prediction accuracy.
7

---

Published as a conference paper at ICLR 2025
5 NUMERICAL EXPERIMENTS
In this section, we numerically test the proposed approach several high-dimensional problems.
5.1 2D CAVITY FLOW PROBLEM
The GMKRR algorithm is demonstrated on cavity flow problem (Callaham al. 2022 that have high ambient dimensions and admit quasi periooii dynamics on a 2D torus. Through this example, the predictive accuracy and long-tern numerical stability of GMKRR are benchmarked against conventional SINDy methods. Two SINDy-type methods are considered for comparison: (1) Vanilla SINDy, denoted simply "SINDy" below, and (2) "Manifold SINDy" (M-SINDy) that is proposed in Callaham et al. (2022 to account for the manifold constraint. trajectory of 30000 steps is available, and all models are trained using the first 4000 steps. The details of the dataset and hyperparameters of the three models are provided in Appendix
10
0.05
SINDY
N SINDy
x 10
GHKAR
10
M-SINDy
SINDY
GMKRR (Ours)
1012
150
154220
224296 300
(a) Comparison of error growth.
(b) Comparison of kinetic energy in long- term prediction.
Figure 3: Results 2D cavity flow.
The prediction errors of the three models are compared in Fig. 3a where the solid lines and the shades show the evolution of the mean and range of RMSE over the test trajectories, respectively. The error in SINDy is large, and maintain at similar level of error through the examined time horizon. The error in M-SINDy grows rapidly; this is attributed the phase shift in the predicted dynamics. In addition, the M- M-SINDy predictions start with a non-negligible error at the beginning, this is due to the error of the manifold representation learned from data. As for GMKRR, the error maintains as low as 1/6 of the M- SINDD errors, with also a smaller standard deviation; such small error is granted from the enforcement of manifold constraint every time step.
Next, the long-ten term numerical stability of the three methods is compared. In the cavity flow problem, the kinetic energy, computed as E |xl fluctuates and maintains constant amplitude and range, which indicates that the flow states stay a now-dimensional manifold. Here, all the models are solved from the first data point for the entire 30000 time steps. The time evolution is shown in Fig. 3b The SINDy prediction gives completely wrong trend, which also indicates that does not satisly the manifold constraints. Both M-SINDy and GMKRR capture the fluctuation amplitudes well; particularly GMKRR follows the truth accurately up to f 150 when M- MINDDy already shows significant deviations. Furthermore, M- M- SINDY shows spurious sigh-frequency fluctuations throughout the simulation, that appears be non- non-physical. Admittedly, GMKRR ends up with phase shift over long-term prediction, and this is arguably numerically unavoidable. However, the comparison above highlights the capability of GMKRR to produce accurate long-term predictions using a limited amount of data, outperforming SINDy and M-SINDy.
5.2
MORE HIGH- BICH-DIMENSIONAL PROBLEMS
Next, we benchmark GMKRR against the digh-dimmnsional examples developed in CANDyMan (Floryan Graham 2022), that have well-defined manifold structures. These examples include the beating and travelling dynamics of the ID Kuramoto-Sivan -Sivashinkky (KS) equation and the 2D reaction-diffuue diffusion (RD) dynamics (Fig. aa)-))). For the KS travelling case, we also consider two additional baselines, NODE and LDNet, that admit oowddimessional (latent) structures but do not exactly account for intrinsic dimensions. The Isomap algorithm (Balasubramanian & Schwartz 2002) employed visualize the intrinsic coordinates of these dynamics (Fig. 4(e-h)); note that
8

---

Published as a conference paper at ICLR 2025
Isomap is only for visualization purposes and its results are not used in the subsequent model training and prediction. Furthermore, the comparison of prediction errors over time between GMKRR and the baselines is shown in Fig. Ii-)) Finally, the computational costs and the error metrics of the models are provided in Table where "RMSE" and "Max Error" refer to the average and maximum of RMSF along the trajectory. Throughout this section, the GMKRR model uses simple local SVD to estimate the tangent space, unless otherwise specified. More details on model hyperparameters, training, and result comparison are provided in Appendices E.3 and E.4
KS beating dynamics. The dynamic oscillates single frequency and lie 1D submanifold in a 64-dimensional ambient space Ke)). This is a relatively easy case for both GMKRR and CANDyMan, but the error of GMKRR one order of magnitude smaller than that of CAn throughout the prediction horizon (Fig. 4i)). Furthermore, the cost for training GMKRR is over 3000 times lower than CANDyMan.
KS travelling dynamics. The dynamic has two timescales that differ by factor of approximately 200. The fast timescale corresponds to the oscillation of the spatial wave shape (Fig. 4b)), while the slow timescale corresponds to the translation of the wave in space (Fig. [e). The combination of the two timescales results in a 3D manifold (Fig. 4(f-g)). Learning the travelling dynamics is more challenging than in the previous case because the complete manifold cannot be covered without including at least full period of the long timescale that spans a large number of data samples (approximately 5000 steps). Here, five models are considered: (1) CANDyMan: The model involves a special data preprocessing step based on prior knowledge (see Appendix E.3 so that the fast and slow timescale dynamics are learned separately, and only the first 100 steps are needed for training. (2) NODE: The modeling strategy and preprocessing procedure are identical to CAnDyMan, except that the fast dynamic is represented by generic NN without assuming any www-dimensional structure. (3) LDNet: The same setup NODE, except that the fast dynamic is learned in a 2D latent space. (4) GMKRR- The model does not consider any prior knowledge as in the CANDyMan, since such knowledge is not always available.
(5) GMKRR-FFT: As an apple-to-apple comparison against CANDyMan, this model is trained using exactly the same data preprocessing strategy and dataset as CANDyMan.
In the initial period, which covers few short timescales, the GMKRR models clearly outperform CANDyMan and L.Net and are comparable with NODE (Fig. 4f).. Notice that LDNet also admits an encoding-decoding error with large RMSE at initial time. In the long term, the baselines NODE and LDNet, that do not account for intrinsic dimensions, perform the worst; this is not surprising as the geometrical constraints are not enforced during prediction. The error of GMKRR- is similar to CANDyMan, although the former does not leverage the prior knowledge like the latter. Lastly, the error of GMKRR-FFF FFTT is consistently 1-2 orders of magnitude smaller than all the other baselines (Fig. (K).. Furthermore, while GMKRR- FFFT and the baselines are trained using the identical dataset, the former is 2-3 orders of magnitude faster than the latter. Clearly, GMKRR beats most baselines even without any prior knowledge, and considering the prior knowledge only expands the advantage of GMKRR over all the baselines.
Reaction diffussonn dynamics. The dynamics oscillates at a single frequency (Fig. 4 d)), and lives on a 1D submanifold (Fig, 4h)). However, the data is 20402- dimeeeiimaiio and hence may pose a computational challenge to both GMKRR and CANDyMan. First we consider CANDyMan and "GMKRR-Ful Full", which are directly trained using the original 00402-dimensional data. When com- pared CANDyMan, GMKRR- Full achieves (Fig. 40). While the training cost of GMKRR-Full is 1/72 of CANDyMan, the former does take longer time prediction, due to the need for GMLS normal correction in all 20402 dimensions. Next, we devise the "GMKRR-PP" PCA" models to reduce the computational cost. The model first reduces the data by PCA into 6 dimensions, and then learn a GMKRR model in the PCA coordinates. In the prediction phase, the GMKRR first predicts the dynamics in the 6 PCA coordinates and then maps back to the original space by the PCA modes. The "GMMRRR-CA-ATT SVD technique to estimate the tangent space, while the "GMKRR-PCA-T4" model uses 4th-order estimation by GMLS. When compared to GMKRR-FFIL Full, the computational costs are significantly reduced in GMKRR-PCA -CA-TO while maintaining nearly identical predictions, hence in Fig. 41), the curves of GMKRR- Full and GMKRR-PCAA overlap with each other. Lastly, GMKRR-PCAL shows improvement over the other GMKRR models, with mild trade-off between cost and error.
lower error in prediction
model uses the local
9

---

Published as a conference paper at ICLR 2025
KS Beating
KS Travelling
Reaction Coactton-Diifusion
Component
Component
(b)
cc)
(d)
space
space
time
time (fast scale)
time (slow scale)
x
x
(e)
(f)
(g)
(h)
Compactly folded 2D
Dominently 1D
(l)
C4 GANAR CAROH PCAJO
20
10-0
Figure 4: Comparison of GMKRR and CANDyMan: (a-d) data; (e-h) visualizations of low- dimensional manifold structure; (i-I) Prediction errors.
Table Model performance for the CANDyMan cases
CASE
MODEL
Training Time (s)
Prediction
RMSE
Max Error
Time (s)
KS
CANDyMan
48.04
0.6703
0.2744
2.130
Beating
GMKRR
0.01367
0.5133
0.02549
0,1972
CANDyMan
70.10
19.44
2.391
15.22
NODE
395.6
4.166
6,690
22.14
KS
LDNet
77.87
13.53
7.600
22.46
Travelling
GMKRR- Full
11.68
20.82
1.604
12.45
GMKRR- FFT
0.04550
24.41
0.3671
3.088
CANDyMan
265.8
14.59
0.1797
0.5517
Reaction-
GMKRR- Full
3.723
175.4
7.394E-4
2.959E-3
Diffusion
MMKRR-PCA-TO
0.01953
4.786
7.399F:-4
2.961
GMKRR-PCA-T4
0.1187
8.5855
1.216E-4
1.339E-3
CONCLUSIONS
We have developed GMKRR, in tandem with a geometry- wwweeery-ppeeevingggg time integrator, for learning dynamics on manifolds; it overcomes the limitations of existing approaches in terms of computa- tional efficiency and prediction accuracy. The GMKRR involves an intrinsic kernel formulation, so that its training and prediction are highly scalable for data high--immnsional ambient spaces. The empirical results show that our models outperform the classical SINDy-based and neural- network- based methods significantly, especially in the long-term prediction accuracy of dynamics on manifolds. We also contribute theoretically by providing the error bound of the prediction by the proposed methods. Future work may focus on improving the scalability to large datasets using, e.g.. pseudo-inpil inntt strategies, extending to dynamics with parameters and inputs, and extending the framework to combat noisy data, leveraging ideas from Candès et al. (001 Messenger &
Bortz 2021 Pasel et al. 2022
10

---

Published as a conference paper at ICLR 2025
REPRODUCIBILITY STATEMENT
We are committed to ensuring the reproducibility of our results. The code used to generate all the figures (except Fig. 2 and the results in T available. At the stage of double-bliig review, the code is provided as a ZIP file in Supplementary Material on OpenReview. Detailed instructions for running the experiments, including the configuration of hyperparameterss and execution of the analysis, are provided in the README. The datasets used in our work are either publicly available or described with clear instructions for downloading and preprocessing the repository. steps required to replicate experiments and analyses are documented to enable seamless reproduction by the reviewers and future researchers.
REFERENCES
Uri M Ascher and Linda R Petzold. Computer methods for ordinary differential equations and differential differential-algerriii algebraic equations. SIAM, 1998.
Peter BBaddoo, Benjamin Herrmann, Beverley McKeon, and Steven I. Brunton. Kernel learning for robust dynamic mode decomposition: linear and nonlinear disambiguation optimization. Proceedings of the Royal Society A, 478(2260):2021883) 2022.
Mukund Balasubramanian and Eric L Schwartz, The isomap algorithm and topological stability. Science, 295(5552):7-7, 2002.
Peter Benner, Serkan Gugercin, and Karen Willcox. A survey of projection-based model reduction methods for parametric dynamical systems. SIAM review, 57(4):483-531, 2015.
Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15):3932-335 3937, 2016.
Jared L Callaham, Steven L Brunton, and Chrrithittep Loiseau. On the role of nonlinear correlations in reduced-ordes modelling. Journal of Fluid Mechanics, 938:A1, 2022.
Maripaz Calvo, Arieh Iserles, and Antonella Zanna. Runge-kutta methods on manifolds. In Numerical Analysis: AR Mitchell 75th Birthday Volume, pp. 57-70. World Scientific, 1996.
Emmanuel] Candès, Xiaodong Li, Ma, and John Wright. Robust principal component analysis.
Journal of the ACM (JACM), 58(3):1-37, 2011.
Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. Foundations of Computational Mathematics, 7(3):331-36) 2007.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Ronald Coifman, Yoel Shkolnisky, Fred J Sigworth, and Amit Singer. Graph laplacian tomography from unknown random projections. IEEE Transactions on Image Processing, 18110)1891-1199, 2008.
Christopher B Croke. Some isoperimetric inequalities and eigenvalue estimates. In Annales scientifiques de I'École normale supérieure, volume 13, pp. 419-435, 1980.
Miles Crosskey and Mauro Maggioni. Atlas: geometric approach to learning high-dimensional stochastic systems near manifolds. Multiscale Modeling & Simulation, 15(1):110-156, 2017.
Peter E Crouch and R Grossman. Numerical integration of ordinary differential equations manifolds. Journal of Nonlinear Science, 3:1-33, 1993.
David L Donoho and Carrie Grimes. Hessian eigenmaps: Locally linear embedding techniques for oigh-dimensional data. Proceedings of the National Academy of Sciences, 100(10):5591-550) 2003.
11

---

Published as a conference paper at ICLR 2025
Karthik Elamvazhuthi, Xuechen Zhang, Samet Oymak, and Fabio Pasqualetti. Learning mani- folds: Universal approximations properties using geometric controllability conditions for neural odes. In Learning for Dynamics and Control Conference, pp. 1-11. PMLR, 2023.
Urban Fasel, Nathan Kutz, Bingni W Brunton, and Steven Brunton. Ensemble-sindy- Robust sparse model discovery in the low-data, high-noise limit, with active learning and control. Proceedings of the Royal Society A, 478(2260) 2022.
Daniel Floryan and Michael D Graham. Data-driven discovery of intrinsic dynamics. Nature Machine Intelligence, ((1))11133112) 2022.
Kai Fukami, Takaaki Murata, Kai Zhang, and Koji Fukagata. Sparse identification of nonlinear dynamics with www-dimmensionalzeddd flow representations. Journal of Fluid Mechanics, 926:A10, 2021.
Bryan Glaz, Li Liu, and Peretz P Friedmann. Reduced-order nonlinear unsteady aerodynamic modeling using surrogate-base based recurrence framework. AIAA journal, 48(10):2418-2429, 2010.
Ernst Hairer. Solving differential equations manifolds. Lecture notes, 2011.
John Harlim, Shixiao Willing Jiang, and John Wilson Peoples. Radial basis approximation of tensor fields on manifolds: from operator estimation to manifold learning. Journal of Machine Learning Research, 24(445):1-85, 2023.
Seungwoo Jeong, Wonjun Ko, Ahmad Wisnu Mulyadi, and Heung-ll Suk. Deep efficient continuous manifold learning for time series modeling. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023.
Shixiao Willing Jiang, Rongji Li, Qile Yan, and John Harlim. Generalized finite difference method on unknown manifolds. Journal of Computational Physics, 502:1 12812, 2024.
Peter Kunkel. cefferenttallalgehne equations: analysis and numerical solution, volume 2. Euro- pean Mathematical Society, 2006.
Ben Leimkuhler and George W Patrick. A symplectic integrator for riemannian manifolds. Journal of Nonlinear Science, 6:367-384, 1996.
X Lin, D Xiao, and F Fang. Data discovery of low dimensional fluid dynamics of turbulent flows. arXiv preprint arXiv:2401.055555 2024.
Alec J Linot and Michael D Graham. Deep learning to discover and predict dynamics on an inertial manifold. Physical Review 101(6)0062299, 2020.
Jean-Chie Christophe Loiseau, Steven Brunton, and Bernd Noack. From the pod- galerkkln method to sparse manifold models. Handhook of Model Order Reduction, Volume 3, 2021.
Ser Aaron Lou, Derek Lim, Isay Katsman, Leo Huang, Qingxuan Jiang, Ser Nam Lim, and Christopher M
De Sa. Neural manifold ordinary differential equations. Advances in Neural Information Processing Systems, 017758--15555 2020.
Romit Maulik, Bethany Lusch, and Prasanna Balaprakash. Reduced-order modeling of advection- dominated systems with recurrent neural networks and convolutional autoencoders. Physics of Fluids, 33(3), 2021.
Daniel A Messenger and David Bortz. Weak sindy for partial differential equations. Journal of Computational Physics, 443:1 10525, 2021.
Davoud Mirzaei, Robert Schaback, and Mehdi Dehghan. On generalized moving least squares and diffuse derivatives. IMA Journal of Numerical Analysis, 3210 2012.
PySINDy. Documentation. Sitmppess/I. Btep:/4___viii____eeeeee http://.awoooo.acaaaaaaaaaaaaaaa accessed 22-Aug
latest/
ht Maw.:66-aainoooddddd 2024. (Online;
12

---

Published as a conference paper at ICLR 2025
Francesco Regazzoni, Stefano Pagani, Matteo Salvador, Luca Dede, and Alfio Quarteroni. Learning the intrinsic dynamics of spatio- temporal processes through latent dynamics networks. Nature Communications, 15(1):1834, 2024.
Werner Rheinboldt. Differentiallalgebbric systems differential equations manifolds. Mathematics of computation, 43(168):473-482, 1984.
Clarence W Rowley and Jerrold E Marsden. Reconstruction equations and the parhunen-loève expansion for systems with symmetry. Physica D: Nonlinear Phenomena, 142(1-2) 12)1-19 2000.
Lei Shi, Xiaolin Huang, Yunlong Feng, and Johan AK Suykens. Sparse kernel regression with coefficient-bastt based l q- regularization. Journal of Machine Learning Research, 20(161):1 2019.
Amit Singer and H-T Wu. Vector diffusion maps and the connection laplacian. Communications on pure and applied mathematics, 65(8):1067-1144, 2012.
Ingo Steinwart, Don R Hush, and Clint Scovel. Optimal rates for regularized least squares regres- sion. Proceedings of the 22nd Annual Conference on Learning Theory, pp. 79-93, 2009.
Pantelis R Vlachas, Georgios Arampatzis, Caroline Uhler, and Petros Koumoutsakos. Multiscale simulations of complex systems by learning their effective dynamics. Nature Machine Intelli- gence, (44)3359-33 366, 2022.
Wen-Xu Wang, Rui Yang, Ying-Cheng Lai, Vassilios Kovanis, and Celso Grebogi. Predicting catas- trophes in nonlinear dynamical systems by compressive sensing, Physical review letters, 106 (15):154101, 2011.
Tao Wu, Xiangyun Gao, Feng An, Xiaotian Sun, Haizhong An, Zhen Su, Shraddha Gupta, Jianxi Gao, and Jürgen Kurths. Predicting multiple observations in complex systems through low- dimensional embeddings. Nature Communications, 15(1):2242, 2024.
Qile Yan, Shixiao Willing Jiang, and John Harlim. Spectral methods for solving elliptic pdes unknown manifolds. Journal of Computational Physics, 486: 12132, 2023.
Lu Yang, Xiuwen Sun, Boumediene Hamzi, Houman Owhadi, and Naiming Xie. Learning dy- namical systems from data: A simple coos-validattin perspective, part Sparse kernel flows for 132 chaotic dynamical systems. Physica D: Nonlinear Phenomena, 460:134070, 2024. ISSN 0167-2789. doi: httpp://oo.oorg///oooili/plaaaa
URL https
www.:.scieeete. www.:ssceeeee...c rect. com wtt.:///cc..cccaaa......////////
Kevin Zeng, Carlos E Pérez De Jesús, Andrew J Fox, and Michael D Graham. Autoencoders for discovering manifold dimension and coordinates in data from complex dynamical systems. Machine Learning: Science and Technology, 5(2):025053, 2024.
Zhenyue Zhang and Hongyuan Zha. Principal manifolds and nonlinear dimensionality reduction via tangent space alignment. journal on scientific computing, 26(1):313-33) 2004.
13

---

Published as a conference paper at ICLR 2025
A A LOCAL PARAMETERIZATION OF SMOOTH MANIFOLDS BY GENERALIZED MOVING LEAST- LEAST-SQUARES METHODS
In this appendix, discuss local SVD method approximate local tangent space, review the Generalized Moving Least Squares (GMLS) method, and its application to approximate a local parameterization of smooth manifolds which in turn gives us higher-order approximation of the local tangent space, and discuss an algorithm for the normal correction for nigh-dimensiona. ambient space.
Let X = | = M c R2 be a set of point cloud data that is available for training. For any base point X, M (Xo can be either in the training data set any new points), we denote its K K-nearest neighbors in X by Sx, t) k=1 c X. By definition, we have Xo,1 Xo to be the base point.
A.1
LOCAL SVD APPROXIMATION TO THE LOCAL TANGENT SPACE AND ITS DIMENSION
The classical local SVD method Donoho Grimes 2003 follows:
Zhang & Zha 2004 approximate T, as
Construct the distance matrix [Di...., Dxl RO*K where d and D, == Xo, Xo-
2. Take a singular value decomposition of = UΣVT. Let To be the leading d' colollmms of which approximates orthonormal basis of 7xp M.
We should point out that this approximation has convergence of order N for K d' n l Kmaxl, where Kmax denotes the maximum principal curvature and for uniformly sam- pled data set size denotes the size of training data (See Remark 9 in Harlim et al. (2023)..
A.2 GENERALIZED MOVING LEAST SQUARES
To attain oigher-order approximation, consider the generalized moving least square method. Generally, consider the regression for a smooth function f (M)) Denote f, ( T We consider using the intrinsic polynomials to approximate smooth func- tions over neighborhood of xg. By intrinsic, we refer to the polynomial defines local chart with base point Xg-
First, let ppdd be the space of intrinsic polynomials with degree up to in d-dimensions at the Xg point Xo, i.e., %0 span((Px). ahaiss)) where (a1...,a)) is the multi-index notation with |a| = 01 + + ag and Pxo.a is the basis polynomial functions defined as
p..d
Pxo. a (x) = 11 (x,) (x (x-x))]
xn)
| l.
1=1
Here, ( are the d tangent vectors at the base point Xo that are estimated by the local SVD
scheme, TO (7) (xp)..., (x)))). By definition, the dimension of the space p!dd is m d
d
%)
For K m, we can define an operator Fp : fxn e - Fplx, E plid such that Ipfy, is the optimal solution of the following squares problem:
Σ (f(xa))) - f(x,,))
min Pepld
(14)
The solution to the least-squares problem
14 can be represented as
FpLN
Σ
ba Px
lais!
where the concatenated coefficients b= (ba:))... ba(m))
satisfy the normal equation,
(Ø b)b o Pxn-
(15)
14

---

Published as a conference paper at ICLR 2025
Notice that here is number of coefficients a such that [a] = and thereafter used the notation {a( j= to denote all possible m multi indices. The above normal equation can be uniquely
identified if o RK*m with components
O = ax) (x0,)), 1sks 1 isjsm.
(16)
full rank matrix.
For uniformly sampled data set X size N, it was shown that the convergence rate of the GMLS approximation with polynomial of order-e to the ath derivative o of order-N - (seeJiang al. 2024). For the tangent space approximation, the function will be the local parameterization of the manifold (as discussed in the next section), and the tangent space will be Jacobian such that 1. Therefore, if we consider GMLS approximation with polynomial ial of order > 1, we effectively have higher-order approximation the local SVD method.
A.3 GMLS APPROXIMATION
LOCAL. PARAMETERIZATION LOCAL TANGENT SPACE
In this section, manifold :: Tr, M = nd M <<R" such that (0) = Xo-
describe the GMLS approximation local parameterization of the
Given the global tangent vectors at the base point Xo, (š) (x))li=.....d position to approximate the normal basis vectors (X,( j= j= Then one can define local coordinate chart for the manifold near the base point using the embedding map
the QR decom-
(x, (0,,0......,00 = X + 1 (X0) (272(xx) +... (X))
2j (0,,02..... ................................ 0 (X,), /=1
Σ
(17)
where (0,,02...,,0f)) E Rd with 0, (xx 1x--xxxx denote the local coordinates and 2 n (Xg)
(x xx-xx Xp) is a smooth function over the (e..0...., 0n) p.lpaane. This local parametrization is known as the Monge-Gaue Gauge parametrization of the manifolds.
We will simply approximate zj R - with the GMLS method function that maps the local coordinates (0),0,....,0)) to point on the manifold M by solving the local least squares problem:
min Σ (zj,e 2, 2)) tycPil ky k=
(18)
where the available data are,
2j.k (Xg) (Xo, -X)).
In this notation, we are fitting model of the following form (intrinsic polynomials),
10 ..00
Σ
bj 00 000000 00.0000.000
laist
where 2y (X,,)) is shorthand notation for & ( Ol) evaluated at the local coordinates of Xo,k. i.e., 0f (xn) (Xo,k Xo) for I,...,..
With this GMISS approximation, we also inherit the approximate tangent space from the columns of the Jacobian Zx defined as in [7 with zj replaced by its estimate 2j.
Dis, ( f......00.l (61 (Xo)
(x0))
(19)
where
(Xx) 7,(x0)
Σ |=1 30; 2, (0,.02...,0aa) (n (x()).
(20)
Effectively, the GMLS approximation improves the approximation (as discussed at the end of the previous subsection) by adding additional terms in the normal components in 200
15

---

Published as a conference paper at ICLR 2025
A.4 GMLS FOR NORMAL CORRECTION HIGH AMBIENT DIMENSIONS
For problems of high ambient dimensions, when the above procedure is applied to the normal correction in the Euler+NC method, computational bottleneck will emerge. That is, at every time step, all the (n d) normal vectors Ñ [ R2, ... nn al R need be computed to form the (z - d) local least squares problems [88 The normal vectors may be obtained by, e.g., QR decomposition when is large. The same issue occurs in the regression problem 12 where one needs Ñ once f
T [2, f2....... ... Žal E m which would cost O(n(n - d)²) and is intractable
attained from the GMLS regression in (20).
To circumvent the explicit computation of the Ñ, employ the following projection-based approach. First we define the projection matrix Q= ÑÑ which maps a point x e to the subspace spanned by N. Equivalently Q=1-10 TT so that the formation of Ñ is avoided, and the matrix-vector product xx=x- T(T x) has low cost of O(dn). Next, the following modified least squares problem defined:
min R l.
I x(xx,)))
(21)
where the data
Q(x0,t xn) ER
and the model x= "",
I 04)= = Σ ba0 00.000000.000 ...002
aa/s/
Problem 21 is essentially n local least squares problems; once it solved, the normal components in 17 are computed collectively as
n-d d Σ
22) Qx(0)....,
In the case where the ambient dimension is extremely high, the solution of 21 can be paral- lelized to further accelerate the computation.
The same approach is used to avoid an explicit computation of Ñ by replacing T with in the formulation above.
ESTIMATION OF THE INTRINSIC DIMENSION OF A MANIFOLD
A critical parameter in GMKRR is the intrinsic dimension of the manifold. Here two algorithms for estimating d from data are provided.
B. LOCAL SVD
Following notation in Appendix let the diagonal components of E, the singular values D, be denoted by 01 02
ox. To determine d, normalize the singular values by defining.
We choose dg = dg())) as an approximation to the local tangent space Tx, M that depends on a fixed threshold 8 (0, chosen such that,
do(0) = argmax
d: 2 2
dzl
isd
This threshold must be tuned to the data, noting that higher thresholds will yield higher dimen- sional tangent spaces and lower thresholds will yield lower dimensional tangent spaces with all the trade-offs that this implies.
16

---

Published as a conference paper at ICLR 2025
B.2 KERNEL-BASH METHOD
An alternative method for dimension estimation that avoids tuning is provided in Coifman et al. (2008). The method relies ad hoc choice of scalar-valued kernel p(x,x';n) with bandwidth parameter 77.
Given dataset of point 12)1-12. NN, define the sum
N
S(n) = N2
One can verify that -o S(n) = I/N and lim S(n) = Coifman et al. (2008) shows that
S(n) x s7 s7d/2 in range of 77, and proposed the following estimation of intrinsic dimension. First define
2dlog(S(n))
d(n) dlog(n)
then the estimated intrinsic dimension is d = dd(n"), where
17 = argmaxd(m)
In the numerical implementation of the above intrinsic dimension estimation method, the deriva- tive d(()) can be estimated analytically, since d(()) / dS an. simply by finite difference. Further- more, we suggest using a modified SE kernel p(x,x';n) = exp (-|1xxx((|(1(2m))
where is the
maximum of pair-wise L2 distances of the dataset; this modification normalizes the distances and makes the search for maximum d more robust.
Finally, as side note, while the above bandwidth parameter 7) is developed for intrinsic dimension estimation, we found empirically that a good initial guess for the bandwidth parameter y in the SE kernel of GMKRR can be computed based on 17""
Y 10(
In the KS and reaction reaction-diffusion cases presented this paper, to perform the cross-validation in GMKRR, we empirically set grid search of, e.g., 0.8y", 1.0y", 1.2y", etc. Yet, it was found that choosing y for the SE kernel is sufficient for accurate GMKRR predictions in some cases.
C COMPLEXITY OF EACH PREDICTION STEP
To attain and Ñ in Step 1, we employ K K-nearest neighbors (knn). generic knn algorithm time complexity OLN(d K)), where denotes the number of nearest neighbors. Subsequently, we solve an SVD of matrix of size n K, where denotes the dimension of ambient space, the time complexity of local SVD O(n max(n, K] min(n, K]3) obtain T and Ñ. Since use the higher-order approximation in 20), GMLS regression is employed in parallel for each co-dimensional coordinate, In each coordinate, we solve linear problem of size m, where m < K denotes the number of monomials in GMLS fitting, so the cost in constructing + a in 15 O(Km²) and inverting O(m²), which is cheaper than the SVD. So, the total cost still dominated by local SVD. For problems with n >> K, see Appendix(A.o A.4 implementation to avoid forming the normal vectors Ñ in the SVD step and the improved normal vector estimates Ñ after we obtained î from the GMLS step. The complexity of these operations is far less than the computations above as reported in Appendix(A.4
A.4Ioo for an efficient numerical
In Step of the algorithm, we employ the GMLS again attain each component of g in parallel. Each component takes O(m² + m²) operations.
The cost of evaluating f, (x) fxr(x) includes the total computational costs of the first step above to attain Tx the meatrix-matrin vectoo vector multiplication, O(Nd² + nd), the evaluation of r(x) in g which is Old? nN), and the cost of evaluating exponentiation in the SE and Matérn kernels.
The total cost of evaluating 13 consists of the total cost of evaluating f (x), the time complexity in multiplying Ng, which n(n d), the O(m² K m²) operations in GMLS fitting from Step 2 to
17

---

Published as a conference paper at ICLR 2025
attain g, and the complexity in evaluating d component of g, which is (n d)dm, since we fit the GMIS d-dimensional model with a total m parameters in each coordinate. we use the implicit approach in Appendix then the multiplication Ng is implicitly done, and the total computation cost reduce from n(n d) + (n - d)dm to ndm.
D PROOFS
D.1 PROOF LEMMA 2.1
Proof. First, for any pair (x,x') M x M, k(x.x)" p(xxx') (P())()) pix xX(x)x k x). Second, for every finite set of points (x) in and (f)) in X(M),
N
Σ ((,,(x
Σ f
=
f/=1
Σ
x
u
Σ 1,j=1
p(x;,,y))
Σ Uik
Σ -1i-11
Σ
2
where u; P(x;)f; [Un, "....., Munll" and the last because that p is positive semi-definite so each term in the summation over k is non-negative. Hence & is an operator-valued kernel. O
D.2 PROOF THEOREM
Consider the following initial value problem involving system of ODEs,
x(r) f(x())), To.
x(1n) Xx,
(22)
where 4 II x C (f) and f: M M TMM R" is Lipschitz continuous. For the discussion below, we assume that (22 has an invariant measure
Assumption D.I. Letf, be an approximation off, such that - [ffffell?? O(c) ase 0. We also denotex, that solves the approximate dynamics,
(1) = (t),
t> for Xe(1a) = Xg-
(23)
Here, fe is assumed to be Lipschitz continuous.
To facilitate the analysis below, we review the following basic result:
Lemma (Gronwall), Let p e L. (10, TI), that is, fo p(t) dt <o.. Also let p z and g, e C((O, 7])
with g nondecreasing. If
(())+((++++ ( p(s)p(s)
then
Φ( (())) [( p(s)dss. (1,
Let's first the error induced by the approximate vector field In the following, the expectation defined with respect the distribution of the initial conditions, X ~ μ.. Particularly, the expectation of the random variable x(r) is understood as integration over the distribution of x(r) that solves 22 with initial distribution Xo - μ.
Proposition D.I. Let the Assumption
D.I be valid. For any >>0.
F [I((T) - =
where C(T) 2L2 2 0, for some constant independent of
18

---

Published as a conference paper at ICLR 2025
Proof. For any fixed Then,
[x(() -x()ll s [f(x(s)) (xx
s ! [f(x(s)) - (ss I [ife((()) (x))x ds
s I [f(((s) (x((s)) ds f" [Ix(s) -x- ds.
Jo
where L>> is the Lipschitz constant for fc. By Gronwall, we have,
(x(() (x))l e [ [[f(x(s)) - - [x
By the Jensen's inequality,
/ ( ff (x((s)) s ' -
(x(1) ((10-xx(ee xx1I x000 = ا
Taking expectation, applying Fubini's theorem and the assumption,
Ep (000)) e L E. [If(x((()) (x((s)) ds ce Pe?
for some c 0.
For convenience of the discussion, let us state the GMLS error bound in terms of the fill distance Mirzaei et al. 2012
hx.M sup min dg (x, x)),,,
where dg. M M RR denotes the geodesic distance.
Lemma D.2. Let X = Xxy c M be a cotusi-unifoomlly distributed data with hx.M ht. Define
M UKEM B(x, C2h a union of geodesic ball over the length C2ho, for some constant C2 > 0. Then
for any real-valued function E C++I(M*) and a that satisfies lal with probability higher than 1 1- N
D f(x) -
fx)( (D) ch f+1-|o1|
X.M
(24)
for all E M and some c>> 0. Here the semi-norm
Ifcenam
I fll 1
max
is defined over M In the error bound above, we used the notation f as the GMLS approximation to D" f using local polynomials up to degree where D denotes general multi iinlti-dimessionall derivative with multiindex a.
Assumption D.2. Let the training data set X be sampled ii.d. with sampling distribution μ. which absolutely continuous with respect to the volume measure, du= Vol, with 0 qmin q 9max-
Then we can show that:
Lemma D.3. Let the Assumption D.2 valid, then the filled distance
Px-p(hx.M > 6) s =xx(xcC))=
where C C(d)/Vol(M) is some constant that depends on the intrinsic dimension of the manifold M.
((() Proof. Suppose hx.M > 8, so there is x M such that min x, dg (x,x;) 2 8. So (x) n Ø. In
other words, eachx; E X M MVB, (x). This has measure M(M 1 Cdd, where the inequality based on the lower bound estimate for the Volume of geodesic ball (Proposition 14 in Croke
(1980)). Hence, ming, dg (x,x;) ) 2 occurs with probability less than - C s exp(-CN8") This completes the proof.
19

---

Published as a conference paper at ICLR 2025
With this estimate, we conclude that with probability higher than 1 x the fill distance is bounded from above,
Vol(M) log(N)
(25)
C(d)
The result above requires the sample data to be independent and identically distributed (i.i.d.).
While in general the data may not be i..d. they time series, the training data for the ID example are drawn i.i.d from the sampling distribution μ (see Appendix
Using the same argument and repeating the proof of Lemma A.2 Yan 2023 one can show that with probability higher than 1- N
2. the separation distance,
log 1/d N2
Ci (d) Omin 2-2 s 9x.M = min dg (x, C2(d) 4max
2 1/)
for some ej (d), C2)))),, so the N,
assumption Lemma (D.2) for a fixed X of size
4x,M S fx, S Cg9.,
remains valid although with constant Cq that grow on the order of Ni/d, So, we conclude that under the Assumption D.2 the GMLS error estimate for any fE C++1 (M*) is given by,
log(N) N
1 - D(x(x)
|
Mct-
(26)
for |a| = and constant Č= c
Vol(M) C(d)
1/d
that depends d and Vol(M) with probability higher
than 1 N-
We now solve the ODE 23 with the following approximatinn (which exactly 13
& = x (&,, ((rn) &x.
(27)
where we have defined := ==
(4;).
Let yr: [0, 1] M < Rd denotes the geodesic that satisfies (0) =x, == X(())) and y = (xe,).. In the discussion below, we denote the true tangent and normal vectors x by Te. and Ne., respectively. Since expx., (f, i)tt) (A(t),, we can employ the following Taylor's expansion the exponential map,
Xe,i1l = expx (fe (x, A1) = (i(At)
1
(100) マ (0) +
x (x,+ A) 2 (fe f (Ar) +0(AA)
(28)
where is the second fundamental form of M, which is quadratic function of (X,,)) that is,
a (f, (X,,i), (x..)) N 5, (x (x,,),
(29)
where V,v y denotes the Euclidean R²" directional derivative of in the direction of and (Xe,)) is an extension of (x,,), such that (x(x)))) (xe,))
On the other hand, it is clear that,
X - xx,- =T, (Xe,/+1 x).)) + N c.1 (Xe,/+1 - Xe.i) = (Xe Ne (Xe,111- (30)
Subtracting 30 from 28 and using 29 and taking NT e,f" replacing with At, we have,
(A1)² N
((,,)))) N
(10++1xx,,)) x.,)) (x,,+11
[xx. (x ((Ar))),
(31)
which motivates the GMLS estimation with polynomial of degree 2.
20

---

Published as a conference paper at ICLR 2025
In general, however, do not have the underlying Te, and Ne.i If we numerically employ the GMLS to approximate the local parameterization :: 12 M with polynomial of degree-s as prescribed in Appendix we have the error bound in 26 with aa = 1. That is, with probability
higher than
Te. Te.illo = (og(N) N
It (ac) :=
where we have defined 5 to simplify the notation. Using the Gram-Schmin formula and basic norm triangle inequalities, one can deduce that the normal vector estimation has the same error bound, Nz./ - Ne.1 ll s 8. This means that,
a (f.(x,,)) (x..))))))))))) = N o,l V (x (x(e.e) No
o, V Lix. (x ) 0(6),
(32)
in the uniform sense. Also,
Te, TT o (Xe,1+1 Xe,)) + Ne, of (Ke,/+1 Xe,1)
xe,+11 - Xe.
(xxx Ne, (Ke,111 xe,i) Te. e,i - -Te, Tei (Xe,1t1 - Xe.).
(33)
Subtracting 33 from 28 using 32 and taking N replacing with At, we will numerically
estimate the following quantity.
(Ar)
e MA))) 0x.
xx,i)
N V F
f (x, ((An) +611).
(34)
Based on the GMLS error bound in 26 we immediately have the following result.
Lemma D.4. Let the Assumption D.2 be valid. Suppose that the solution generated by the system of
ODEs in 22 lie on smooth manifold, ofat least such that (88 is valid. Let be at least C*+ (M*), where s 2 2. Furthermore, let g be the GMLS polynomial of order estimate for E C (M²) in
with f 22. Then, with probability higher than -2. K.
(lgg(N) N
1BCT. ((x,,))) - (e
Cg
(g(e
log(N)
+C,Ar
(4c (x
N
where Gg and C, are the two constants of the error bounds of GMLS estimates for, respectively, g andı, as defined in 26
Proof. By GMLS error estimate, is clear that, by probability greater than 1 N'
log(N) 4
f -gt (e co
Cg
18l0
On the other hand, with probability greater then 1
N.
f Ff fx compa-d
(Ne.. - Ne.1) T (xe,1-1 - xx, 1
No.1 - -N.. no X X loo
C. log(N)
(35)
N
where have used the bound in 26 and the first order Taylor expansion in 28 these two bounds, the proof complete.
Combining
Proposition D.2. Letx; be the solution of with Lipschitzfe. We approximatex, numerically through the solver defined in (27 with in 34 to be the GMIS approximant of g as defined in 31 Let the assumption in Lemma D.4 he valid. Then, with probability higher than 3
N
2min f
log(N) N
Ep [|x 11, CAP2
21

---

Published as a conference paper at ICLR 2025
as N and Af 0, where Č depends qmax Vol(M), T,L, cf M4c (xx, [g
Proof. To establish the convergence, let us state the local truncation error by inserting the true solution of 23 to the discrete solver in (27)
Ti = xc.++
%e.
(x(x,)))
A 1 (xx.i)))
Af
At Ne. NN., (e
(36)
where we have used (88) with quadratic and higher order terms expressed in terms of g using (29 and 34. We first note that this scheme is indeed consistent, that is, with probability higher than
s Ar( No.1 fe (x..)) (e
4
s sr(Ce
log(N)
(glcen(r))
N
(log(N)
+C,
[ (L, [1
where we have used the bound in (D.4 and the fact that Ne.: - -N. = with probability 1 N This bound effectively means that [T; loo = O(áz) as Ar - 0, with constant stated as above.
Following the standard convergence analysis by subtracting 27 from 36 and denoting the difference as E; == Xe.i &e,i, we obtain
1 F + At [f. (x) (Kt.))) + (Ar) Ne./ -
+T
(37)
and using the fact that & and & are Lipschitz continuous and denoting L as the largest Lipschitz constant among these two functions, have
E lloo
LAt 1.3 (A.)) (A)) Elylo + +IT; Al,,
where Č = člgl (rr) Solving this discrete equation, we have,
| (1 LA++ (2 + At Σ (1+ (1+ ((0D)++
For time interval [0, 7] where nAr T, we have,
(1+ LAr+ (2 L fui-Att s Lnndr eLT
With this identity, it is clear that, with probability higher than 1
f [+lle)
|| ||Elll+ CTAI log(N)
| LT
Since no error is committed the initial condition, Eg and we conclude that,
(+1s
log(N) N
At
as Ar and N Taking the expectation of square under μ and using the upper bound of the density in the Assumption D.2 the proof complete.
With these results, the proof of Theorem 1 is rather straightforwaar, owing to the fact that,
/EE. -
V Ep [|x,, Xe, 112] EEp ||x & [1]
Ep ||x, Xe, 11 [
Ep [lxx, &e,
(38)
and the upper bounds in Propositions Proposition D.I we note the notations X,, x(t,, x(T) and x == (t,),, = ( where f, nAr=T.
D.I and D.2 with e When using the result from
22

---

Published as a conference paper at ICLR 2025
Table Model parameters of GMKRR for the 1D cases
CASE
SE Kernel
Order of Tangent
Band- width,
GMLS for Space Es- NC,
timation,
==
1.0
3
D 0.1
1.0
5
D = 0.3
1.0
5
D
1.0
4
E ADDITIONAL EXPERIMENTAL DETAILS
All the reported results, including computational costs, are generated on a MacBook Pro Intel Core 17. All the training and prediction
CPU.
E.1 A ID DYNAMICS EXAMPLE
Dataset.
The solution the scalar ODE,
3 2
0 cos(O), 0(0) =0
has closed-form solution
0(t) 2 tan
(tan(V/5/4)//
and using the given embedding
(x),x) (r(0) ((s))))))))
one can evaluate the exact trajectory (Xi, x2) at any / 0, that lies on the manifold.
For given (K, D) combination, the trajectory is sampled over [0, 128] with step size of At = 0.04, which ends up with 3201 steps. The first 200 steps, which covers approximately 1.5 periods, are used for training; the test starts from = 0 and covers all the 3201 steps.
Hyperparameters of the models. In Fig. the GMKRR model uses the SE kernel. The hyperpa- rameters include the kernel bandwidth y, order of GMLS for normal correction É, and order of GMLS for estimation of tangent space s; is only used for Euler+NC. A validation procedure is used to identify the best hyperparameters for the GMKRR model combined with a particular ODE solver as follows. grid search is performed over (y. E 1,.1, 1.0.5.0, 10.0) x (3, 4,,44556,, 7) (3,4,5, 7, and for each combination (y, C, s), the dynamics is solved using the chosen solver for 50 steps starting from 0 and the RMSE is recorded. When initial value of y were selected, more rounds of refined grid search over (y, s) would be performed until no further improvement in prediction is observed. For example, if y = 1.0 is selected from (0.1, 0.5, 1.0, 10.0), then narrower range (y, s) E (0.8, 0.9, 1.0, 1.2) x (3, 4,5, 6, 7) x (3,4,5,6, 6, would be considered next. The hyperparameters associated with the lowest RMSE are used for the GMKRR model. practice, GMKRR model prediction is found to be insensitive to the hyperparameters when Euler+NC is used, Hence, for (K, D) combination, the GMKRR model trained using the optimal y and for RK2 and the optimal for Euler+NC. The detailed choices are provided in Table 2 This choice of hyperparameters gives a slight disadvantage on Euler+NC. However, as one has seen in Fig, [b) and Fig. for more complete prediction horizon, RK2 still fails miserably for longer prediction time all cases, while Nuler+NC consistently predicts the periods and amplitudes even though only first-order time accurate.
Comparison with other baselines. In Fig. 1b (and Fig. 5 for the entire prediction time and other values of D), the NODE method (Chen et al., 2018 and the LDNet method (Regazzoni et al. 2024) are included the baselines. The codes for these methods publicly available
23

---

Published as a conference paper at ICLR 2025
D=0.0
1
0
D=0.1
2.5
0.0
0.3
Euler+
1
W
RK2
NODE
x
0
LDNet
Truth
-1
D-0.5
1
*
0
1
0
20
40
60
80
100
120
Figure 5: Predictions time for various values of D (A more complete depiction of Figure 1b
online (NODE and and LDNe F. For NODE, the network architecture consists of two inputs and two outputs, with two hidden layers of 32 neurons each and PReLU activation function. The NODE is inteeinterrtted by RK4. The model is trained using Adam optimizer with learning rate I x 10 over 1000 iterations, using batch size of 20 trajectory segments; in each batch, the segments are randomly selected from the same training dataset used by GMKRR (one trajectory of 200 steps). For LDNet, the network architecture consists of dynamics network and a reconstruction network. The former has two inputs (number of states, i.e., ambient dimensions) and two outputs (latent variables), with a hidden layer of 20 neurons and tanh activation function; the latter has three inputs (latent variables and index ambient dimensions) and one output (the value of one state), with hidden layer of 20 neurons and tanh activation function. The LDNet is designed to be time- integrated by the forward Euler method. The model is trained using first Adam optimizer with learning rate of 10- over 200 iterations, and then BFGS optimizer for 1300-4000 iterations, until the loss drops to below x 10 10-9. However, is worth noting that we failed to train LDNet using the same size of training dataset used by GMKRR, and the reported prediction results based on the training of twice the size of the dataset used in training GMKRR.
Besides the predictions of the last periods in Fig. 1 and the entire time in Fig. the test RMSEs are provided in Table[3 While the RMSEs of Nuler+NC and NODE are comparable for the most of D's, the performance of the Euler+NC method clearly better than the NODE when D 0.1. For this latter case, NODE fails to maintain the geometrical constraint over long time prediction. The LDNet behaves in an irregular manner too: it achieves the best accuracy for D 0, 0.3, 0.5, but fails to predict the true trajectory after about 50 model time unit for D = 0.1. Furthermore, that when only one latent variable used (i.e., the same as the intrinsic dimension) LDNet could not produce any useful predictions. Furthermore, the average training and prediction costs for all the methods are given in Table While the prediction costs of the methods are comparable, the training costs of NODE and L.NNet are 3-4 orders of magnitude higher than the GMKRR methods.
note
Effect of kernels. Besides the results in Fig. the effect of kernels on the GMKRR model is also examined. We consider the most challenging case of D 0.5, and four kernels:
SE kernel: p(x,x') =xpp=]xx>x>
2. Matérn 5/2 (M52): p(x,))))(xxxx=)))___)))))))))))) (1+ /5d 5/3d (0=0::/))=(0155555/saaaeeeee////
hacpp::/ggthhu...ccmr/gg//////rr ffeq
24

---

Published as a conference paper at ICLR 2025
Table 3: RMSE values for different cases
Case (D)
Nuler+NC
RK2
NeuralODE
LDNet
Case 1 (D = 0.0)
0.0184
1.5178
0.4372
0.0001
Case 2 (D 0.1)
0.0333
1.9530
1.5649
2.8000
Case 3 (D = 0.3)
0.0224
1.6325
0.2349
0.0002
Case 4 (D 0.5)
0.3842
2.0433
0.4450
0.0157
Table Average training and prediction times for the 1D example.
Method
Training Time (s)
Prediction Time (s)
Euler+NC
5. 5.9
0.1077
RK2 Neural ODE
3.363
580.9
1.003
LDNet
84.81
2.574
3. Matérn 3/2 (M32): p(x,x') = (1 exp)x(vvdddd
Matérn 1/2 (M12): p(x,x') = exp(-d)
where d = d=|xxxx|11
After a cross- cross-validation similar to the one outlined at the beginning of this section, the hyper- parameters for different kernels are listed in Table Unlike the case in Fig, [ the bandwidth parameter for each kernel also obtained by cross- vross-valiattio,, to ensure fair comparison among the kernels. The prediction is shown in Fig. and the RMSE's are listed in Table Clearly, there no significant differences in terms of prediction skills among the different kernels; this confirms that the framework is not restrictive only the SE kernel. This has smooth vector fields. However, we do note that is very likely that there are some examples of certain regularity where it is advantageous to employ the Matérn kernels.
expected since this problem
Table Hyperparameters and RMSE's for models with different kernels.
Kernel
SE Kernel
Order of GMLS for NC, timation,
Tangent Space Es-
RMSE
Band- width,
SE
1.90
4
0.4068
M52
0.28
4
0.5861
M32
0.56
0.7860
M12
0.65
3
5
0.5559
Effect of accuracy of tangent space. Among the examples considered this study, the example is the only one where all the information is known exactly; this includes the manifold parameterization and tangent space as well the dynamics. Hence, this example serves good case for the numerical verification of the convergence bound for the Euler+NC method in which repeated below,
Theorem 3.1
log(N) 5 N
V Eu [Ixx
||x
& I c+Ar
where Co eLT
The verification is designed as follows. First, to avoid the effect of c, i.e., the error in estimated dynamics, the true dynamics in ambient space for the 1D example is used. In this case, the error bound should be O((log(N)/N) /N where the constant depends time step size At = 1 and a short prediction horizon T 10-4, this to avoid the dominance of errors due time discretization and constant C the bound. Subsequently, series of datasets
At. Next, choose small
25

---

Published as a conference paper at ICLR 2025
SE
M32
Truth
M52
M12
-1
95
100
105
110
115
120
125
Figure 6: Model predictions using different kernels.
of size N; 160 i= 1,2,3,4 generated. For each N₁ data points sampled from the manifold, and the dataset is used for the tangent space estimation and the normal correction in the Euler+NC method. Lastly, 100 test trajectories are considered, whose initial conditions are randomly sampled on the manifold. The RMSE of prediction errors at the last step is recorded for each N₁. The convergence plot shown in Fig. which shows the cases 2,3. The two cases closely follow the theoretical convergence rates.
10-5
- ==2
(log(N)NN>2
- f-3
(log(NYNN>33
10-6
RMSE
10-7
10-8
320
640
1280
2560
5120
N
Figure Convergence of Euler+NC method given true dynamics.
E.2 A 2D CAVITY FLOW PROBLEM
Dataset. The 2D cavity flow dataset is publicly available from PanINyySDDy 2024 and is generated from one Direct Numerical Simulation (DNS) solution of a 2D laminar cavity flow problem. According to Callaham al. (2022), the original DNS solution has 30000 steps, with step size of Ar= 0.01s. The first 4000 steps are used to perform a linear dimension reduction by Proper Orthogonal Decomposition (POD). The leading 64 POD modes are retained and all the 30000 steps of solutions are projected onto the 64 POD coordinates. Hence, the available dataset consists of one single trajectory of 30000 steps, and the ambient dimension 64.
26

---

Published as a conference paper at ICLR 2025
For fairness in comparison, the SINDy, M My, and GMKRR models are trained using the first 4000 data points. For testing, the remaining 26000 data points are split into 13 trajectories of 2000 steps.
tlyperparameters of the models. First, the setup of the GMKRR algorithm is discussed. The intrinsic dimension of the dataset is estimated to be 2.38 using kernel-based algorithm (see Appendix() ; this indicates that the dataset approximately lives on wwo-dimensional manifold, and hence d 2 is chosen. The rest of the hyperparameters in GMKRR are chosen via validation procedure similar to that of the ID example. The kernel bandwidth y = 10, polynomial order of GMLS for NC ( 3, and the tangent space estimation is done by local SVD and denoted s 0.
The SINDy algorithm used as baseline case. Due to the limitation on the curse dimension- ality, the SINDy model is only trained to predict the first POD coordinates, using third-order polynomial features.
The M-SINDy algorithm is an existing approach proposed in Callaham et al. (002)) to learn dynamics on manifold. This algorithm first applied Dynamic Mode Decomposition (DMD) to represent the original 64 o4-dimensional dynamics by cemporally-coheree coordinates z(t), x(t) = Uz(r). Then 4 DMD coordinates, denoted 2, are manually selected by physical intuition to represent the manifold which the data points reside. Subsequently, standard SINDy with third-order polynomial feature is applied to learn the dynamics à = f(2), while another fifth-order polynomial model g(²) is trained to estimate the full DMD coordinates using& 2. In the prediction stage, the ODE 2 = f(2) is solved to generate the dynamics in 2, then the POD coordinates are recovered by DMD modes Ug(2). In this paper, the open-source implementation PySINDy (0244) of M-SINDy for the cavity flow problem is used to generate the results for comparison.
modes U and DMD
0.04
0.04
0.02
0.02
Mode
-0.02
-0.02
-0.04
-0.04
19.2
DNS
M:S/NDy
SINDy
CMKRR (Ours)
0.02
0.02
0.02
Mode
0.00
0.00
-0.01
-0.02
-0.03
0.0
0.2
0.4
0.6
0.8
1.0
39.0
19.2
19.4
19.6
19.8
20.0
Figure& 8: Comparison of predictions for a representative trajectory.
Additional results. The predicted dynamics from SINDy, M M-SINDy, and GMKRR are compared for one typical trajectory Fig. g which shows the 1st and 3rd POD coordinates the first second and the last second of test trajectory. The Mode I results show typical characteristics of the three methods. The SINDy starts to deviate from the truth from the beginning, as does not satisfy the manifold constraints. Both M- M-SINDy and GMKRR predict the trajectory of the solutions accurately in the beginning, but over time a discernible phase shift observed the M-SINDy predictions. The Mode 3 results show that, using the same amount of training data (4000 points), GMKRR captures the oscillation amplitudes much better than M- M-SINDy.
The convergence study of GMKRR is performed on the number of samples, N. The correct functioning of GMKRR algorithm requires the training dataset to capture the entire manifold, yet for the quasi-steady cavity dynamics, the complete coverage of the manifold requires over 3000
27

---

Published as a conference paper at ICLR 2025
6.0e-02
Test errer
Reference, order-1.05
2.9e-02
Nax RASE
6.4e-03
500
1000
4000
Figure 9: Convergence study of GMKRR: Max RMSE as function of the number of training samples N.
steps. Hence in the convergence study the following strategy is employed. The original training dataset of 4000 data points are split into 100 shorter trajectories of 40 steps. dataset of N 100K data points is formed by choosing the first K steps of each shorter trajectory. For convergence study, a series of datasets are generated with K 10, 15, 20, 30, and for each model the average егтог of the test trajectories is recorded. The results are shown in Fig. 9 indicating RMSE- O(N }); this convergence rate aligns with standard kernel regression methods.
E.3
KURAMOTO-SIVASHINGG DYNAMICS
Datasets.
The KS dynamics is governed by the following partial differential equation (PDE),
U, + uu, Uxx + Vilxxx= = 0, x E [0, L), 2 0.
Here we consider periodic boundary condition u(r,0) = u(t, L) and a tunable parameter v. When L 2x and given appropriate initial condition =(0,x),
discretized grid of 64 points, one
obtains beating dynamics for v 337
and travelling dynamics for = 4 (Rowley & Marsden 2000 B7
In Floryan & Graham 2022 both the beating and travelling dynamics are considered; the models are trained using short time intervals of training data (e.g.. 1-2 periods oscillation) and predicted the dynamics over long time intervals (e.g.,
10 periods of oscillation). However, in their open- source repository, only the training datasets are provided. To quantitatively assess the long-tern term prediction errors, we have to regenerate the complete long-term datasets. Furthermore, since the numerical methods employed by Fllorya & Graham 2022) and us differ, we take the following procedure to ensure the consistency of datasets.
Our numerical solution of the KS equation is obtained using the standard Exponential Time Differencing method. Next,
1. Beating dynamics: We choose step size of 0.0001, and the initial condition to be the last step in CANDyMan's training dataset of beating dynamics. The simulation run for 75000 steps, and the initial 25000 steps of transients are discarded. The remaining 50000 steps are subsampled with a step size of 0.01, so that eventually 500 steps are saved as the dataset. The period of beating dynamics is approximately 45 steps.
2. Travelling dynamics: We choose a step size of 0.001, and take the same strategy as above for the initial condition. The simulation is run for 180000 steps, and the initial 40000 steps of transients are discarded. With a sampling step size of 0.01, 14000 steps are saved the dataset. The long period of travelling dynamics is approximately 5000 steps.
Both the GMKRR and CANDyMan algorithms are tested on these newly-generater datasets. For CANDyMan, since the physical behaviors of the new datasets are the same as the original ones, we employed the identical model architectures and all hyperparameters (e.g., number neurons and training epochs) from the original paper. The hyperparameters of GMKRR models are found through a cross-vee validation procedure and listed in Table Note that here the initial guess of the bandwidth parameter determined using y defined in App. B.2
28

---

Published as a conference paper at ICLR 2025
Beating dynamics. The GMKRR and CANDyMan models are trained using the first 100 steps of dataset, which correspond to approximately 2 periods of oscillations. The prediction is computed by taking the first step of the dataset as the initial condition and performing time stepping for 500 steps, i.e., approximately 11 periods.
The predicted spatiotemporal dynamics are shown in Fig. 10 visually similar results, as shown in Fig. [0a] GMKRR does show one order of magnitude lower error than CANDyMan. The sources of error in CANDyMan are explained via the comparison of the time histories of u at two spatial locations; in Fig. Tob the last two periods of time history are shown. Clearly, CANDyMan shows some phase error after 11 periods of prediction. Furthermore, in the second panel, while the amplitude of oscillation is low, CANDyMan shows non- physiccll tighly-oscillatoryy response; this is attributed to the errors in the autoencoder of CANDyMan.
While both models seem produce
Travelling dynamics. As explained in Sec. 5.2 five models are considered for this case: CANDY- Man, NODE, LDNet, GMKRR-Full Full, and GMKRR- FFT. The GMKRR-Ful model is trained using every other step in the first 6000 steps of data (so total of 3000 steps) without assuming any structure in the model. The subsampling of every other step is to reduce the computational cost of the kernel regression. The CANDyMan, NODE, L.Nett, and GMKRR- KFT-FFF employ the following special data preprocessing step. Denote the data sample at time f as u(f) E @ The preprocessing starts with a discrete Fourier transform to decompose the data u(r)= = Re(û()) exp(ic(t))) into a standing wave (()) and its phase shift (1). The standing wave and phase dynamics capture the fast and slow timescales, respectively, and both have an intrinsic dimension of 1. After the preprocessing step, the four models only needed the first 100 steps to train for the standing wave dynamics u ú= SNN (ū) as well the change in phase per time step AΦ- = gN)(
The hyperparameters for CANDyMan are the same as the original paper, where its atlas-of-charts approach applied to fNN and gNN is a vanilla fully-connected NN (for each chart). For the NODE, fNN consists of 64 inputs and 64 outputs, with two hidden layers of 128 neurons each and tanh activation function, and the rest the same the 1D case, except that 2000 iterations were used in training: gNN consists of one hidden layer of 32 neurons and tanh activation function, and trained by Adam optimizer with learning rate 10 3. The LDNet is the same as NODE, except the fNv is replaced by the latent dynamics. The dynamics network has 64 inputs and two outputs (latent variables), with two hidden layers of 8 neurons each and tanh activation function; the reconstruction network three inputs (latent variables and x-coordinate) and one output (the value of one state), with two hidden layers of 8 neurons each and tanh activation function. The rest is the same the ID case, except that 6000 iterations the same dataset as other models. Note that LDNet would fail to produce reasonable predictions if only one latent variable (i.e., the intrinsic dimension) is used.
used in training. In this case L.DNet use
Figure lla shows the initial short period of dynamics. All models perform reasonably well, and all captures the short time-scale with low errors. Yet, the GMKRR models either outperform or are on par with other baselines in terms of prediction accuracy. Figure] shows the long-term dynamics. There is clear difference in the error level across the three models, that is attributed to the mismatch in the long time-scale. GMKRR- FFT almost captures exactly the long time scale; this results the lowest error, that is almost one order of magnitude lower than all other baselines, even though these models leverage the same prior knowledge. The comparison highlights that, while the models that do not account for geometrical constraints (i.e., NODE and LDNet) can capture short-term dynamics well, their accuracy quickly deteriorate over time. The GMKRR outperforms other models, which is attributed to its capability to model the dynamics using exactly the intrinsic dimension.
E.4 REACTION- DIFFUSION DYNAMICS
The reaction-difffssion dynamics is governed by the following coupled 2D PDEs,
[1 - 1111(42+4 + + flu² +2)++ (Uxx + Uyy)
V; -p(u² v² (n² + d2(Vxx Uyy)
on square domain (x, E 10) 1 10) with homogeneous Neumann boundary conditions
and parameters d = d2 = 0.1, p= 1. We follow the same numerical solution method as in CANDyMan by directly transcribing the MATsAd- solver in its public repository into Python
29

---

Published as a conference paper at ICLR 2025
Truth
CANDyMan
GMKRR
6
12
2
-3
-6
-9
0
-12
RMSE 2.744e-01
RMSE 2.549e-02
-3
-6
-9
-12
(a) Full spatiotemporal dynamics.
x/L=0.39
X/L=0.78
8.78
-7
8.76
8.74
8.72
-9
Truth
CANDyMan GMKAR
4.0
4.5
5.0
4.0
4.5
5.0
(b) Dynamics at select spatial coordinates.
Figure 10: The KS beating dynamics.
30

---

Published as a conference paper at ICLR 2025
Truth
GMKRR Full
GMKRR FFT
15
RMSE 2.519e-02
RMSE 5.787e-03
CANDyMan
NODE
LDNet
12
2
9.909e-02
RMSE 1.251e-02
RMSE 3. 3.669e-01
15
12
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
1.5
(a) Initial spatiotemporal dynamics showing the short time-scale.
Truth
GMKRR Full
GMKRR FFT
Ch 0001
RMSÉ 1.604e400
RMSE 3.671e-01
6
-9
CANDyMan
NODE
LDNet
12
RMSE 2.391e+00
RMSE 6.690e+000
RMSÉ 7.600e
15
12
100
(b) Full spatiotemporal dynamics showing the long time-scale.
Figure 11: The KS travelling dynamics.
31

---

Published as a conference paper at ICLR 2025
Table Model parameters of GMKRR for the CANDyMan cases
CASE
MODEL
Intrinsic Dimen- sion, d
SE Kernel
Order of
Tangent
Band- width, y
GMLS for
Space Es-
NC, l
timation,
KS Beating
GMKRR
440.163
4
0
KS
GMKRR-Full
3
98.0816
3
0
Travelling
GMKRR-FFT
1
61.47
4
0
GMKRR-Full
1
31000.6
4
0
Reaction-
EMKRR--CA-TIE
31319.5
Diffusion
GMKRR-PCA-T4
1
31319.5
4
4
using NumPy/SciPy. The domain is discretized into uniform 101 101 grid, the solution one time step is of dimension 20402. A step size of 0.05 is used, and 10200 steps are simulated; the first 5200 steps are discarded to avoid transients. A period of response approximately 72 steps, and hence the full dataset contains nearly 70 periods. All the models are trained using the first 200 steps and tested on the full 5000 steps of data.
Figures 12a and 12b show results at the last time step. While visually all spirals seem similar, CAN does have errors that are two orders of magnitude larger than GMKRR methods. In addition, the high-order tangent space estimation shows further and dramatic improvement over the other GMKRR models. The large error in CANDyMan is explained by Fig. 12c The Man model does not accurately capture the period, and hence results in accumulation of error over time.
32

---

Published as a conference paper at ICLR 2025
Truth
CANDyMan
GMKRR-FUll
GMKRR-PCA-TO
GMKRR-PCC-T4
1.00
0.75
5
0.50
0
0.00
-5
-0.50
-0.75
-10
RMSE 3.690e-01
RMSE 1.398e-03
RMSE 1.
RMSE 9.310e-05
0.75
0.50
0.00
-0.50
-1.00
10
0
10-10
10-10
10-10
10
(a) Solution at = 250, component
Truth
CANDyMan
GMKRR-Full
GMKRR-PCA-TO
GMKRR-PCA-T4
0.75
0.25
0-
0.00
-0.25
-5
-.50
-10
-1.00
RMSE 3.690e-01
RMSE 1.398e-03
RMSE 1.399e-03
RMSE 9.310e-05
1.00
0.75
0.50
0.25
00.25
-0.75
-1.00
10
10-10
10-10
10-10
x
(b) Solution at = 250, component
1.0
0.5
Tutth CANDyMen SMKRR-Fall SMKRR-PCA POATTU
-0.5
GMKRR-PCA
-1.0
35
45
50
130
135
140
1.45
150
230
235
240
245
250
(c) The u dynamics at (x,y) (10;- all GMKRR models overlap.
Figure 12: Reaction coaction-diffusion dynamics.
33


Reviews:
[{"review_id": "BBlNtPlrmO", "reviewer": "Reviewer_pB4e", "summary": "The authors propose a kernel method for learning and integrating ordinary differential equations on manifolds. The main contribution is the method for learning ODEs, called Geometrically constrained Multivariate Kernel Ridge Regression Method (GMKRR), a method that relies on the theory of Reproducing Kernel Hilbert Spaces. As for integrating, the author propose a new Euler method with a Normal direction correction.", "questions": "Instead of the Euler-NC method, have you considered a second order Euler method for the ODE, by differentiating the kernel operator?", "limitations": "", "rating": 10, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The paper is very well written. The topic relies heavily on results from operator kernels, but I found that the introduction to the topic by the authors was adequate.\nThe contributions of the authors are novel and have promising applications in data-driven discovery of latent dynamics.\nThe experimental evaluation is good.", "weaknesses": "The method involves solving very large systems of equations at each step, and thus is computationally expensive. Nevertheless, the initial training time, (if existent, for instance when solving the system with a Cholesky decomposition) is much less compared with a PINN.\nFurthermore, the timings are superior when compared with other data-driven methods, such as CANDyMan.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "ZlfDFPw5Jf", "reviewer": "Reviewer_u6B3", "summary": "Authors propose a new method for solving ODEs on manifolds in the multidimensional case which is based on the geometrically constrained operator-valued kernel and a geometry-preserving ODE solver. The efficiency and advantages of the method are demonstrated experimentally for a number of multidimensional problems and in comparison with alternative approaches.", "questions": "1. In the abstract you note that \"machine learning methods ... struggle with ... lack of interpretability\". It is not entirely clear from the text of the manuscript whether your method solves this particular problem.\n\n2. Can you comment on whether NC Euler method (from sec. 3.2) is new and what alternative approaches exist for predicting dynamics on complex manifolds (if as you note RK2 method behaves badly here)?\n\n3. Section \"Related work\" lacks a description of the essence of CANDyMan approach.\n\n4. Is it possible to consider more baselines, e.g., method from the work \"Neural Manifold Ordinary Differential Equations\", described in \"Related work\"?\n\n5. There are probably a few typos. Line 074: \"We supplement the papers with...\" (probably should be singular, i.e., \"paper\"); Line 446 and 475: \"submanifold manifold\" (probably word \"manifold\" is redundant here); Line 449: \"the cost ... faster\" (maybe cost is lower or training is faster).", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The presentation of the material in the text of the manuscript is consistent and neat.\n\n2. The problem of solving such a class of equations is relevant for a wide range of applications.\n\n3. It has been empirically demonstrated that the proposed approach significantly outperforms known baselines (SINDy methods, CANDyMan) in both accuracy and speed on several classes of problems.", "weaknesses": "Overall, the work appears to be of high quality and relevant, however I would like to clarify the issues listed below in the \"Questions\" section.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "HEjwa3OxtZ", "reviewer": "Reviewer_7GGN", "summary": "This paper definitely proposes an interesting method to learn the dynamics of a dynamical system that lives on a low-dimensional manifold based on point cloud-type data. The \"geometrical constraint\" refers to that the bases for building the kernel regression are projected to (an approximate) tangent space of the observations. The procedure to obtain the kernel regression's solution is modified to leverage of the (estimated) intrinsic dimension $d$. Then a \"geometry-preserving\" integrator is applied to get the trajectory based on the $\\mathbf{f}_{\\epsilon}$ learned in the kernel regression.", "questions": "N/A", "limitations": "", "rating": 5, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- Efficient compared with the baseline.", "weaknesses": "- Why just 1 baseline?\n- Line 136: the writing on why the kernel regression opts for the L2 regularization and not the L1 regularization SINDy used is not very convincing.                     \n- The paper only compares the end results with SINDy-based method. I suggest doing a more lateral comparison within the method itself.\n  - Why not comparing the performance of different kernels such as Matern?\n  - One of the most important aspect in this method is to compute the projection (matrix) onto the tangent space. Yet, this part is shoveled into the appendix. For example, the robustness is never mentioned, how the sensitive the approximation to the projection matrix affects the accuracy of the whole method? I say this because for multiple toy models in this paper, the representation of the manifold is known, then the analytic tangential projection matrix can be formed, which sets a perfect playground for such a study.\n  - Since one has $\\mathbf{f}_{\\epsilon}$ once trained, why not comparing Euler+NC with a symplectic integrator (such as semi-implicit which shares similar spirit with this correction method)?\n\n\n### Typos\n- line 191: I am not sure if it is a typo but the squared exponential should not look like that, that is the Laplace kernel.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "pez7n0ia80", "reviewer": "Reviewer_zhfF", "summary": "This paper overcomes some of the limitations of SiNDY and CandyMan to learn ODEs with solutions represented manifolds. The authors show that by enforcing geometric constraints, they can improve efficiency in training as well as accuracy compared to other methods on toy examples by ensuring solutions lie close to a lower dimensional submanifold that can express the data.", "questions": "Are there specific types of out-of-distribution scenarios where the method is expected to perform well or poorly?\n\nCan this method be used on domains or solutions that do not have or allow periodic behaviors that reduce the intrinsic dimensionality of the original solution manifold? It would be beneficial if the authors could provide a discussion on how the proposed method might be adapted for non-periodic systems.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper is well written and provides detailed theoretical justification for their framework with clear numerical experiments highlighting the proposed method’s superiority to state-of-the-art works in this area. They benchmark with several new methods that are suitable for the problems they focus on and outperform in both training efficiency and accuracy.", "weaknesses": "The authors mention that their approach overcomes limitations of other algorithms in regard to making predictions on out of distribution data, yet they don’t detail the extent of which their method can practically handle problems where the manifold doesn’t exhibit properties or yield constraints that allow you to represent your solutions on tangent planes of low dimensional submanifolds. It would be beneficial to provide specific examples of manifold types or problem domains where their method might struggle.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "NY6LqUQtwR", "reviewer": "Reviewer_YaJd", "summary": "The paper presents an approach for learning ODEs on manifolds using kernel methods such that learnt ODE solutions are constrained to lie on the learnt manifold. They also present a method for solving ODEs on manifolds more accurately than standard RK based approaches in terms of how closely solutions lie on said manifold.", "questions": "1. How do you view your approach as being more interpretable than neural-network based methods?\n2. How do you suspect you suspect your approach will scale for problems with higher dimensional submanifolds?\n3. How do suspect your method will perform in the presence of significant observation noise given that you are required to estimate the manifold and the time derivatives on the manifold from data?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "- Learning ODEs constrained to evolve on manifolds using linear methods is a beautiful idea.\n- The geometry-preserving time integrator is well-motivated and could have impacts well beyond learning differential equations from data. \n- The numerical studies convincingly demonstrate the proposed approach is better than the SINDy-style approaches they compare to.", "weaknesses": "- While theoretically interesting, I think the authors could have done a better job of motivating their approach in the introduction from a practical perspective. In particular, it wasn't clear to me that the standard approach of learning latent differential equations, is limited in practice by the fact that learnt differential equations will deviate from the learnt manifold. Since the standard approach is is so widely used, providing concrete references / studies which demonstrate these issues in practice would have been helpful for motivating the work. \n\n- I think it's arguable that the one of the biggest challenges with SINDy-style algorithms is that they struggle in the presence of noise, for example:\n\n    - Messenger, Daniel A., and David M. Bortz. \"Weak SINDy for partial differential equations.\" Journal of Computational Physics 443 (2021): 110525\n\n    - Fasel, Urban, et al. \"Ensemble-SINDy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control.\" Proceedings of the Royal Society A 478.2260 (2022): 20210904.\n\n    Since your approach also requires estimating derivatives using finite differences, it would have been helpful to include some discussion on this issue. \n\n- You mentioned that one of the goals of your work as compared to standard approaches for learning latent differential equations is to learn interpretable models. As I understand it, SINDy-style models are interpretable because they learn sparse approximations to differential equations. Since you relax the sparsity requirement, I struggled to understand how your proposed approach is more interpretable than neural-network based models.\n\n\n*Additional feedback. To be clear, these points are here to help and will not be a part of the decision assessment.* \n\n- When introducing SINDy, I think it would be helpful to include a reference to this work since I think they also proposed to learn sparse dynamics using dictionaries of basis functions in 2011.\n\n    - Wang, W.X., et al. Predicting catastrophes in nonlinear dynamical systems by compressive sensing. Phys. Rev. Lett. 106, 154101 (2011).\n\n- The results from your numerical studies would be made stronger by including comparison to learning latent ODEs with neural networks, for example:\n\n    - Toth, P., Rezende, D. J., Jaegle, A., Racanière, S., Botev, A., & Higgins, I. (2019). \"Hamiltonian generative networks.\" arXiv preprint arXiv:1909.13789.\n    - Chen, Ricky TQ, Yulia Rubanova, and David Duvenaud. \"Latent ODEs for irregularly-sampled time series.\" Advances in Neural Information Processing Systems 32 (2019): 3.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Spotlight)