Nyström-Accelerated Primal LS-SVMs: Breaking the $O(an^3)$ Complexity Bottleneck for Scalable ODEs Learning
Authors: Weikuo Wang
Yue Liao
Huan Luo
Venue: NeurIPS 2025 poster ()
Keywords: Closed-form approximate solution
Nyström-accelerated learning theory
least squares support vector machines (LS-SVMs)
ordinary differential equations (ODEs)
Primary area: theory

Abstract: A major problem of kernel-based methods (e.g., least squares support vector machines, LS-SVMs) for solving linear/nonlinear ordinary differential equations (ODEs) is the prohibitive $O(an^3)$ ($a=1$ for linear ODEs and 27 for nonlinear ODEs) part of their computational complexity with increasing temporal discretization points $n$. We propose a novel Nyström-accelerated LS-SVMs framework that breaks this bottleneck by reformulating ODEs as primal-space constraints. Specifically, we derive for the first time an explicit Nyström-based mapping and its derivatives from one-dimensional temporal discretization points to a higher $m$-dimensional feature space ($1< m\le n$), enabling the learning process to solve linear/nonlinear equation systems with $m$-dependent complexity. Numerical experiments on sixteen benchmark ODEs demonstrate: 1) $10-6000$ times faster computation than classical LS-SVMs and physics-informed neural networks (PINNs), 2) comparable accuracy to LS-SVMs ($<0.13\%$ relative MAE, RMSE, and $\left \| y-\hat{y}  \right \| _{\infty } $difference)  while maximum surpassing PINNs by 72\% in RMSE, and 3) scalability to $n=10^4$ time steps with $m=50$ features. This work establishes a new paradigm for efficient kernel-based ODEs learning without significantly sacrificing the accuracy of the solution.

Nyström-Accelerater Primal LS-SVMs: Breaking the O(an²) Complexity Bottleneck for Scalable ODEs Learning
Weikuo Wang",2, Yue Lia Huan L.
'College of Civil Engineering and Architecture, China Three Gorges University, Yichang. China
Hubei Geological Disaster Prevention and Control Engineering Technology Research Center, Yichang, China College Basic Medical Sciences, China Three Gorges University, Yichang, China "Hubei Key Laboratory of Tumor Mieneenvirommett and Immunothcrapy, Yichang, China {ukwang, liaoyue, hluo)@ctgu. weriggu.ddu...m edu. cn
Abstract
A major problem of kernel-based methods (e.g., least squares support vector chines, LS-SVMs) for solving linear/nonlinear ordinary differential equations (ODEs) is the prohibitive O(an²) (a 1 for linear ODEs and 27 for nonlinear ODEs) part of their computational complexity with increasing temporal discretiza- tion points n. We propose novel Nyström--accelerated LS-SVMs framework that breaks this bottleneck by reformulating ODEs as primal-space constraints. Specifically, we derive for the first time an explicit Nyström- based mapping and its derivatives from one-dimensional temporal discretization points to higher m-dimensional feature space (1 m S n)), enabling the learning process solve linear/nonlinear equation systems with m-dependent complexity. Numer- ical experiments on sixteen benchmark ODEs demonstrate: 1) 10 6000 times faster computation than classical LS-SVMs and thysics-informas neural networks (PINNs), 2) comparable accuracy to LS-SVMs (< 0.13% relative MAE, RMSE, and || ~ xe difference) while maximum surpassing PINNs by 72% in RMSE, and 3) scalability to n 10* time steps with m 50 features. This work establishes a new paradigm for efficient kernel-based ODEs learning without significantly sacrificing the accuracy of the solution.
Introduction
Ordinary differential equations (ODEs) are foundational tools for modeling dynamical systems across scientific domains, including physics, engineering, and biology [ Classical numerical methods. such as Runge-Kurt schemes and finite difference discretizations. have been the comerstone of ODE solving due to their rigorous error analysis and convergence guarantees. However, these methods face limitations in achieving satisfactory accuracy with large time steps and incur high computational costs in long-time time integration, particularly for stiff ODEs 3 Recent advances in machine learning (ML) have introduced data-driven paradigms for solving ODEs, such as Physics-Informer Neural Networks (PINNs) 44, neural ODEs (5), and Gaussian process-based solvers (6. These approaches demonstrate unique advantages in addressing inverse problems, enabling adaptive resolution, and directly incorporating observational data. These capabilities are often challenging for classical techniques.
Among ML-driven methods, kernel-based strategies--prriculaat Least Squares Support Vector Machines (LS-SVMM) 7 - emerged compelling alternatives for solving ODEs. By refor-
*Corresponding authors
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

mulating the ODEs as constrained optimization problem in a reproducing kernel Hilbert space (RKHS), LS-SVMs leverage kernel functions to implicitly capture nonlinear dynamics while ensuring regularization against overfitting. This framework has shown success in solving initial/boundary value problems g and parameter estimation tasks [9, benefiting from the inherent flexibility and theoretical soundness of kernel methods. However, critical limitation persists: the computational complexity of kernel-bacc -baed ODE solvers scales as O(n + p))) when solving linear ODEs, and as O((3 + p 2)º) per Newton iteration in the case of nonlinear ODEs (8) where rt is the number of temporal discretization points and p denotes the order of the ODE. This O(are) (a = for linear ODEs and 27 for nonlinear ODEs) scaling constitutes major computational bottleneck, rendering traditional LS-SVMs impractical for long-time time simulations or fine- grained temporal resolutions, and thereby undermining their utility in large-scatl acaee scientific applications.
To address this challenge, we propose a new Nyström-acceleran the computational complexity to (((m+ + p) ) for linear ODEs and ( + p+ +n)) per Newton iter- ation for nonlinear ODEs, where m < n represents the number of subsampled landmark points. The key innovation of our approach lies in its primal-domain formulation, which strategically integrates the Nyström method tradionltaa used for low-rank kernel matrix approximations [[0] into the LS-SVMs optimization procedure. Unlike dual-domain kernel methods, our framework leverages the Nyström approximation to construct an explicit Lnite-dimensional feature map, enabling efficient computation of high-order derivatives essential for ODE operators. This not only preserves the theoretical benefits of cernel-based regularization but also facilitates a scalable surrogate model for high-resolution ODE systems. The main contributions include: (1) systematic integration of the Nyström approximation within the LS-SVMs ODE solver, ensuring stability through error- controlled subspace selection; (2) a novel application of the Nyström method to obtain explicit feature mappings and their derivatives, enabling accurate representation of differential operators; and (3) comprehensive numerical validation on sixteen benchmark problems-including stiff and nonlinear vystems-demonstrall
L-SSVMs framework that reduces
order-of-magniuude speedups without sacrificing solution fidelity.
This work bridges the gap between the expressive power of kernel methods and the scalability demands of modern scientific computing. By mitigating the O(an3) complexity barrier, our method unlocks the potential for LS-SVMs to tackle large-scale ODE problems prevalent in multi-physte simulations, biological network modeling, and real-time control systems, where traditional kernel solvers were previously deemed infeasible.
1.1 Related Works
The development of solvers for ordinary differential equations (ODEs) navigates trade-off between computational efficiency, theoretical robustness, and application flexibility. This section reviews the landscape of numerical, neural, and kernel-based methods, outlining the distinct challenges that motivate our work.
Classical and Enhanced Numerical ODE Solvers: Classical solvers, such as Runge- Kutte schemes [ and linear multistep methods [12, form the bedrock of ODE simulation. Their principal strength lies in well-established convergence guarantees and explicit error control through adaptive time-stepping. However, their performance is intrinsically linked to the temporal discretization. Stiff ODE systems usually necessitate exceedingly small time steps to maintain stability, leading to prohibitive computational costs in long-time integration. While probabilistic enhancements 13 and Bayesian filters [14] introduce valuable uncertainty quantification, they often exacerbate computational burdens. Even MI-enhanced controllers [15] 16, 17 primarily optimize existing parameters without overcoming the fundamental constraints of discretization, leaving scalability challenges for high-resolution systems largely unresolved.
Neural Networks for ODE Systems: Neural network- based approaches represent a paradigm shift towards ODEs solving [8 19. Neural ODEs (5) offer a continuous-depth framework that adapts flexibly to irregular time series, and Physics- Informed Neural Networks (PINNs) [4] enable the seamless solution of inverse problems by incorporating physical laws as soft constraints. Despite their flexibility, these methods are prone to convergence issues due to non-convex optimization, often resulting in local minima and unpredictable training outcomes 20]. Furthermore, they can exhibit numerical instability when applied to stiff systems (2 and suffer from spectral bias [22] which impedes the accurate resolution of high-frequency solutions. Although differentiable solvers improve gradient flow, they inherit the iterative cost of the underlying numerical schemes.
2
2

---

Kernel Methods in ODE Solving: Kernel methods, including formulations based on Least-Squares Support Vector Machines (LS-SVMs)
and Gaussian processes 2 24 25,, provide mathe-
matically rigorous alternative. Operating within Reproducing Kernel Hilbert Spaces (RKHS), they offer convex optimization landscapes that guarantee convergence to global minimum, avoiding the pitfalls of non-convex training endemic to neural networks. This framework yields closed-form solutions with strong theoretical error bounds, as evidenced in kernel collocation techniques [26]. The primary limitation, however, is their computational scalability: solving the resulting dense linear system requires matrix inversion with (an") complexity, which becomes prohibitive for large-scale or long-ting time ODE simulations.
Scalability in Kernel-Based ODE Solvers: Accelerating kernel methods has been a central focus in supervised learning, where techniques like the Nyström method [TO] (77 (28) and its advanced variants 22 30 successfully reduce complexity to (m)")) for m < ni landmark points via low- rank matrix approximations. However, these innovations have seen limited translation to the ODE context. Exceptions, such as reduced-rank Kalman filters 31. exploit low-rank structures but remain disconnected from kernelized optimization frameworks like LS-SVMs. Consequently, significant gap persists; no prior work has systematically reduced the ((an)) bottleneck of LS-SVMs for ODE solving while preserving their regularization benefits. Our method addresses this by integrating Nyström approximation directly into the LS-SVMs optimization process.
1.2 Novelty and Contributions
Our work bridges the gap between kernelized ODE
Given on ODE
fome. /U...3' /
solvers' theoretical strengths and their practical scal- ability constraints. In contrast to classical Nyström
("()- A0
Harhh 1- a.ydle.
weethdsswhiich are typically applied in the dual formulation of kernel-base gression problems using the kernel trick [10] our method operates entirely within the primal formula- tion. This provides explicit representation of the nonlinear feature mapping and its high-order deriva- tives, which is specifically tailored to the temporal structure of ODEs and ensures stability in stiff sys- tems (Figure M. This aligns with trends in data- driven scientific computing but specifically targets the under-explored challenge of O(an)) complex- ity in kernel- based ODE solvers. Our framework could be extended to large-scale problems in systems biology, control theory, and multi-agent dynamics, Figure Flowchart of Nyström-accceleraten L.S- where traditional kernel-based ODEs solvers were SVMs for Efficient ODE Learning. previously impractical.
The Legraeggon of the grodene
opeenddyyy
classification and re-
1
N00 optimnally conditos
76
Explict feature
vodel L.ns var.te P...
18 woon purcoo
PV 10.0 f (.)))))) [nr
Salving the primall problem
J = e +'x))x) +
enthedd
closed form appr osimation solution
2 Preliminaries
This section describes the problem statement, a short introduction to LS-SVMs (Appendix A. for ODEs solving is given to highlight the prohibitive (an²) complexity problem considered in this paper. Consider general p-th order linearnonlineaa ODE with the following form:
Cly) = y (p) = [(,, v. V' ... (P--)), I E [41,
where L represents an p-th order linear/nonlinear differential operator depending on the lin- ear/nonlinear function f(t, y, v'.... , -(p-11), [n] is the problem domain, t is the input signal, and (()( (1) denotes the ε-th derivative of y with respect I, E [0.p). The p necessary initial or boundary conditions for solving the above ODE are
IVP: Iyl = Un BVP: BValv] = qas = 0. ......p
(2)
where IV,, are the initial conditions (all constraints are applied at the initial value of the independent variable i.c.. t = t1) and BV,, are the boundary conditions (the constraints are applied at multiple
3

---

values of the independent variable t. typically at the ends of the interval [t1, tn] in which the solution sought). Up and 9N are given scalars.
When the function ftt, y: v. -)(=11) in (I is linear, the p-th order ODE to be solved is linear. The optimization problem based on LS-SVMs for learning the p-th order linear ODE can be formulated by the method in Ref [8] Its solution requires solving a system of (n + p) linear equations with ([ + p) computational complexity. When the function is nonlinear, the p-th order ODE is nonlinear and written in the following form (initial value problems for illustration):
y ( vt. ... ))
y(11) = Vo, ypp-1)
Up-1, & S tx
(3)
Assume the approximate solution: (t)) w" (t() b. Additional unknowns y, are introduced to keep the constraints linear in w. This yields the following nonlinear optimization problem:
200 w w +222200 2
minimize W.b.
s.t. yd w² (1))) (t,) f (t (tr Hi. . VVi 1) + +ei ==2.....!
v(41) = w 40 (41) + U1 (1)) wn (11-11)(1) = ue, 2= y: = w (1,) + +b+ E.. i= 2, nn.....
By using classical kernel trick, Eq. can be effectively solved. For detailed mathematical derivation, Appendix A.2 The nonlinear system, which consists of (3n+p-- equations with (3n (3n+pp- unknowns, is solved by Newton's 's method (32.. For each Newton iteration, solving for unknowns requires the computational complexity of O((3n - p 2) In cases of long-time interval with more than 10a time steps, the prohibitive O(27n²) part of the computational complexity causes the expensive computations even may stop the run.
3 Proposed Nyström--cceleratee LS-SVMs (NLS-SVMs)
We consider the explicit model. j((t) = w (t) + b. as a closed-form approximation solution to the ODE (i.e., [D). A key innovation of the NLS-SVMs approach is its application of the Nyström method to directly approximate the high-dimensional nonlinear feature map (t) and its higher-order derivatives, enabling the direct substitution of the approximate solution ((t) and its derivatives into the governing equation. Therefore, the NLS-SVMs method bypasses the kernel trick, allowing the model parameters be efficiently estimated in primal form, with computational complexity dependent only on the dimension of the approximated features. Leveraging Mercer's theorem, the derivatives of the feature map pt() can be analytically expressed using derivatives of the kernel function. To this end, we define the following differential operator:
[ (t) = (2, v) w= [92k,
V. v) v)ll
gota K K(u,e) 0u=
(5)
where [0])y denotes the (i, (i.j))th entry of matrix D Given long-time interval He the Nyström method, as defined in Eq. 34
B4) (Appendix (A.3). is employed to construct an explicit finite- dimensional feature mapping (tt) E R with (m < n). This allows the model parameters w Rm
and E R to efficiently estimated in the primal space, circumventing the computational burden associated with the dual formulation. Throughout this paper, we adopt the radial basis function (RBF) kernel K (u, v) = erp( (u 2)² /a²). Using the differential relations established in Eq.(5). the first- and second-order derivatives of the kernel function are given by:
Vi(K v) 2(u ( v) r²
K(u, u), (2(KK (u,v)) 4(u v)²
2 o²
K K(u,v)
(6)
By extending Eq. (34 the p-th order derivative of the feature map (pt) can be explicitly written as:
ok (p) (t.) = - ( (t,) = vm c
"m K (tk, , t.)
(7)
These explicit expressions for the feature map and its derivatives enable the direct numerical treatment of y-th order linear and nonlinear ODEs introduced in Ref B and Eq. entirely within the primal framework. As a result, the dominant computational complexity substantially reduced from O(a3)) to (m²) for linear ODEs and 3(2)) for nonlinear ODEs, where m <

---

3.1 NLS-SVMs for learning p-th Order Linear ODE
Consider the general p-th order linear ODE with initial value problem (IVP) in
8
p ( Σ fx(t)//o ((t)))(-))))) r(t), E tal y(t1) vi, (
Uk:
k- ......
(8)
The approximate solution can be obtained by solving the following optimization problem:
minimize wa
2 -ww +22224
2 Ei=
(9)
s.t. (t,) = w" (p( (ti) = -- (1(
+ (t.) + Eiss i= 2, 2.....
y(t1)= = =" (t1) + +b= y1. (( (t)== f= 2
Lemma 1: Given differentiable positive definite kernel function k R R IR landmark points subsampled from dataset of time points, where m < n. and regularization parameter y Rt, the solution optimization problem Eq. is obtained by solving the following primal problem.
A ( (m) (1-((=1)
-f
(fp 1)
8 (t1)
(0(n-1)
1) (t1))
(p))) f1 ((=1))
-1 (fp 1)
1
(t1)
0
0
((0-1)(t1)
where E R is identity matrix; & (e((2) ..... (T
R(--1)xm w =
E R(n-1)xm (p)
1)xm
Ar
(1())) (12) (L,) R(n-l)xm f. (12))) ..... (Ln))² E
Vr
@
A1
t1
R" fp (fo(t))..... fo(tn)) E R"
: = [1, .............................. E R 1.
40 (t1) E R (1) (th) E R Rlxm
A = 7( (m()) fi (x-1)
X
Up
RI E mmx(a-1) VV = f1-f,11 E R (m gr
(10)
(r((2), rtt)) E 1=-l: e gmx1 See Appendix A.4 for proof.
3.2 NLS-SVMs for Learning p-th Order Nonlinear ODE
We establish a new paradigm based on NLS-SVMs for solving Eq.
Lemma 2: Given differentiable positive definite kernel function K R R m landmark points subsampled from dataset time points, where m & n, and regularization parameter
E I the solution to optimization problem Eq. (4) is given by solving the following primal problem.
ab
0h
0Y.
w
Yi Y
ØY
0Y
0Y2
b
dY:
BA
dYa
8Y:
Å1
Y3
(11)
Ow
db
dy
By.
Xp
Ypt2
0Y
OY
0Y,
96
0Y,43
aY
dy.
3/
Y p+a
where Y E R Y2 E 1x11, E gixl. .....Y,++.. E 1
[(n-1)x1 The non-
p+3
linear system, which consists of equations with unknowns (w, b, A solved by Newton's method. While more efficient iterative schemes 33
Ap. yi
E R (m++++)) is
34
35 exist, their
exploration falls outside the scope of this paper. It worth noting that when learning a class
of p-th order nonlinear ODEs, the changes in the ODEs primarily affect the (T))())(L
[(4,, Ve. , (p the Lagrangian loss function. Therefore, the Jacobian matrix needs to be updated based on the derivation of the four functions
(WTTp()(L)) "
f(44, Me, V6. "....."",")))
Wt. ")) function, indicated by
oY 0Y, 0Y and 2Y See Appendix A.5 for proof. Oy:
de Oy.

---

Convergence analysis: Since the tunable parameters (i.e., regularization parameter and kernel parameter o²) in the nonlinear system Eq.(1) influence the structure of the Jacobian matrix, we theoretically analyze their impact on the convergence behavior of the Newton solver.
Theorem: Consider the system of nonlinear equations in Eq 11 represented by F(a;; ; = where F' : IR-VV + R (N m + n + p) depends on a kernel parameter o² and regularization parameter y. Let Jr(æ; o² 7) denote its Jacobian matrix. Assume solution x* exists such that
Fa=";"" 7) = 0. The Newton iteration is given by:
Fk+1 = Zy - [Jr(# - FFxx;
(12)
Suppose the following conditions hold in a ball B(x",6) of radius > around x::
1. Lipschitz continuity of the Jacobian: There exists a constant L(o², 7) > 0 such that
|||x(æ; o² Jo(y; 7))l S L(o", ,))l=> yll, Vx, y E ((x', (13)
2. Bounded inverse of the Jucobian: The Jacobian matrix (a o², 7) is non-singular for all a E Bix", 5), and there exists a constant 3(a², ,)) > 0 such that
|| 02; , -|| s 3(02 Ve E B(x",6)
(14)
Then, for any initial guess To ((x", ,8) satisfying
1 5, K where K= 3= 3(0²: [
ax x || < min
(15)
the iteration converges quadratically to x=", with the error Ft x² satisfying:
||e+i|l < K||ex||
(16)
Proof: According Eq (2). the error the next step Ck+1 1*+1 x= can be expressed as (parameters 7 and x² are omitted for clarity):
Ck-1 ex (J=( = [Jr( (Jr(an)) (F(xx) F(x))))
(17)
The term F((x)) F(x") is estimated using Taylor expansion with F(x*) = 0:
Jr J= + tex) F(xx) F(x") S''
(18)
Substituting this back yields:
ex+1 = [Jf(EN)] ( (1(eett)) - JF(Ø" tex)]exdt
(19)
Taking norms and applying the Lipschitz condition (3) to the integrand gives:
(xxx) - J + tex) S L|eex tex || = L(1 +)]lee||
(20)
The integral f2(( t)dt = 1/2. Then, applying the bound on the inverse of the Jacobian (Eq. 14)).
[[]J+(x))] S 3, we obtain the final result:
|ec+1 || S -/2
3. 11e K|lexil?
(21)
Implications: The convergence rate is governed by the product L)) which is jointly influenced by the parameters y and o². Here, the regularization parameter y ensures the non- singularity of the Jacobian, while the kernel parameter " controls the smoothness of F. A larger o² generally reduces the Lipschitz constant L, thereby promoting convergence the potential cost of oversmoothing. Consequently, a careful balance between o² and is critical to achieving optimal trade-off among convergence speed, solution accuracy, and numerical stability.
Computational complexity reduction: The computational complexity of the proposed NLS-SVMs method, as indicated in (1) and Eq.(IT). substantially reduced compared traditional LS- SVMs ODE solvers. It scales with the dimension m (m < n)) of the nigh-dimensional feature (t,)). For linear ODEs, it decreases from O((n + p)²) O(I + )) while for nonlinear ODEs, the per-New Newton-iteration complexity is lowered from O(33 + 2) to O((m + + n))3). Empirical analysis demonstrates that the proposed approach significantly enhances computational efficiency while preserving prediction accuracy. Moreover, the memory requirement for linear ODEs decreases from O(n + p)²) to ((m + 2)² and for nonlinear ODEs decreases from O(33 + p 2)²) to O((m + + 2) Appendix A.6 covers the initial value problem (IVP) for the first-order linear ODE, Appendix A.7 presents the IVP and boundary value problem (BVP) for the second-order linear ODE, and Appendix A.8 provides the derivations for the first-order nonlinear ODE with IVP.

---

4 Numerical experiments
This section evaluates the performance of the proposed method using sixteen benchmark ODE problems (see Appendix [. These include six first-order, nine second-order, and one fourth- order ODE, covering a wide range of common types such as stiff, linear, nonlinear, and singular ODEs, as well as ODEs with time-varying input signals, undamped free vibration, and higher-order dynamics. Prior to the performance evaluation, we conducted a preliminary analysis involving kernel function selection and sampling strategies across all benchmark problems. Convergence behavior was further examined using Problem 4, a representative nonlinear ODE. Additionally, a comprehensive comparative study will be performed against suitable baseline models to thoroughly assess the efficacy of the proposed approach. All numerical experiments in this study were conducted on a computer system equipped with Intel()) Core TM) i9- 14900HX processor, 32 GB of RAM. Source code is available at: httpp:: thub. com/ Catp/AAAACCComLa///iiiiii SatpappaLLNNNS-SWWWW
The selection of the kernel function plays a critical role in the proposed methodol- ogy, as it directly governs the model's capacity capture nonlinear patterns. To illustrate this, we examine Prob- lem I as representative case. As shown in Figure 2 the radial basis function (RBF) kernel projects samples into an anfinie-dimmnnsionis feature space, offering universal approximation capabilities and robustness to variations in data scale. In comparison to linear and polynomial kernels, the RBF kernel demonstrates superior stability and accuracy when applied to complex nonlinear systems (see results for Problems 4, 5, 6, 14 and 15 Table Consequently, the RBF kernel is adopted in all subse- quent experiments. Both the NLS-SVMs and the classical
RBF kernel justification:
12.5
Andiyical Solsion
Licer Somel P Pdlctin
10.0
Keym Precitionn
RaF-Haral Predition
7.5
5.0
2.5
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Figure 2: Numerical results for Problem I using the Linear, Polynomial and RBF kernel within the NLS-SVMs frame- LS-SVMs tuned via two hyperparameters: the kernel work.
bandwidth o² and regularization coefficient y. A comprehensive quantitative analysis of kernel performance across all benchmark problems provided in Appendix B.3 (Table
Sampling strategy anal- ysis for Nyström land- marks: To assess the in- fluence of sampling strat- egy on the performance of the proposed NLS-SVMs method, we evaluate three distinct approaches that are equidistant sampling, ran- dom sampling, and leverage score sampling. The analysis for Problem 1, summarized in Table I reveals that all three sampling methods achieve a comparable level of prediction accuracy. How- ever, equidistant sampling demonstrates notable advantages in computational efficiency and error reduction. Specifically, it improves computational speed by approximately 16% compared to lever- age score sampling and reduces the RMSE by about 10.66% relative to simple random sampling. Comprehensive numerical results evaluating these sampling strategies across all benchmark problems are available Appendix B33 (see Table These results demonstrate that the equidistant sampling strategy delivers robust and superior performance. It either matches or surpasses the other methods in the majority of test cases, particularly in challenging scenarios such as stiff systems (c.g., Problem 3), second-order ODEs (c.g., Problems 8, 9, 10, and 12), and problems with singularities (e.g.. Problem 13). Furthermore, equidistant sampling exhibits exceptional stability, maintaining consistently low error levels across the entire test cases. While the computational efficiency is generally comparable across all methods, the equidistant strategy frequently attains superior accuracy without incurring substantial computational penalty. In many cases, its runtime is nearly identical to that of the fastest alternative. These collective results indicate that the equidistant sampling strategy provides an optimal balance between numerical precision and computational expense.
Table 1: Performance comparison of NLS-SVMs with different sam- pling strategies Problem
ODEs
Model
MAE
RMSE
1V 0l &c
Time/s
Random
8. 8.70x10
1.06x10
3.06x10
0.52
P1
Leverage score
7.94x10-
9
1
0.64
Equidistant
7.94x10
9.47x10
1.87x10
0.55
Convergence analysis and numerical validation: Building upon the theoretical convergence analysis presented in Section 3.2, which established that the convergence of the Newton-type solver
7

---

is governed by the proper selection of the regularization parameter y and RBF kernel parameter o². this section provides an empirical validation of these theoretical findings. The analysis demonstrated that local quadratic convergence is achieved when these parameters are chosen appropriately. The numerical experiments confirm the theoretical predictions. For the representative nonlinear system (Problem 4), the solver exhibits stable residual convergence, typically within 10 to 50 iterations, as depicted in Figure(). 3). This robust convergence behavior is consistently observed across a wide range of tested problems, with detailed results provided in Appendix (Table The empirical evidence strongly corroborates the theoretical reliability of the solver for the considered class of problems. Furthermore, the numerical results validate the theoretically indicated parameter configurations. Robust and efficient convergence is achieved with 7 = 10(6/T) for all sixteen benchmark problems. The kernel parameter effectively set to x² = 1 1/8/10 for first-order ODEs, and x² = for second-order and higher-order ODEs. See Appendix [B5 for more details, These results collectively demonstrate that the proposed solver delivers numerically stable and efficient convergence across diverse problem types, firmly aligning with the established theoretical framework.
Relative Resioual Norm
2.0
10
Relative Residual Norm
2.0
10
2.0
10
101,
1.5
1-
107, VAM
RASE
10, eaice
1DP, JoMA
Relative Residual
1.5
21
1452
1.0
110
RMSE
RMSE
1.0
0.5
0.5
0.0
20
30
40
50
10
20
30
50
9.0
10
20
30
40
50
Number iterations
Number of iterations
Number of iterations
10
22
30
Figure 3: Numerical results of convergence analysis for Problem 4
Comparative analysis with classical numerical methods: RK4 and EAB: To validate the effectiveness of the pro- posed NLS-SVM method, we conducted a comparative study against two classical numerical techniques: the fourth-order Runge Kutta method (RK4) and the Explicit Adans-Bashforth (EAB) method. The evaluation was performed on three represen- tative test cases: Problem (stiff ODE), Problem 6 (nonlinear ODE), and Problem 11 (singular
Table 2: Performance comparison for solving Problems 3, 6, and 11 using NLS-SVMs, RK4, and EAB
ODES
Model
MAE
RMSE
Train/s
Predict/s
SVS
2
2.77x10
0.460
0.001
P3
RK4
3.87x10
3,41x107
0.880
0.880
EAB
4.97x10
2.90x10
0.030
0.030
NLS-SVMs
8.89x10*
1.01x10
1.230
0.001
P6
RK4
1.59x10'
1.21x10
0.021
0.021
EAB
1,28x10
1.11x10
0.036
0.036
SLSS
2
2
0.100
0.001
P11
RK4
3 6
2.87x10
0.010
0.010
EAB
1.19110
3.
0.020
0.020
On Problem 3, NLS-SVMs achieved an MAE of 2.06 x 10-8, 8, which 18.7 times lower than that of RK4 (3.87 x 10-7) and 2.4 times lower than that of EAB (4.97 x 10-8), while being 880 and 30 times faster in prediction time, respectively. For Problem 6, the proposed method delivered an MAE of 8.89 x 10-4, reflecting a 17.9-fold improvement over RK4 (1.59 x 10-9), along with a 21-fold speedup (0.001 vs. 0.021 s) in predictions. In the case of the singular ODE (Problem 11), NLS-SVMs maintained high accuracy (MAR=2.36 x 10-9), whereas EAB failed with an error nearly 5 million times larger (MAE=1,19 x 10 10-3), establishing the proposed approach as a high-accuracy solver. Additional numerical results are provided in Appendix B.3 (Table 6 It is important to note that the proposed method is not universally superior but offers practical advantages in specific scenarios: (1) it rescues simulations where EAB may fail on singular systems; (2) it enables real-time control applications to stiff systems where RK4 may be computationally prohibitive; (3) it maintains reliability at practical step sizes where other methods require extreme refinement; and (4) it preserves accuracy across range of step sizes, unlike the step-sensitive errors observed with RK4 and EAB.
ODE). The results are summarized in Table 2
Informee Neural Network Baseline Selection: To establish a rigorous baseline for compar- ison with the proposed NI NLS-SVMs method, a systematic model selection procedure was conducted for Physics-Informed Informed Neural Networks (PINNs). The selection process comprised two main stages:
8

---

first, an architectural optimization within the Vanilla PINN framework, followed by a comparative evaluation against Fourier Feature PINN to determine the optimal PINN variant. Vanilla PINN Architecture Selection: An ablation study was performed to identify the most effective network archi- tecture for Vanilla PINN. We compared a standard three-layer configuration against an eight-layer architecture using Problem 1 as benchmark case. The results indicate that the deeper architecture (MAE = 2.32 10 RMSE = 2.93 x 10 -3 Time = 1275 $ provides no significant accuracy im- provement over the three-layer counterpart (MAE 4.45 x 10 RMSE = 5.05 10-3, Time 672 s ), while incurring a substantial computational overhead (1275 s vs. 672 s). Consequently, the three-layer architecture
selected as it offers the optimal balance between accuracy and efficiency for the class of problems under consideration. This model utilizes identity activation functions, a learning rate of 0.001, and the Adam optimizer, with gradients computed via PyTorch's 's automatic dif- ferentiation. Comparative Evaluation with Fourier Feature PINN: The optimized three-layer Vanilla PINN was then evaluated against. Fourier Feature PINN across all 16 benchmark ODEs. Fourier Feature PINN demonstrates superior accuracy in the vast majority of cases. For instance, in Problems 3, 11. and 15, it achieves dramatically lower errors (often by several orders of magnitude) in MAE, RMSE, and Loo norm compared to the Vanilla PINN. Notably, it can also lead to substantially faster training times, as seen in Problems 11 and 15. However, this advantage is not universal. The Vanilla PINN outperforms its counterpart in specific cases like Problems 4, 7, and most notably Problem 16, where the Fourier Feature version fails completely with exceedingly high errors. Furthermore, the Fourier Feature PINN often requires longer computational time, though this is not always the case. Detailed comparison of numerical results are provided in Appendix B.3 (Table 7 Based on this comprehensive analysis, the three-layer Vanilla PINN is selected as the representative PINN baseline for subsequent comparisons with the proposed NLS-SVM method. The chosen model provides a meaningful benchmark for evaluating the relative performance of the proposed method, while Fourier Feature PINN remains recommended for more complex multi-phe multi-physies or iigh-gradient scenarios. All subsequent comparisons with PINNs in this study refer to this optimized three-layer Vanilla PINN configuration.
1.0
Pu Probiier
1.0
21
-25
50
03
05
10
15
28
23
30
10
05
155
20
Problem 1
Problem 2
Problem 3
Problem
-2.
Anglle Saire
322
15
-12
315
0.2
6.0
Iapikil Staton
-11
84
116
1.03
1.18
1.39
120
0.0
0.5
1.
1.5
26
Problem 5
Problem 6
Problem 7
Problem 8
2.5
6:50
6:3
6.
012
37
CA
422
6.2
2.0
60
8.5
Problem 9
Problem 10
Problem 11
Problem 12
3.50
36
25
-1.1
144
28:
4
15
025
03
-33
03
04
1.
38
19
-10
1.
Problem 13
Problem 14
Problem 15
Problem 16
Figure 4: Comparison of model predictions against the analytical solution for Problems 1 to 16

---

5 Results and discussions
Figure presents the comparision of solutions obtained by LS-SVMs, PINNs, and proposed NLS- SVMs for all 16 benchmark ODEs, with corresponding analytical solutions serving as ground truth. By observation, the proposed NLS-SVMs demonstrate significant improvements computational efficiency while maintaining competitive solution accuracy relative to classical LS-SVMs. Table 3 quantifies the significant computational speedups (10-6000 times) under comparable accuracy regimes for key benchmarks. Across all test cases, NLS-SVM, achieved comparable solution accuracy to LS-SVMs (< 0.13% relative MAE, RMSE, and ly - اl difference). Notably, the computational time (including training and prediction time) of NLS-SVM, was reduced by a factor of 10 to 3355 compared to LS-SVMs, depending on the number of time steps n. In comparison to PINNs, the proposed NLS-SVMs model exhibited consistently superior accuracy with maximum reduction by 72% in RMSE, while runs 49-6426 times faster than PINNs (see Appendix B.8 for detailed results).
The scalability of NLS-SVMs is further validated on problems scaled to 100 seconds with 5 x 10 time steps. Under these large-scale conditions, conventional SSVMS and PINNs failed to execute due to memory constraints on the tested hardware, whereas NLS-SVM produced stable and accurate solutions with execution times between 17.5 and 460 seconds (see Appendix [.6 for details). This performance highlights the method's ' practical utility in resource-limited settings. While the framework is validated on ODEs, the underlying Nyström acceleration is readily extendable to large- scale partial differential equations by reformulating the differential operators and to other kernel-based ML tasks requiring efficient approximation. It is important to acknowledge a primary limitation of the current method: although it achieves a significant acceleration for nonlinear ODEs, the computational complexity remains at (n²)) Our future work will focus on breaking this complexity barrier.
Table 3: Comprehensive performance comparison of NLS-SVMs against PINN and LS-SVMs for solving ODEs: accuracy and efficiency
Model
A
AM/%
A ARMSE/%
A lly U /%
Speedup
P4
NLS-SVMs vs PINN
+0.05
L 5.34 x 10"
7.200 x 10"
4 1.81 x 10²
t 19
NLS-SVMs vs LS-SVMs
+0.00
1 4.94 x 10-2
16.14 x 10-2
t 1.29 x 10-1
t 370
NLS-SVMs
à 1.51
1.000 10
1.4.0 10
1
P5
NLS-SVMs LS LVMs
+0.00
1.00
10
44400 10
1.7.0 x 10
10
P16
NLS-SVMs vs PINN
NLS-SVMs vs LS-SVMs
+0.00
+0.00
1.33 x 10-2
1.24 x 10"
1166 10-2
41.35 10
4.333 x 10-2
4 1.779 x 10"
1 58
1 6426
6
Conclusion
This study introduces a novel Nyström--accelertted Least Squares Support Vector Machines (NLS- SVMs) framework for solving ODEs, establishing significant improvement in computational efficiency and scalability over classical ML methods such as standard LS-SVMs and physics-informed neural networks (PINNs). The method successfully overcomes the O(an²) (a 1 for linear ODEs and 27 for nonlinear ODEs) complexity bottleneck inherent in classical LS-SVMs for ODE problems, achieving speedup factors ranging from 10 to 6000 depending on the problem size and characteristics while consistently preserving solution accuracy. Extensive numerical experiments across 16 benchmark problems confirm the model's versatility and robustness, with relative errors remaining below 0.1% across all tested cases. Even for large-scale systems with fine temporal discretization, the solution time is bounded within 460 seconds, demonstrating the practical viability of the framework in real-time engineering applications such as dynamic system control and multi- physics simulation, where conventional kernel methods or neural models are often limited by memory or latency constraints.
Acknowledgements
The authors acknowledge the funding support provided by the National Natural Science Foundation of China (No. 52208485), the Key Research Project of the Educational Commission of Hubei Province of China (Grant No. D20241203), and the Talent Research Startup Fund of China Three Gorges University (No. 2023RCKJ013).
10

---

References
[1] Edward L. Ince. Ordinary differential equations. Courier Corporation, 2012.
[2] Peter Deufthard and Folkmar Bornemann. Scientific computing with ordinary differential equations, volume 42. Springer Science & Business Media, 2012.
[3] George D Byrne and Alan C Hindmarsh. Stiff ode solvers: A review of current and coming attractions. Journal of Computational physics, 70(1):1-62, 1987.
[4] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informe neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
[5] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[6] Markus Heinonen, Cagatay Yildiz, Henrik Mannerström, Jukka Intosalmi, and Harri Lähdesmäki. Learning unknown ode models with gaussian processes. In International confer- ence on machine learning, pages 1959-1968. PMLR, 2018.
[7] JAK Suykens, T Van Gestel, J De Brabanter, B De Moor, and J Vandewalle. Least squares support vector machines. World Scientific Publishing, Singapore, 2002.
[8] Siamak Mehrkanoon, Tillmann Falck, and Johan AK Suykens. Approximate solutions to ordinary differential equations using least squares support vector machines. IEEE transactions on neural networks and learning systems, 23(9):1356-136)) 2012.
[9] Siamak Mehrkanoon, Saeid Mehrkanoon, and Johan AK Suykens. Parameter estimation of delay differential equations: an integration- free Is-svm approach. Communications Nonlinear Science and Numerical Simulation, 19(4):830-8) 841, 2014.
[10] Christopher Williams and Matthias Seeger. Using the nyström method to speed up kernel machines. Advances in neural information processing systems, 13, 2000.
[11] John Dormand and Peter J Prince. family embedded runge-kutta formulae. Journal of computational and applied mathematics, 6(1):19-26, 1980.
[12] Ernst Hairer, Gerhard Wanner, and Syvert P Nørsett. Solving ordinary differential equations I: Nonstiff problems. Nonstiffrroblems. Springer, 1993.
[13] Patrick Conrad, Mark Girolami, Simo Särkkä, Andrew Stuart, and Konstantinos Zygalakis. Statistical analysis of differential equations: introducing probability measures on numerical solutions. Statistics and Computing, 27(4):1065-11) 2017.
[14] Filip Tronarp, Hans Kersting, Simo Särkkä, and Philipp Hennig. Probabilistic solutions to ordinary differential equations as nonlinear bayesian filtering: new perspective. Statistics and Computing, 29(6): 297-1315, 2019.
[15] Philipp Hennig, Michael A Osborne. and Hans P Kersting. Probabilistic Numerics: Computation as Machine Learning. Cambridge University Press, 2022.
[16] Nicholas Krämer, Nathanael Bosch, Jonathan Schmidt, and Philipp Hennig. Probabilistic ode solutions in millions of dimensions. In International Conference on Machine Learning, pages 11634-11649, PMLR, 2022.
[17] Jonathan Schmidt, Nicholas Krämer, and Philipp Hennig. A probabilistic state space model for joint inference from differential equations and data. Advances in neural information processing systems, 34: 13744-12388, 2021.
[18] Andrew J Meade Jr and Alvaro A Fernandez. The numerical solution of linear ordinary differential equations by feedforward neural networks. Mathematical and Computer Modelling, 19(12):1-25, 1994.
11

---

[19] Hyuk Lee and In Seok Kang. Neural algorithm for solving differential equations. Journal of computational physics, 90-11:100-131 1990.
[20] Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal Scientific Computing, 43(5):A3055-A300 2021.
[21] Suyong Kim, Weiqi Ji, Sili Deng, Yingbo Ma, and Christopher Rackauckas. Stiff neural ordinary differential equations. Chaos: An Intendisciplinary Journal of Nonlinear Science, 31)9), 2021.
[22] Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. International conference machine learning, pages 5301-5310 PMLR, 2019.
[23] Amir Gholami, Kurt Keutzer, and George Biros. Anode: unconditionally accurate memory- efficient gradients for neural odes. In Proceedings of the 28th International Joint Conference on Artificial Intelligence, pages 730-736, 2019.
[24] Oliver Hamelijnek, Arno Solin, and Theodoros Damoulas. Physics-infrmed variational state- space gaussian processes. Advances in Neural Information Processing Systems, 37:98505- 98536, 2024.
[25] Frank Dondelinger, Dirk Husmeier, Simon Rogers, and Maurizio Filippone. Ode parameter inference using adaptive gradient matching with gaussian processes. In Artificial intelligence and statistics, pages 216-228. PMLR, 2013.
[26] NR0960 Aluru. A point collocation method based on reproducing kernel approximations. International Journal for Numerical Methods in Engineering, 108611108-11 2000.
[27] Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco. Less is more: Nyström computa- tional regularization. Advances in neural information processing systems, 28, 2015.
[28] Andrea Della Vecchia, Ernesto De Vito, Jaouad Mourtada, and Lorenzo Rosasco. The nyström method for convex loss functions. Journal of Machine Learning Research, 25(360):1-60, 2024.
[29] Petros Drineas and Michael Mahoney. Randnla: randomized numerical linear algebra. Communications of the ACM, 59(6):80-90) 2016.
[30] Cameron Musco and Christopher Musco. Recursive sampling for the nystrom method. Advances in neural information processing systems, 30, 2017.
[31] Philippe Moireau and Dominique Chapelle. Reduced-order unscented kalman filtering with applicationto parameter identification in aage-dimensiona, systems. ESAIM: Control, Optimisa- fion and Calculus of Variations, 17(2):380 405, 2011.
[32] Carl T Kelley, Iterative methods for linear and nonlinear equations. SIAM, 1995.
[33] Ron S Dembo, Stanley C Eisenstat, and Trond Steihaug. Inexact newton methods. SIAM Journal on Numerical analysis, 19(2):400-408 1982.
[34] Stanley C Eisenstat and Homer F Walker. Globally convergent inexact newton methods. SIAM Journal on Optimization, 4(2):393-422, 1994.
[35] loannis K Argyros and Said Hilout. Weaker conditions for the convergence of newton's method. Journal of complexity, 28(3):364-387, 2012.
[36] Susmita Mall and Snehashish Chakraverty. Comparison of artificial neural network archi- tecture in solving ordinary differential equations. Advances in Artificial Neural Systems, 2013(1):18189) 2013.
[37] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 95):987- 1000, 1998.
12

---

[38] Yonghyeon Jeon, Kyung Ryeol Baek, and Sunyoung Bu. Deferred correction neural network techniques for solving ordinary differential equations. Engineering Applications of Artificial Intelligence, 135:108777 2024.
[39] A Prothero and A Robinson. On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations. Mathematics of Computation, 28(12)):145- 162, 1974.
[40] Ivan Ryzhikov, Eugene Semenkin, and Ilia Panfilov. Evolutionary optimization algorithms for differential equation parameters, initial value and order identification. In International Conference on Informatics in Control, Automation and Robotics, volume 2, pages 168-176. SCITEPRESS, 2016.
[41] Nirmal Kumar, Ali Akbar Shaikh, Sanat Kumar Mahato, and Asoke Kumar Bhunia. Applications of new hybrid algorithm based on advanced cuckoo search and adaptive gaussian quantum behaved particle swarm optimization in solving ordinary differential equations. Expert Systems with Applications, 172:1 14646, 2021.
[42] Ali Sadollah, Hadi Eskandar, Joong Hoon Kim, et al. Approximate solving nonlinear ordinary differential equations using least square weight function and metaheuristic algorithms. Engineering Applications of Artificial Intelligence, 40:117-132, 2015.
[43] Luma NM Tawfiq and Ashraf AT Hussein. Design feed forward neural network to solve singular boundary value problems. International Scholarly Research Notices, 2013(1):650467, 2013.
[44] Hongli Sun, Muzhou Hou, Yunlei Yang. Tianle Zhang, Futian Weng, and Feng Han. Solving partial differential equation based on bernstein neural network and extreme learning machine algorithm. Neural Processing Letters, 50(2):1153-1172. 2019.
[45] Shagun Panghal and Manoj Kumar. Optimization free neural network approach for solving ordinary and partial differential equations. Engineering with Computers, 37(4):2989-3002 2021.
[46] William E Boyce, Richard DiPrima, and Douglas B Meade. Elementary differential equations and boundary value problems. John Wiley & Sons, 2021.
[47] Endre Süli and David F Mayers. An introduction numerical analysis. Cambridge university press, 2003.
[48] Alaeddin Malek and R Shekari Beidokhti. Numerical solution for high order differential equations using a hybrid neural wwtwork-optimizatiinn method. Applied Mathematics and Computation, 183(1):260-22) 2700022711120 2006.
[49] Yanfei Lu, Qingfei Yin, Hongyi Li, Hongli Sun, Yunlei Yang, and Muzhou Hou. The Is-svm algorithms for boundary value problems of high-order ordinary differential equations. Advances in Difference Equations, 2019(I):195 2019.
13

---

NeurIPS Paper Checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: The claims in the abstract and introduction regarding the efficiency and accuracy of our proposed NLS-SVMs are empirically validated by the comprehensive experiments in Section 4. Furthermore, the discussion in Section 5 accurately reflects the scope and limitations of our work.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
It is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: See Section 5. It is important to acknowledge a primary limitation of the current method: although it achieves significant acceleration for nonlinear ODEs, the computational complexity remains at (nn)).
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory assumptions and proofs
14

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: Theorem (including proof) and Lemmas are provided in Section 3, and the proofs for Lemmas are provided Appendices A.4 and A.5.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification:
Reproducibility notes are present in Appendices A and B.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
specific model and empirical evaluation, may
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
15

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
The python code implementing the proposed NLS-SVMs method is available at https:///ittun com/ http//AAAcCCooa//Maaaa... SVMs
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co cc/
public/ggiidee/odeSuuuiii httpii//auuaaeeeeodeeeeiiiisiioo wwipbbmmiisiooooiccii.... for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https for more details.
//nips //ips.cc/
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: See Section 4 and Appendix B.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We visualize experimental results with informative figures and tables. These include key metrics like R², RMSE, and MAE.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoo bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
16

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
Question:
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification: See Section 4. All numerical experiments in this study were conducted on a computer system equipped with Intel(R) Core(TM) 19-14900HX processor, 32 GB of
RAM.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics htttp:://eerrpp::
Answer: [Yes]
Justification: The research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- cration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Potential positive impacts include accelerating scientific discovery and engi- neering design, discussed in Sections 5 and 6.
17

---

Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (c.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., denloymemt of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: The answer NA means that the paper poses such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
The source code implemented in this study is proprictary and developed entirely by the authors. No external assets such as third-party code, data, or pre-trained models were used. Therefore, all materials are appropriately credited and utilized in accordance with their original terms.
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
18

---

The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, pwwpsssipppsseeeee thcode com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: Full documentation provided to ensure the complete reproducibility of our work. All newly introduced code is well-documenter and accompanied by usage examples.
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details the dataset/code/moddl as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and sereenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
supplemental material is fine, but if the main contribu-
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
19

---

Justification:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
Declaration LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy (ttppe///orrrpscccoooffrrrrr//// for what should or should not be described.
20

---

A Technical Appendices
A.1 Least Squares Support Vector Machines (LS-SVMs)
The learning objective of LS-SVMs is to find a nonlinear mapping relationship 7 between the input
variable æ; E Rd and the output variable y/ R from the given training dataset ( yo 12=1- The
goal in regression problem is to estimate a model of the form y(ι) w ((x)) + b. A simple description of the mathematical model follows:
minimize 1p (w,ee) w,
1 wwT
w
7 2 _
(22)
subject to = w" (((x))))))) + bb te.,
i I2...
(23)
where 6; E R denotes the error variable, 7 regularization parameter, b E R is the bias term, w E Rh is the weight parameters, <0x;;) E R?th is a high-dimensional feature variable of input data
T; where up.): : Rd + Hth is the feature map and is the dimension of the feature space, and y: indicates the output data. The dual solution is then given by
Ω + I Inxn/1 /
1n
y
(24)
1
0
where Rig = (x 2zy) = (pxx)) (x(xx;, is the ij-th entry of the kernel matrix. 1,, = [l..... 1 e
R", a = [ai..... Onl y [M...... Vn T and Ixx is the identity matrix. The model in the dual
form becomes: y(x) 2im1 (2, + b.
A.2 LS-SVM for p-th order nonlinear ODE
The Lagrangian of the constrained optimization problem becomes 7 E
L (w, b, ei, Eis Vc, 04, Tc B.) 2 ww #2 -
N
Σ ai (wT)(()) (p) (t:) f (t.. vh. (p-1) ...., Uhas Ve. vo.) -ci)
(25)
31 (wI 40 (41) + - v1) -..- 33 s0 (11) - (41) up))
Th (y w" 40 (4:) - b- - Ei)
After obtaining KKT optimality conditions and making use of Mercer's theorem, the solution is obtained in the dual by solving the following nonlinear system of equations:
Ω
hT
0
0 1 (1xxn))(-11)
f
0n-1
S2g
hg
1n-1
In-1
3
U1
h1
hg
12
0 0
(26)
1.1
Bp-1
Op-1
0 __-
0
0
n1
D(y)
1n
0p-1
0z-1
(nn-1)x(n-1)
3
where
2 02
In-1/7,
4n-1/7 D(y) diag (f'
@f(t,v) v)
af v)
f(y) - [S (12,, (12,82), ,...,,, (Ins va .f'(z)
Oy
Oy
k=12.v=w2
a = (02, an) T , = =Im2, Iy,,..., mm] T y = Iy2...., Un T
58
[2 2:n.2:n
2 [1]]
526 [s2]]
ho [1a]]]
[ias] 1.n
2:
2m..2m
1,2 ......
h1 = [ 1,2 i....., [2 i.m
.00-1 = 10, E RR--
21

---

A.3 Nyström Method for Feature Formulation
This appendix section presents the nonlinear feature mapping using Nyström method 7 Firstly, given the input dataset (4+]==1 the corresponding kernel matrix consisting of Gaussian kernel functions is denoted as R(n, n) R=*" where 94 = (44, 1;), i, 1,.... n According to Mercer's theorem, in high-dimensional h(n h S oo) feature space, the following relation can be obtained:
((t;, ( = E (to
(27)
where 2 $2 2 ... 2 0 represent the eigenvalues; Φ1, 22.-.' 0h represent the corresponding eigenfunctions. The eigenvalues and eigenfunctions associated with the following integral equation:
f K [((,, sxdx(z)
(28)
where p(t) represents the continuous-type probability density function of the input space In case of a large-scale training dataset {4, vi w it is assumed that the input space {4:]==] is distributed independently and identically between them, and all of them are sampled from p(t). To approximate the integral equation of the above eigenfunction, m sample data {t., Ms s==1 (m < n)) can be selected from the large-scale training dataset (tz. y/ 1 this case, p(t) can be approximated by the sampled input dataset {t. 1=1 so that the continuous p(t) can be characterized by the discrete {tx [t.]== Then, Eq. can be written as follows:
Σー KZz,t) xo (K0)(2)
(29)
m
For the selected sample data {1.)== the corresponding kernel matrix 9um, m) and the decomposi- tion of its eigenvalues and eigenvectors can be described by the following relation:
Ω mm,m) (T (m)
(m)) (mxm)
(30)
where (m)m)) E Rm.xm represents the small-scale kernel matrix approximating the large-scale orig- inal kernel matrix Dmn, n) by mmm,m) ~ (m,,m) ((m.m)(9(m.mm) [m
- ((m.)) with elements Q. K (tt,, fz), $ (mxx) Rm represents matrix consisting of the eigenvectors of the kernel 8, 2 L...,, m;
matrix 2 (mmm)) where a @x,t), s, k = 1,..., m. and Å (mxm) diag((si;...; &m]) E Rm represents the diagonal matrix consisting of the eigenvalues the kernel matrix St [m, Since S m)) is a symmetric semi-positive definite matrix, each element in A A(mxm) is no less than 0.
Then we can utilize the eigenvalues and eigenvectors of the small-scale kernel matrix Q 1m,, m)) to extract the features based on the proposed NLS-SVMs. By replacing x in Eq (2) with t. and substituting s, z, = 1, ..... q all into Eq. (9) and writing the summary equation in the form of Eq. 30 the following equation is derived:
(x ~ / /m2.as
(31)
Substituting Eq. 31 into Eq. 29 leads to the k-th eigenfunction expression:
(x(tz)
& Σm
(32)
The relation between the high-dimensional vector (tt,t) and the eigenfunction ((t;) can be derived as follows:
Since op (1))((t)))) = (itl))))===VVVeeettl// (k)
TT (())((1)) Em
(33)
substituting this into Eq. 32 and combining with
(1). get the following:
V _m K
px((4) (epp(tt)
(34)
V
The high-dimensional feature vector (p(4) Rm (m < n) can be obtained from Eq. 34 where
(ptt.) = (t((t)), .... (t.)
22

---

A.4 NLS-SVMs for Lemma 1
Proof 1: The constraint conditions include n time-discrete equations where each represents the discretized form of the ODE equation itself at each time point. The Lagrangian of the constrained optimization problem (Eq. (9) becomes:
P Σ 1n
L(w,b, (w,,b,,1,". )1:-.; wnw 7 (w)) (w)
+ (+)(++)++())
( (t:) (w= Σ fl (t))
+ +,g(()))+((t)))))
+ t (w" + 6 v11)
+ + Ap wt (lp--)))(t) Up)
(35)
Then the Karush Kuhn Tuuker (KKT) optimality conditions are as follows:
@L(w,,, b, b, (....... = E A( - f1 (1-11-1)
} w
dw
- [A(Fp 1)]6+ 600 (ti)d +... (1ol0-1) '>00 Ar
b, (....,,,.........
[V(W(()) -f -1 -.. -fo
Ob
(Vp 1)]b+ - Vr
(36)
OL(W, b, AI b,>1,...,;p)) = (p(t)) )" + --
Qh1
aL(w, b, b,,,,,,,,.,,,)))) (vl(=1) Up
Odp
Therefore, these equations in the matrix from gives the linear system in 1q. The model in the primal problem: ((t = wn (p((t+++ +
A.5
Nyström Method for Lemma 2
Proof 2: The Lagrangian of the constrained optimization problem Eq. in its primal form becomes:
2 w w + 7 ( ( (1)) [(1s Mc. (p-1) ,..,, Vi, w1.)) T
L(w, b, Ap. vi)
2
(F) f(le, Ve (p-1) .... Vb, Vo.))
(y:. - w (((.))
2
(v,, w ((((tt)) - b) + Xp (wT (0 (t2) -Up
+ +... A1 (w p(t)) + + v1)
(37)
Then the KKT optimality conditions are follows:
23

---

@L(w, b.A.......y yi) Y
Jacobian
0Y,
Ow
aww b, A1 ,.....,,,pp
,p: yi)
@L(w, b,A1,...,,, yi)
Jacobian
0Y2
0b
Iw. b, A1 A1...... Ap: Yi)
@T(Ww, b, b,,,,,,,......................,,
:..11...;Ap.. Ap.
Jacobian
0Y3
Y
011
Oww, b, A1....... Ap: y1)
(38)
aL(w. b.A1...... Bi)
Jacobian
0Y p+ O(w, b, A...... Ap. yi)
yp+2
AAp
OL(w, A......., Ap Ui)
Jacobian
p+
p+3
Oy:
Mw b, X1,
Ap,
Therefore, these equations the matrix from gives the nonlinear system Eq. primal form becomes ÿ(t) w (t()) + b.
The model
A.6 NLS-SVMs for Solving First-Order ODE with IVP
As first example, consider the following first-order IVP:
(1) (t) r((), = 4 4st5 t..
(39)
Then start by assuming the approximate solution to of the form j(t) = wn w In the NL.S- SVMs framework, the approximate solution can be obtained by solving the following optimization problem:
_ w 2 5 D
minimize whb.e
v1(1,) Wn ' (t.) (1) +b (r,) + i= 2.2.....
s.t.
u(t1) = wf +b=
(40)
The Lagrangian of the constrained optimization problem becomes
L(w, w 2 (" fi(t)) ( (tt)) +b - r(ti))
2
(w" (1: [w' ppt:)) + -r(tt)
+ + ((1)) + -b- v1)
(41)
Then the Karush- Kuhn -ucker (KKT) optimality conditions are as follows: (L(,,b
Ow
O, @L(w.b.)) = 0, (L(w))))) 0, therefore, the solution to Eq. 41 is obtained by solving the
2h
82
following primary problem:
E + A (w' ff
Af (ff.) 1) (f.1)
(t))
w
Ar Vr
V v.
b
(42)
4 (ti)
0
01
where E R9x9 is identity matrix; V [x((1), m.. p((n))
)T E (fn-1)xm, T
(t2)) (T,))) R(n-l)xm,
[f fittt) f(,t)) R=-l, 1 [1...., 1|
Rn-1, 10 (4) plxm, 7/ fU-l" E R (r(12), ((()))".)((())))"
f1j E Rix(n-1),p
E ER=L: E mxx1, 7 is a regularization parameter.
The model in the primary form becomes ((t) T7 (t() +

---

A.7
NLS-SVMs for Solving Second -rdde ODEs with IVP and BVP
IVP Case
Let us consider a sedond-order IVP of the form:
(t) = (t) + (t())x(() + r(t), v(t1)
y((1)) = 02 t 5
(43)
+ bb then obtained by solving the following optimization The approximate solution, ((t) = ((t) problem:
1
w" 7 -
minimize
whb,o
2
2
(t,,,) w= (t,) = f((tt)/"'''' (t)) + f2(4) [ 4 (t+) + +6] +r (t,) e,
i == .....,
5((() w= ((1)) 122
y(11) - (t1) +b= 01
(44)
The Lagrangian of the constrained optimization problem becomes L(w, (1.
1 www + 2 (0 f((tt)="''(ttt) (
fa(ti) ((t,)) +b -((.)) T
w (t))"(t - fi (ti (' - f2(t.) (w" ((t)) +b - rr(tt)
+ +1 (w ((1)) + b v1) + A2 (w (t())) (ti) -22)
(45)
Then the Karush KoaisshKKhnTTeee Tucker (KKT) optimality conditions are as follows: 01...b.d1.do) dw
d. 0, (L(w.bAA,,)) 0, therefore, the solution Eq. 0, (L(..b,A1,2)) 0,
45.
is
8A1
8A2
obtained by solving the following primary problem:
E A Ag" 4/
P2
-A (f2 1)
6"' (t1)
(t)(t)))
E
Ar
f1 V
f22
-V (ff-1)) 1)
b
Vr
(t1)
0
0
A1
(t1)
d2
U2
(46)
where E E Raxm is an identity matrix:
= (p(12), ..., ((m)))
(n=-ijxm;
W (2'(12)..... (',))) E R(n-tixm.qT"
(12( (t (Ln)) T (n-1))mm
fi
[f1(12),... [i(t..)) E RR--1; fz = (f2(t2) .... )1 E 11-1;1
11 [1..... 11 E Rn-1,
((1) E R2xq; (41) E A 7/1"" S. S2 R= = 1 121 RIx(a-1);r =
E ER"-1; mx11 y is regularization parameter.
The model in the primary form becomes (() w (t)) b. BVP Case Consider the second-order BVP of ODEs of the form
(47) Then the parameters the closed-form approximation of the solution can obtained by solving the following optimization problem:
(t( fi(t''((t) + fa(t)y(t) + r(t), y(t1) a1 y(t)) on ti t
5 ww 7 En
minimize w.b.e
s.t.
v'(t.) T _T (t.) = ((()) (wT
(t.) + +r((tt)) e, 2, ....,
v(t)) -- ( +b=
y(t,)) = (t.) 4n
(48)
The same procedure be applied to derive the Lagrangian and afterward the KKT optimality conditions. Then, one can show that the solution to Eq.48) linear system:
(88) is obtained by solving the following
E + A A V' vf -A (fz 1) -V
t1
(tw)
w
Ar
V "T"
fi
f2
(fz 1)
b
Vr
4 (4)
0
0
A1
Q1
(1,))
Ag
In
(49)
25

---

where E E Rmxm is an identity matrix; yf = (x(t2), 0 11 )1" R(n-2)xm, W
( (tn-1))
R(n=2)xm..p"
(t((
T E p(n-2)xm, f1 =
(f((2),.., (1 E R==- f2 (f2(1)) , ..... f2 mn-
n 1 [...... 11" E R=-2;
30 (8)) E 1ammm;pp ;p (tx) E A ? _ f. f
y [m 2V 1 f21
g1x(a-2))r= (rtt))...,rtt..))) r(t.))) e I 2 R. mmx1, y is regularization parameter.
The model in the primary form becomes ((t) = w (t)) + b.
A.8 NLS-SVMs for Solving First-Order nonlinear ODE
Here, we take Eq. as the object of study. The Lagrangian of the constrained optimization problem becomes
1 ww + [WWWTp (')))((1)) -1(1, vi) (')))(())
L(w, b, X, y:)
- [(11,
2
2 (y: - p(t.) b) (y: w² (t,)) b) + x (w- (p(tt) + b- vz). (50)
2
Then the Karush KuhhhhKKun---hhhhh -ucker (KKT) optimality conditions are as follows: aL(W.b,). Y1, BL(w.b.A.... Y2
aw
b.w.w)) Ya,
aL(w,b.Au. Y4, therefore, the solution to 50.
8b obtained by solving the following primary problem:
BA
By
8Y1
8Y
8Y
8Y1
0Y:
0Y;
0b
0Y2
012
w
b
Y
Y2
8X
Dy:
(51)
0Y;
0Ys
0Y3
0Ys
À
Y3
Bw
8Y4
8b
0Y4
8A
8Y4
y;
Y
8
8b
8A
8y
where Y1 E Rmx1 Y2 E RIxL,Y3 E R
RIxI,Y4 E R(a-1)x1, The nonlinear system, which
consists of equations with unknowns (w, b, X, y:). solved by Newton's method. The model in the primary form becomes ((t = w" w='(((t)== + +b.
B
Supplementary Material
B.1 Sixteen benchmark ODE problems
Problem 1: Consider the first-order ODE with time varying coefficient g
36
37
(1(
+
(t+
1+3(" y= 1 2 + 142+13
143r
e 3], (0) =
(1+--13
Problem Consider thd Dahlquist's problem non-stiff linear problem 38 given by y'())
-y(1), E (0, v())
Problem 3: Consider well-known stiff test problem, the Prothero-Robinson equation
38
g((t) = -y0²(t((t) - sin 1) + cos t, E (0, 10]: y(0) =
Problem 4: Consider the following nonlinear first-order ODE, which is studied by Junaid et al
40 (t) yt) >((r)+ E (0, 3), (0) =1
Problem 5: Consider the following nonlinear first-order ODE [38 defined by st()) = (t),
(0, 1], =() -1
Problem 6: Consider the following nonlinear first-order ODE
defined by y'(t)
2 v()) In7 + 1-,,/ E (1,1 =0
Problem 7: Consider the following second-order IVP 42 defined by y"(() + 01)) v() =
1, (0, 20], =0) = 0, y (0) = O, the interval solution for this test problem is considered from to 20, which is assumed to be wide range for solutions.
Problem 8: Consider the second ODE having the initial conditions based on the laws of physics as follows 42 (y"(() + 18y(t) = 0, t E (0, 2). y(0) = - - y (0) =
Problem Consider the following second order damped free vibration equation 36
4y'(t) + 4t()) = 0, 10). v(0) (0)
26

---

Problem 10: Consider the following second-order ODE with time-varying input signal 37 :f)
(t((t) + y(t) = - (4/5) cos t, & E [0, 10), v(0) = 0, (0) =1
Problem 11: Consider a second-order singular ordinary differential equation 43 : +
cos(t) + sin(t) 0, E [o, 11, (0) 1, (1)) = cos(1)
t
Problem 12: Consider the following second-order ODE
3/ + y(t) = 2, [0. 11, y(0)
1. v(1) =
Problem 13: Consider singular differential equation 45 43 ky +1((+++++)) //)=)//////++++++++)" (()) 0,
[0, 1]. v(0) 1,y(1) =
Problem 14: Consider one of the most famous second order ODE. It is given as follows 42
/y (t) (v'( - 2e (t) = 0, t E (0, 1), v()) = 0, y(1) =
Problem 15: Consider the following nonlinear second-order ODE
defined by y"(t)
-y(t) + ( v(t)
t E 11, 1- = y(1) = 0.324027137
Problem 16: Consider fourth-order linear ODE
48
49
defined by (1)) 120t,
1-1,1), 1, = 5, v1) = 3, (1)) =
B.2 NLS-SVMs for Solving Second-order ODE in Terms of Problem 15
Problem15: 5²(1) = -v( + 2(v(())² v(1) , E (111,11,,((1)) = y(1) = 0.324027137
1 minimize-w - w+ + - ン w.b.er...11.
(52)
s.t. v²(+) = -- (t""( (ti) = vv(4) + 22 (t:)t (t:)))/u(tt) te,, i= 2= n2 -1 ((44 - 4 (t)). i= 22 -
v(4:) ((4) +b+ Eo, i= .
(11) +b= !1 (t.) + 4n-
The Lagrangian of the constrained optimization problem becomes
w" w + 2 w - (1gg)))))))()
2
w T (1,) (-v(t+) 2(0((4))) + (y - ((.)) - b)
(u: wF - b) + x x(1) + n1) Xo (w= (((x) + b - qn)
(53)
Further simplified:
2 w Tw + 2 2w T (t,) - (- + +2( ('t))) T
L(w, b, A1, X2, H.)
2
w 6" (1,) - (- + 2 Φ
+ 3 (y: - s(t,) - b)T
( NF - b) A (w) ((t)) + q1) Xe (w² pp(x)) + b qn)
(54)
Then the Karush-Kai -hucerr (KKT) optimality conditions follows:
27

---

@L(w, b, λ1, 22,
Y1.
( (5% - 2(4/), w)² y + y:ll + + 3 - k.. b1]
+ 7 [(((")))) [40 (( /y: ysl 4(()1, " /
dw
+ A1s + X24pt.
(55)
Y2: @.(Ww, b, \1, X2 y1) 71 - [1 m. bTl x1 + X2
(56)
Øb
aL(w, A,,1, X2, v:)
Y3:
11:, + -b-00 41
(57)
0d1
Y4: OL(W, b, Å1 A2, Vi)
pix + qn
(58)
0A2
@L(w. b, X1, A2, yi) = 7 (2((0' w)2 /2/// + 2 w00" w) 33
Y5:
8y:
+y/: 4i, w bf
(59)
Construct the Jacobian matrix and establish the matrix equation:
8Y1
8Y
8Y1
8Y
8Yi
8Y2
8w
8Y2
8Y2
8A1
01g
8Y2
872
w
Y
8A1
0A2
Py.
b
Y2
0Y3
0Y3
8Y3
0Y3
Y3
(60)
Ow
0X
0A
0Y1
0Y1
0Y1
0Y
DY
A2
8Y,
8Y,
8Y
aY,
YBS
8w
6b
0A1
0Ag
Fill the results of the Jacobian matrix calculation into the equation matrix:
OY, dw > II
ITI
(4P"))
9Y
Pr..
Y
III
1
HY2
b
A1
(61)
0
A2
8Yss
I
Yi
Ys
dw
Ty:
where
0Y1
+ +>(""")" [49", -4diag 3//,w
[p/,w
5 w
1/'
47(4)%))²
diag
49 4 diag
4pl.
dw
yi
yi
(p" 2(1p1,w)2
+ diag
+yi] diag
s0',
+ x((p)) (T T
(62)
V
y!
0Y1
12(xk),w) )2
= 76p" diag
+1
ay.
V
2(10), w)2
-474' diag
T
u.
y
[2 w)
1
-44t T
Yi
y
(63)
28

---

0Y5
2(4p/, w)
- diag
50%
diag
'4',
Ow
2(1),,w)
+ 7
diag
kp,Ww
+ yi
ddia
7,p,
(64)
yi
Yi
0Y5
2((/, w)2
2(4)),w)
2((4), w)2
(2(9/,w) w)2
y diag
Vi
+1
+1
+1
Oy.
yh
y
yi
y¹
(65)
Use Newton's 's iteration to solve for w and b, and then construct the prediction equation in the original space: (tt) w (p())))))))))) b.
B.3 RBF kernel justification, Sampling for Nyström landmarks and PINN baseline comparison
Table 4:
Comparative performance evaluation of Linear, Polynomial, and RBF Kernels for NLS-
SVMs in Solving ODEs
ODEs
Model
MAE
RMSE
1V - Úll
Linear
1.70x10
2.00x10'
5.05x10
Problem
Poly
1.92x10
2.45x10²
1.29x10'
RBF
7.94x10*
9.47x10
1.87x10
Linear 1.54x10"
2.15x10'
1.13110'
Problem 2
Poly 3.37x10
3.91x10'
7.44x10'
RBF
3.16x10
3.77x10
9.75x10
Linear 6.34x10*
7.15x10"
1.24x10"
Problem
Poly 1.19x10-1
1.43x10
4.80x10
RBF
8.06x10
2.77x10
1.40x10
Linear 3.29x10
4.67x10'
1,24x10'
Problem
Poly 3.16x10
4.45x10
1.16x10'
RBF
5.47x10
6.92x10
1.5011
Linear 1.0xx10'
1.44x10'
5.82x10"
Problem 5
Poly
1.53x10
1.76110
2.84x10
RBF
8.79x10
1.07x10
1.67x10-
Linear 4.01x10'
5,05x10"
1.6x100'
Problem
Poly
9.67x10
1.10X10
1.85x10
RBF
8.89x10
1
3.70x10
Linear 4.68x10-1
6
1.56110"
Problem 7
Poly
3
1.47x10
RBF
1.43x10
3.02x10
Linear
4.30x10
5.23x10
2.01x10'
Problem 8
Poly 3
3.93x10'
6.70x10
RBF
5.66x10
7.05x10
2.23x10
Linear 1.91110
2.60x10"
1.59x10'
Problem
Poly
8.05x10
9.67x10
1.73x10
RBF
5.74x10
6
2.21x10
Linear 1.17110'
1.6x10"
8.50x10
Problem 10
Poly
5.23x10-1
6.25x10-'
1.3xx10
RBF
4.29x10-
6.06x10
3.09x10
Linear 9.62x100
1.27x10'
5.79x10
Problem 11
Poly 7
3.00x10
RBF
2.36x10
2.97x10
9.43x10
Continued on next page
29

---

ODEs
Model
MAE
RMSE
11 5
Linear
9.75x10'
1.27x10
5.54x10
Problem 12
Poly RBF
2.75x10
3,24x10
5.02x10
3.15x10
3.94x10"
1.3x10
Linear
4.55x10"
6,58x10"
2.84x10'
Problem 13
Poly RBF
1.27x10'*
1.59x10*
2.91x10
3.48x10**
4.36x10
1.45x10
Linear 2.12x10'
2,63x10'
7.68x10'
Problem 14
Poly
2.38x10
2.82x10
4.42x10
RBF
1.17x10
1,57x10
3.85x10
Linear 1.08x10*
1.22x10'
1.76x10
Problem 15
Poly RBF
2.93x10'
3.30x10'
7.96x10'
1.68x107
1.98x10"
3,46x10*
Linear
4.27x10'
6.18x10
3.53x10
Problem 16
Poly 1.67x10"1
1.92x10'
2.86x10'
RBF
1.72x10
2.06x10
5.07x10
Table Comparative performance evaluation of NLS-SVMs utilizing Equidistant, Random, and Leverage-based sampling strategies for solving ODEs
ODEs
Model
MAE
RMSE
1ly - Āl
Time/s
Random
8.70x10
1,06110*
3.06x10*
0.52
Problem
Leverage score
7.94x10*
9.47x10
1.87x10
0.64
Equidistant
7.94x10*
9.47x10
1.87x10
0.55
Random
3.41x10
4.06x10
8.14x10"
0.52
Problem 2
Leverage score
3.22x10
3.82x10
7.5
0.57
Equidistant
3.16x10
3.77x10
9.75x10
0.55
Random
1.31x10'
1.75x10"
8.99x10'
0.53
Problem
Leverage score
3.12x10
4.02x10
1.79110
0.53
Equidistant
2.06x10*
2.77x10
1.0x10
0,56
Random
5.47x10
6.92x10
1.50110
9.53
Problem 4
Leverage score
5.47x10
6.92x10
1.50x10
Equidistant
5.47x10
6.92x10-1
1.50110
10.2
Random
8.79x10
1.07x10
1.67x10
1.07
Problem 5
Leverage score
8.79x10
1
1
1.02
Equidistant
8.79x10
1.0x10
1.67110
0.97
Random
8.89x10*
LO1x10
L.7x10
0.80
Problem
Leverage
L01x10
1.70110
1.45
Equidistant
8.89x10
L01x10'
L
0.69
Random
1.43110
166x10
3.02x10*
0.45
Problem 7
Leverage score
1
1.66100'
3.02x10
0.54
Equidistant
1.43x10
1.67x10'
3.02x10
0.50
Random
9,92x10
4.87x10
0.44
Problem 8
Leverage score
6.90x10
8.79x10
3.02x10
0.46
Equidistant
5.66x10
7.05x10
2.23x10
Random
8
111x10
3
0.43
Problem
Leverage
6.07x10
7
2
0.45
Equidistant
5.74x10
6.92x10
2.21x10
0.47
Random
1.63x10'
2.25x10'
L
0.44
Problem 10
Leverage score
2.29x10
2.65x107
8.12x107
0.43
Equidistant 4
6.06x10
3.09x10
0.49
Random
9
1
2
0.43
Problem 11
Leverage score
2.46x10
3.06x10
8.66x10
0.42
Continued on next page
30

---

ODEs
Model
MAE
RMSE
1 ا
Time/s
Equidistant
2.36x10*
2.97x10"
9.43x10
Random
2.87x10*
3.59x10
1.13110*
0.42
Problem
Leverage
1.35x10
L68x10"
6.31x10²
0.42
Equidistant
3.15x10
3.94x10
1.3x10
0,46
Random
8.76x10"
1.1x10
4.29x10
0.42
Problem 13
Leverage score
5.48x107
6.89x107
2.28x10
Equidistant
3.48x10²
4.36x10
1.45100*
0.48
Random
2.07x10
2.67x10*
7.13x10
1.38
Problem 14
Leverage score
1.95x10
2.63x10
6.66x10
1.39
Equidistant
1.17x10
1.57x10
3.85x10
1.39
Random
1999100
2.35x10"
4.40x10*
3.54
Problem 15
Leverage
1.73x10?
2.04x10"
3.60x10²
3.19
Equidistant
1.68x10?
1.98x10
3.46x107
3.25
Random
3.57x10
4.42x10*
1.5410
0.43
Problem 16
Leverage score
2.09x10*
2.57x10
8.33x10
0.45
Equidistant
1.72x10
2.06x10
5.07x10
0.42
Table 6: Comparative performance evaluation of NLS-SVMs, RK4, and EAB methods for solving ODEs
ODEs
Model
MAE
RMSE
Train/s
Predict/s
NLS-SVMs
7.94x10
9.47x10
0.550
0.001
Problem
RK4
2.20x10
1.80x10*
0.011
0.011
EAB
3.23x10
1.74x10
0.013
0.013
NLS-SVMs
3.16x10"
3.77x10
0.480
0.002
Problem
RK4
1.42x10
8.96x10
0.019
0.019
EAB
2.93x10'
1.36x10
0.021
0.021
NLS-SVMs
2.06x10
2.77x10
0.460
0.001
Problem
RK4
3.87x10
3.41x10
0.880
0.880
EAB
4.97x10
2.90x10
0.030
0.030
NLS-SVMs
5,47x10
6.92x10*
14.530
0.001
Problem
RK4
3.79x10
2.78x10
0.022
0.022
EAB
9.41x10
6.45x10
0.016
0.016
NLS-SVMs
8.79x10
1.0xX10
4.590
0.001
Problem
RK4
3.58x10
3.38x10
0.013
0.013
EAB
1.79x10
1.41x10
0.011
0.011
NLS-SVMs
8,89x10
1.01x10
1.230
0.001
Problem
RK4
7.94x10
2.61x10
0.021
0.021
EAB
1.28x10
1.11x10
0.036
0.036
NLS-SVMs
1.43110*
1.67110*
0.150
0.001
Problem
RK4
1.66x10-
1
0.012
0.012
EAB
4.01x10
3.49x10
0.013
0.013
NLS-SVMs
5.66x10*
7.05x10*
0.066
0.001
Problem
RK4
3.61x10
2.81x10
0.023
0.023
EAB
1.26x10
7
0.018
0.018
NLS-SVMs
5.74x10
6.92x10
0.170
0.001
Problem
RK4
8.21x10
3.08x10
0.022
0.022
EAB
2.32x10
1.11X10
0.026
0.026
NLS-SVMs
4.29x10
6.06x10*
0.075
0.001
Problem 10
RK4
1.00x10-
1
0.029
0.029
EAB
5.48x10
0.018
0.018
Continued on next page
31

---

ODEs
Model
MAE
RMSE
Train/s
Predict/s
NLS-SVMs
2.36x10*
2.97x10"
0.100
0.001
Problem 11
RK4
3.11x10
2,87x10*
0.010
0.010
EAB
1.19110
6.86x10
0.020
0.020
NLS-SVMs
3.15x10"
3,94x10"
0.043
0.001
Problem
RK4
5.75x10
5.20x107
0.016
0.016
EAB
5.38x10
4.61x10
0.023
0.023
NLS-SVMs
3.48x10
4.36x10*
0.098
0.001
Problem 13
RK4
1.57x10
1.21X10
0.011
0.011
EAB
4.79x10
4.15x10
0.024
0.024
NLS-SVMs
1.17x10
1.57x10
1.390
0.001
Problem 14
RK4
9.42x10'
8.18x10
0.025
0.025
EAB
9.75x10
7.57x10
0.019
0.019
NLS-SVMs
1.68x10'
1.98x10'
1.690
0.001
Problem 15
RK4
1,0x110
9.10x10
0.023
0.023
EAB
1.35x10
9.89x10
0.020
0.020
NLS-SVMs
1.72x10
2,06x10*
0.410
0.001
Problem
RK4
3.75x10
4
0.019
0.019
EAB
3.57x10''
4.36x10
0.023
0.023
Table 7: Comparative performance evaluation Vanilla PINN versus Fourier Feature-enhanced PINN for solving ODEs
ODEs
Model
MAE
RMSE
|1v - úl oc
Time/s
Vanilla PINN
4.45x10
5.05x10
9.91x10
672.00
Problem 1
Fourier Feature PINN
1.02x10
1.55x10
3.41x10
549.00
Vanilla PINN
2.03x10
2.09x10
3.26x100
403.00
Problem 2
Fourier Feature PINN
1.33x10
1.64x10
6.43x100
762.00
Vanilla PINN
1.82110"
2.88x10
1.03x10"
1209.00
Problem 3
Fourier Feature PINN
8.51x10
1.15x10
3.72x10
1556.00
Vanilla PINN
5.35x10'
7.21x10
L81x10"
496.00
Problem 4
Fourier Feature PINN
5.57x10
6.67x10
1.50x10'
956.00
Vanilla PINN
L.03110²
L00x10
L.13110
166.00
Problem 5
Fourier Feature PINN
1.90x10
2.26x10
4.41x10
166.00
Vanilla PINN
1.4x110
1.76110
3.57x10
75.60
Problem 6
Fourier Feature PINN
2.32x10
2.47x102
3.86x10'
217.00
Vanilla PINN
6.58x10
8.06x10
1.63x10'
1021.00
Problem 7
Fourier Feature PINN
7.88x10'
8.35x10'
1,42x10
1279.00
Vanilla PINN
5.23x10
5,92x10
L.0x100
1392.00
Problem 8
Fourier Feature PINN
L.13110
1.26x10
L87x10
2170.00
Vanilla PINN
8.66x10-
9.96x10-
776.00
Problem 9
Fourier Feature
8.35x10'
8.72x10'
9.55x10'
1534.00
Vanilla PINN
3.12x10
4.07x10*
8.62x10'
636.00
Problem 10
Fourier Feature PINN
2.61x10
3.31x10
7.02x10'1
2364.00
Vanilla PINN
L66x10
L89x10
3.32x10
817.00
Problem 11
Fourier Feature PINN
2.81x10
3.48x10
1.16x10
16.30
Vanilla PINN
2.28x10"
2.29x10
2.60x1003
969.00
Problem 12
Fourier Feature PINN
1.56x10
1.93x10*
5.91x10
16.70
Vanilla PINN
4.22x10
5.72x10
1475.00
Problem 13
Fourier Feature PINN
2.77x100
3.42x10
9,56x100
33.60
Vanilla PINN
8.22x10
8.77x10
112.40
Problem 14
Continued on next page
32

---

ODEs
Model
MAE
RMSE
11 Ø
Time/s
Fourier Feature PINN
7.38x10"
9.17x10²
2.02x10
36.00
Vanilla PINN
2.42x10'
2.62x10*
3.62x10*
230.00
Problem 15
Fourier Feature PINN
2.14x10
2.66x10*
6.68x10
85.40
Vanilla PINN
1.24x10
1.35110
180x10
2699.00
Problem 16
Fourier Feature PLNN
7.51x10'
8.74x10'
2.38x10'
697.00
B.4
Convergence analysis of iterative solver
Table 8: Convergence behavior and computational efficiency of across varying iteration counts for solving ODEs
ODEs
Index
50
100
200
300
400
500
Residuals
5.1x10**
3.9x10
3.5x10
1.6x10
6.4x10
1.1x10
Problem 4
RMSE
6.92x10*
6.92x10
6.92x10
6.92x10*
6,92x10
6.92x10
Time/s
10.2
24.1
53.0
81.2
107.4
135.0
Residuals
5.1x10'
5.1x10'
5.1x10'
5.1x10'
5.1x10
5.1x10'
Problem 5
RMSE
1.07x10
1.07x10
1.07x10
100x10
1,07x10
1.07x10
Time/s
1.8
12.0
30.5
48.8
65.5
84.0
Residuals
6.0x10
6.0x10
6.0x10
6.0x10
6.0x10
6.0x10
Problem 6
RMSE
1.0x10 0.7
100x100>
1.01X10 1.7
1,0
1.0x10
1.0x10 7.7
Time/s
0.9
2.9
5.4
Residuals
5.1x10
5.6x10
6,4x10
7.3x10'
8.3x10"
9.5x10
Problem 14
RMSE
4.97x10 0.6
5.36x10 0.7
1.84x10
1.58x10
1.57x10 1.3
1.57x10 1.4
Time/s
0.8
1.0
Residuals
1.9110
1.410
8.6x10
2.9x10
2.5x10
7.1x10"
Problem 15
RMSE
1.98x10
1.98x10
1.98x10
1.98x10 2.2
1,98x10" 3.3
1.98x10
Time/s
0.8
0.9
1.3
6.6
B.5
Hyperparameters for Training
Table 9: Hyperparameter configurations for PINN, LS-SVMs and NLS-SVMs in solving ODEs
ODEs
Model
Total sample
Sub-sample
Fpoch/Iteration
Y
PINN
1000
10000
Problem 1
LS-SVMs
1000
10
107
NLS--VMS
100
150
10
107
PINN
10000
10000
Problem 2
LS-SVMs
10000
10
107
NLS-SVMs
10000
50
10
10
PINN
10000
20000
Problem 3
LS-SVMs
10000
10
107
NLSS-SVMs
10000
10
107
PINN
1000
10000
Problem 4
LS-SVMs
1000
50
10
10c
NLS-SVMs
1000
50
50
10
10c
PINN
600
10000
Problem 5
SS-S
600
25
8
105
NLS-SVMs
25
8
10
PINN
10000
Problem 6
LS-SVMs
300
50
1
10
Continued on next page
33

---

ODEs
Model
Total sample
Sub-sample
Epoch
a²
7
NLSS-SVss
50
10
PINN
1000
10000
Problem 7
LS-SVMs
1000
1
10"
NLS-SVMs
1000
10"
PINN
1000
20000
Problem 8
LS-SVMs
1000
1 10
NLS-SVMs
1000
1
10
PINN
1000
10000
Problem 9
LS-SVMs
1000
1000
1 107
NLS-SVMs
1000
50
1
107
PINN
1000
10000
Problem 10
LS SV-S
1000
1
107
NLS-SVMs
1000
50
1
107
PINN
1000
10000
Problem 11
LS-SVMs
1000
1 10
NLS-SVMs
1000
20
1 107
PINN
1000
10000
Problem 12
LS-SVMs
1000
1 107
NLS-SVMs
1000
20
1
107
PINN
1000
10000
Problem 13
LS-SVMs
1000
1
NLS-SVMs
1000
20
1
107
PINN
200
10000
Problem 14
LS-SVMs
200
500
1
107
NLS-SVMs
200
10
500
1
107
PINN
300
20000
Problem 15
LS-SVMs
300
400
107
NLS-SVMs
300
10
400
1
107
PINN
1000
10000
Problem 16
LS-SVMs
1000
107
NLS-SVM SVMS
1000
20
1 107
Table 10: Hyperparameter configuration for NLS-SVMs in large-scale ODEs solving
ODEs
Model
Total sample
Sub-sample
a²
?
Problem
NLS-SVMs
[0.100]
50000
10"
Problem
NLS-SVMs
[0,100]
50000
200
10"
Problem
NLS-SVMs
[0.100] J
10000
100
8
105
Problem
NLS-SVMs
[0,100]
50000
100
10"
Problem
NLSS-VMs
[0,50]
20000
10"
Problem
NLS-SVMs
[0,100
50000
200
10"
Problem 10
NLS-SVMs
[0,100]
50000
200
10"
Problem NLS-SVMs
10,100
50000
200
10"
Problem 12
NLS-SVMs
10,100
50000
200
10
34

---

B.6 Assessment of NLS-SVMs for long-time integration of ODEs
10000
Analytiod faiution
Predition
6000
-0.4
0.0
4000
-100
-1.0
20
40
60
80
100
20
40
100
100
20
108
Problem
Problem 3
Problem
1.0
0.4
-0.4
0.8
100
10
20
30
40
50
0
20
40
60
60
100
Problem 7
Problem 8
Problem
0.5
0.00
-0.5
-025
80
100
00
100
100
Problem 10
Problem 11
Problem 12
Figure NLS-SVMs predictions against analytical solutions for twelve benchmark problems
B.7
Performance comparison of PINN, LS-SVMs and NLS-SVMs
Table 11:
Comparative performance evaluation of PINN, LS-SVMs, and NLS-SVMs in solving ODEs: accuracy metrics and computational efficiency
ODEs
Model
R²
MAE
RMSE
11y - ~ll
Time/s
PINN
0.99999
4
5.05x10*
9.91x10
672.00
Problem 1
LS-SVMs
0.99999
4.81x10
5.72x10
1.17110
9.79
NLS-SVMs
0.99999
7.94x10
9.47x10
0.55
PINN
0.99989
2,03x10"
2.09x10"
3.26x10*
403.00
Problem 2
LS-SVMs
0.99999
6,87x10
8.55x10
3.07x10
1587.00
NLS-SVMs
0.99999
3.77x10
9.75x10
0.55
PINN
0.74911
1_82100
2.88x10'
1.0x100'
1209.00
Problem 3
LS-SVMs
0,99999
2
3.23x10
1.36110
1879.00
NLS-SVMs
0,99999
2.06x10
2.77x10
1.4010
0.56
PINN
0.94554
5.35x10
7.21x10'
1.8x10
496.00
Problem
LS-SVMs
0,99999
5.26x10
7
2.13x10
3774.00
NLS-SVMs
0,99999
5
6.92x10
10,20
PINN
0.99995
1.03x10
1.04x10*
1.13100
166.00
Prohlem 5
LS-SVMs
0,99994
8
1.1x11
9.71
NLS SVMs
0.99994
8.79x10
1.07x10
1
0.97
PINN
0.89021
T
L.76x10
3
75.60
Problem 6
SSSS
0.99993
4.47x10
5.11x10
8.63x10
13.20
NLS--VMM
0.99974
8.89x10
1,01X10
1.70x10
Continued page
35

---

ODEs
Model
R²
MAE
RMSE
11 - ~l
Time/s
PINN
0.91449
6,58x10
8.06x10
1.63x10-t
1021.00
Problem 7
LS MS-S
0.99997
1,43x10
1.66x10*
3.02x10*
22.90
NLS-SVMs
0.99997
1.4x110*
1.67x10
3.02x10
0.50
PINN
0.99973
5.23x10'
5.92x10*
1.04110*
1392.00
Problem 8
LS-SVMs
0,99999
7.14x10
8.93x10
2.86x10*
20.80
NLS-SVMs
0.99999
5.66x10
7.05x10
2.23x10*
0.44
PINN
0,99999
80
9.96x10
1.89x10*
776.00
Problem 9
LS-SVMs
0,99999
1.60x10
1.92x10
6.20x10
20.20
NLS-SVMs
0,99999
5.74x10
6.92x10
2.21x10*
PINN
0.98291
3.12x10
4.07x10*
8.62x10
636.00
Problem 10
LS-SVMs
0,99999
9,40x10
1.34x10
3.99x10*
24.70
NLS-SVMs
0.99999
4.29x10
6,06x10*
3.09x10*
0.49
PINN
0.99999
1,62x10
1.89x10
3.32x10
817.00
Problem 11
LS SV-S
0.99999
2.81x10
3.48x10*
116x107
16.30
NLS-SVM
0.99999
2.36x10
2.97x10
9.43x10
0.45
PINN
0,99999
2.28x10
2.29x10'
2.60x10*
969.00
Problem 12
LS-SVMs
0,99999
1.56x10
1,93x10*
5.91x10
16.70
NLS-SVMs
0,99999
3.15x10
3.94x10
1.31X10
0.46
PINN
0.99987
4,22x10
5.72x10*
1.05x10*
1475.00
Problem 13
LS-SVMs
0.99999
2.77x10
3.42x107
9.56x10
33.60
NLS-SVMs
0.99999
3.881107
4.36x10
1.45x10
0.48
PINN
0.99999
8.22x10
8.77x10
1.34x10
112.40
Problem 14
LS-SVMs
0.99999
7.38x10
9.17x107
2.02x10
36.00
NLS-SVMs
0.99999
1.17x10
1,57x10
3.85x10
1.39
PINN
0.99774
2.42x10*
2.62x10'
3.62x10
230.00
Problem 15
LS-SVMs
0,99999
2.14x10*
2.66x10*
6.68x10
85.40
NLS-SVMs
0,99999
1.68x10"
1.98x10"
3.46x10²
3.25
PINN
0.99812
1.24x10
1.35x10
L
2699.00
Problem 16
LS-SVMs
0,99999
1.50x10
1.88x104
6.84x10
24.30
NLS-SVMs
0.99999
1.72x10
2.06x10
5.07x10
0.42
Table 12: Performance evaluation of NLS-SVMs for long-time integration of ODEs: accuracy and computational efficiency
ODEs
Model
R²
MAE
RMSE
11 عl
Time/s
Problem
NLS--VMs
0.99999
1.38x10
4.23x10
4.40x10
34.80
Problem 3
NLS-SVMs
0,99999
5.95x10
1.1x100'
1.76x10
24.20
Problem 5
NLS--VMs
0.99999
9.96x10
2.36x10
2.70x10
413
Problem 7
NLS-SVMs
0.99911
1.59x10*
4.01x10"
2.36x10
299.00
Problem 8
NLS-SVMs
0.99999
4
5.08x10
1.71x10
24.70
Problem 9
NLSSSVMs
0.99996
2.
6.31x10
6.33x10*
459.00
Problem 10
NLS-SVMs
0.99999
2.72x10
5,47x10
3,49x10*
460.00
Problem 11
NLS-SVMs
0.99981
7.85x10*
9.62x10
1.55x10
372.00
Problem 12
NLS-SVMs
0.99709
6.55x10
7.27x10
1.05x10'
379.00
36

---

B.8
Performance improvements
Table 13: Benchmarking NLS-SVMs against PINN and LS-SVMs on ODE problems: accuracy and speed
ODEs
Model
AR²
AMAE/%
ARMSE/%
A lly
Vll /%
Speedup
NLS-SVM: PINN
+0,00000
3.66 10
4.1.0 10-1
48.01 x 10
1222
NLS-SVMs vs LS-SVMs
+0,00000
3.1.3
10
1.775 10
+7.00 x 10
1 18
NLS-SVM: vs PINN
+0.00010
2.03 10
2.09 10
3.225 x 10
1 733
Pro2
SVMs
1074
1 1.92 10
1.668 10-4
1
NLS-SVMs vs PINN
+0,25088
1.82 10*
2.88 10"
11.03 x 10
2169
Pro3
NLS SVMs vs LS SVMs
+0,00000
1 4.4.0 x 10
4.660 x 10-7
1.4.00 x 10-7
1 3355
NLS-SVMs vs PINN
+0,05445
1..34
7.20
11.81 x
4.49
Pro1
NLS-SVMs LS-SVMs
+0,00000
1.911 10
1
10
11.29 10
373
SVMs vs PINN
0,00001
4 1.51 x 10-4
1 3.00 x 10
5.40 x 10
171
Pro5
NLS-SVMs vs LS-SVMs
+0,00000
1.000 10
4 4.000 10 3
1.700 x 10
t 10
10" 4.44
1.66 10"
4.3.40 10"
f 190
Pro6
5LS-SVMs SSSSVMs
0.00019
1 4.99
1.3.7 10
Pro7
NLS-SVMs vs PINN
+0.08548
J. 6.41 10"
4 7.89 x 10"
1.66 x 10"
1 2012
NLS- va SV-
+0,000000
0.00
10
1 1.00 10-3
0.00 x
1
NLS-SVMs
5.17
10"
Pro8
NVS-SVM: LS S-SVMs
+0,000000
10 4.95
1 6.16 10-3
1.94 x 107
1 4747
NLS-SVMs vs PINN
+0.00000
2.92
10
3.3.4 x 10
1.3.20 x 10
1 1651
Pro9
NLS SVMs va L.S. SVMs
+0,000000
1 4.14
1 5.00 10
1.59 x 10
1 43
NLS-SVMs PINN
+0.01708
3.112
10"
4.00 10"
J.6.2 10"
1298
Pro10
NLS-SVMs vs LS-SVMs
+0.00000
1 3.35
10
1 4.72 10
9.9.0 x 10
t 50
Proll
NLS SVMs vs PINN
+0,000000
1.62
10"3
J 1.90
10
4.332 x 10-2
1 1816
4
1.07 10
f
Pro12
NLS SVMs PINN
+0,00000
$ 2.28
x 10-3
4 2.29
10
4.2.0 x 10-3
f 2107
S vs I.S. SVMs
+0,00000
1.25
10
4 1.54
10
4.460 x 10
1.36
NLS SVMs vs PINN
+0.00012
$ 44.22
10
5.7.2
10"
41.05 x 10"
1 3073
Pro13
SVSs
1 7.10
10
1 99.40 10-6
1.9.4 10
f
Pro14
NLS-SVMs vs PINN
+0.00000
4 7.05
10
4 7.20 10
499.5 x 10-4
8181
NLS-SVM, SVMs vs LS LS-SVMs
+0.000000
4.32 x 10-3
1 6.53 x 10-5
1.83 x 10-1
1 26
Pro15
NLS-SVMs va PINN
+0.00225
22422 10
$ 2.62
3.662 x
t
LS-SVMs
1.47
1075
1
10-5
43.22 10
1.26
Pro16
NLSS-VMs vs PINN
+0.00187
4 1.24
10"
1.35 16"
1.779 x 10"
1 6426
NLS-SVMs vs LS-SVMs
+0.00000
1.33
10
4 1.64 x 10
6.33 x 10
1.58
37


Reviews:
[{"review_id": "hK4N40YpqG", "reviewer": "Reviewer_AWxF", "summary": "The authors present an efficient method to solve ordinary differential equations by using Nystrom approximation. There are several machine learning based methods for solving ODEs, and least square support vector machine (LS-SVM) is one of most successful among these ML methods. However, LS-SVM reformulate the ODE problem to a optimization problem, which involves solving a matrix of dimension n x n, where n is the number of discrete time steps. This incurs a O(n^3) time complexity, which is prohibitive for long time ODE solution/simulation. The authors use Nystrom method, which is a low-rank matrix approximation to reduce the time complexity to O(m^3), where m<<n is the number of subsampled points. This leads to 10-3000 times faster computation with little compromise on accuracy. Experimental results are provided to support the authors’ claim.", "questions": "I am curious why the proposed method is particular for ODEs. Reducing the O(n^3) complexity will be super helpful for many applications, what makes the proposed method only specific for solving ODEs? Or maybe its not, in that case could you clarify where your method will be a good fit in addition to solving ODEs?", "limitations": "Mentioned in strength and weakness section", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength:\nThe motivation for the authors work is explained well, detailed experimental results are provided. \n\nWeakness:\n1. The problem of computational complexity needs to be discussed in more detail. Currently, it is unclear where the number of equations n and O(n^3) complexity comes from. My understanding is that due to finite differences applied on n discrete points in time, there are n variables, which all need to be solved for, along with the p initial/boundary conditions. If this is the case, add a small description about this for readers unfamiliar with ODEs and numeric solvers.\n2. For a paper with Nystrom in the title, there could be more details about the Nystrom methods in the main manuscript. The authors at present put the details in the supplementary, however it is difficult to understand the methodology without an intuitive/conceptual idea about the Nystrom method. Maybe give small overview in a few sentences in section 3, before 3.1\n3. In equation 15, what are f_k, f_p?\n4. Is there a missing summation in equation 16? Summation over i maybe\n5. The authors can think about rearranging their experimental results. As the major claim of the paper is computational efficiency, charts and tables in the main manuscript should support that claim. Currently I think those are in the supplementary. From figure 1,2 the authors can select a subset of ODEs to show that the proposed method achieves comparable accuracy, while using the rest of the space to show evidence for 10-3000 times computation efficiency. \n6. The figures legends at present need to be updated/enlarged, it is almost impossible to read at present, figure 1,2. Figure captions should also be updated, highlight the key takeaway in the caption, for example: The results show that the proposed method achieved comparable accuracy with baseline methods …", "was_revised": true, "final_justification": "The authors responded in details with all reviewer suggestions/comments/questions for clarification. The updated manuscript will be much easier for readers to understand.", "rebuttal": "Thank you very much for your time and effort in reviewing our manuscript. We sincerely appreciate your insightful comments and constructive suggestions, which have been invaluable in improving the quality and clarity of our work. We have carefully considered all the points raised and have revised the manuscript accordingly. Below, we provide a point-by-point response to each Weakness/Question.\n\nThis document is organized in the following manner: Each comment is numbered and noted in bold, and the response is written directly afterward, in normal font, with the quoted text directly from the revised manuscript or revised sections noted in \"...\".\n\n**Response to Weakness 1 (Computational Complexity Clarification):** We appreciate this insightful comment. We agree with the reviewer’s understanding. \"The parameter $n$ represents the number of discrete time points introduced by the finite-difference discretization of the ODE. This leads to solving for $n$ variables alongside the $p$ initial/boundary conditions.\" The above concise explanation has been added to \"Subsection 2.2\" to clarify the origin of the $O(n³)$ complexity for readers less familiar with ODE solvers.\n\n\n**Response to Weakness 2 (Nyström Method Overview):** Thank you for this valuable suggestion. As suggested, a brief conceptual overview of the Nyström method has been added at the beginning of \"Section 3 (before subsection 3.1)\". This overview succinctly introduces its core principle: \"it constructs a low-rank kernel matrix approximation by strategically selecting a subset of $n$ discrete time points, enabling efficient solution of a system of $n$ equations.\"\n\n    \n**Response to Weakness 3 (Definition of $f_k, f_p$ in Eq. (15)):** Thank you for your careful review.In equation 15, *f_k* is a function of time $t$ defined in a general $p$-th order ODE represented by \"equation 14\" for k = 1, ..., p, and *f_p* coresponds to the function when k = p.  To address the reviewer's comment, the authors have corrected equation 14 to clarify by changing *i* to *k* \"($y^{(p)}(t)-\\sum_{k=1}^{p} f_{k}(t) y^{(p-k)}(t)=r(t) ,\\quad t \\in[t_1, t_n] , \\quad y(t_1)=v_{1},\\quad y^{(k-1)}(t_1)=v_{k}, \\quad k=2, \\ldots, p $)\".\n\n\n**Response to Weakness 4 (Missing Summation in Eq. (16)):** Thank you for your careful review. We have carefully checked and confirm that no summation symbol is missing in Equation 16. This equation describes the constraint condition at each time-discrete point $i$. The index $i$ iterates implicitly over all discrete points ($i = 1, ..., n$), representing the discretized form of the ODE equation itself at each point.In the revised manuscript, the authors have added text to clarify. The added text is: \"the constraint conditions include $n$ time-discrete equations where each represents the discretized form of the ODE equation itself at each time point\". \n\n\n**Response to Weakness 5 (Reorganization of Results):** We agree with the reviewer's suggestion that computational efficiency is the primary claim. To emphasize this, \"Table 5 (previously in the Appendix) has been moved directly into the main results section (Section 4)\". This table quantifies the significant computational speedups (10–3000 times) under comparable accuracy regimes for key benchmarks. Selected ODE cases demonstrating comparable accuracy to baselines remain in \"Figure 1\", while efficiency results now dominate the primary presentation.\n\n\n**Response to Weakness 6 (Figure Legends and Captions):** Thank you for your suggestion. The figures and captions have been comprehensively revised:\n\n- *Legends*: \"Legends in Figures 1 and 2 have been significantly enlarged and redesigned  for immediate readability across all elements (lines, labels)\".\n- *Captions*: \"All figure captions now explicitly state the key conclusion\". \n\n\n**Response to Question (Method Applicability).** We appreciate this insightful question regarding the applicability of our method. Please allow us to clarify its scope and extensions: \n\n- *Current Focus on ODEs - Fundamental Constraint.*\nThe method's design in this work explicitly incorporates ODE equations as constraint conditions (e.g., \"Eq.7 and Eq.14\"). This specific formulation requires the target function to satisfy the ODE form. Therefore, the current implementation is indeed specialized for ODE problems where such constraints define the solution space. \n- *Framework Transferability - Beyond Ordinary Differential Equations.*\nWhile the constraint structure is currently ODE-specific, the core computational framework—low-rank approximation via Nyström subsampling—is universally applicable. For non-ODE problems, extending the method requires: Replace ODE constraints with domain-specific constraints and Maintain efficient kernel approximation framework. \n\nThe following text is added in the revised mansucript to address the reviewer's concern. The added text is: \"While demonstrated on ODEs, the Nyström acceleration framework can extend to arbitrary optimization problems by reformulating constraints. This includes large-scale PDEs (replacing differential operators) and ML training (using optimality conditions), subject to low-rank kernel structure\"."}, {"review_id": "V8u7ARhUGK", "reviewer": "Reviewer_EPWB", "summary": "This paper proposes a novel method for solving linear and nonlinear ordinary differential equations (ODEs) based on kernel-based approaches. Although kernel-based strategies have been employed to solve ODEs, the computational complexity of such solvers scales as $O((n+p)^3)$ where $n$ is the number of discretization points and $p$ is the order of the ODE. This makes them impractical for settings where the simulation is long-term or has a fine-grained temporal resolution. To address this challenge, authors propose a new Nyström-accelerated LS-SVMs framework that uses the Nyström method to reduce the computational complexity to $O((m+p)^3)$ where $m<<n$ is the number of subsampled landmark points. The authors test the new approach empirically and report huge speed-ups compared to traditional LS-SVMs and PINNs.", "questions": "- How does this approach compare to Runge-Kutta schemes?\n- What are some of the technical challenges of applying the Nyström method to this setting?\n\nIf those two points are clarified, I'll be happy to consider increasing my score. Currently, it is not yet sufficiently clear to me how impactful and novel the proposed method is.", "limitations": "Yes", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Strengths:\n- Important problem - quick ODE solving can have a huge impact on many areas\n- Clearly-written\n- Empirically validated\n- Huge reported speed-ups on the tested datasets with no or insignificant increase in error\n\nWeaknesses:\n- My main weakness concerns lack of comparison with standard approaches such as Runge-Kutta schemes. Is the approach also quicker or more accurate than those? Are there settings where kernel-based methods need to be used because those schemes are too inaccurate or computationally intensive? The paper mentions that \"these methods struggle with high-dimensional systems, noisy data assimilation, and adaptive mesh refinement\", but I could not find any reference or empirical result that would show it. \n- I think the technical novelty could be better emphasized. Was it challenging to apply Nyström method to LS-SVMs for ODE learning, or is it a straightforward combination of previous ideas? \n- Minor:\n\t- Typo in line 225. I think it should be \"second-order\".", "was_revised": true, "final_justification": "The authors have answered my questions and provided extensive empirical comparisons with current methods. In light of this additional evidence, I have increased my score to 5, and I am happy to see this paper accepted. I understand that some reviewers are concerned about technical novelty. However, I believe the solution's empirical performance makes up for it sufficiently.", "rebuttal": "Thank you very much for your time and effort in reviewing our manuscript. We sincerely appreciate your insightful comments and constructive suggestions, which have been invaluable in improving the quality and clarity of our work. We have carefully considered all the points raised and have revised the manuscript accordingly. Below, we provide a point-by-point response to each Weakness/Question.\n\nThis document is organized in the following manner: Each comment is numbered and noted in bold, and the response is written directly afterward, in normal font,with the quoted text directly from the revised manuscript or revised sections noted in \"...\".\n\n**Response to Comment 1: Lack of comparison with standard approaches (e.g., Runge-Kutta) and Is it quicker or more accurate?** We sincerely thank the reviewer for emphasizing the importance of comparing our approach with established numerical methods like Runge-Kutta. To comprehensively address this point, we have performed extensive new numerical experiments comparing the proposed efficient kernel-based method against the classical 4th-order Runge-Kutta method (RK4), implemented using variable time step sizes. These comparisons were conducted across all 15 benchmark ODE problems presented in the original manuscript. Key results are summarized in \"revised Figures 1-2 and Tables 1-5 in Section 4 and Appendiices of the updated manuscript\".\n\nThe results reveal a nuanced performance landscape dependent on the integration step size. The following text regarding the performance comparison of RK4 and proposed method is added in \"Section 4\". The added text is: \"For large step sizes (e.g., time step $\\geq$ 0.03 s), RK4 demonstrated computational efficiency advantages over proposed method (typically 3-10 faster computation times). However, this efficiency came at a significant cost in accuracy. RK4 solutions often suffered from low precision (RMSE degraded by ~90% on average compared to proposed method) and occasional divergence. For small step sizes (e.g., time step $<$ 0.03 s), Both RK4 and KBM achieved comparable high accuracy, approaching the analytical solution. But their computational efficiency showed strong dependence on integration duration. For short-duration problems, RK4 could be slightly more efficient.For long-duration problems, KBM demonstrated superior efficiency compared to RK4 with small time step size. Moreover, the proposed approach maintained consistently high accuracy across small-large step sizes (0.001s - 0.1s) with minimal variation in the solution precision. This highlights that the proposed method exhibits significantly lower sensitivity to step size selection\".\n\nTherefore, our proposed method is not universally quicker or more accurate than RK4, but offers distinct advantages in scenarios demanding robust accuracy at moderate-to-large step sizes and/or efficiency in long-duration integrations.\n\n**Response to Comment 2: Are there settings where kernel-based methods need to be used because those schemes are too inaccurate or computationally intensive?** Yes, the experiments conducted in response to Comment 1 identify specific settings where proposed efficient kernel-based methods offer crucial advantages over standard schemes like RK4:\n\n- *Scenarios requiring large integration step sizes:* In fields like Earthquake Engineering, seismic ground motion input data is often recorded at coarse time intervals (e.g., $\\Delta t$ = 0.01-0.05s or larger). Using RK4 with a step size $h$ dictated by the input data resolution ($h = \\Delta t$) often leads to inadequate accuracy or instability, as our results in Section 4 demonstrate (h $\\geq$ 0.03 s). The proposed method’s robust performance at these larger step sizes ($h$ = 0.01s - 0.1s) makes it essential for obtaining reliable solutions where standard schemes fail under coarse temporal resolutions.\n\n- *Long-duration integration problems:* When high solution accuracy mandates the use of small step sizes in RK4 for a long integration interval, the computational burden becomes obvious. Our results show the proposed method outperforms RK4 in efficiency for long simulations under these strict accuracy requirements.\n\n- *Problems sensitive to step size tuning:* The proposed method showed significantly lower sensitivity to the precise choice of time step size (unlike RK4, which demands careful tuning), which simplifies its application and reduces the risk of user-induced errors in scenarios where step size selection is non-trivial.\n\n**Response to Comment 3: Mentions shortcomings of standard approaches (e.g., Runge-Kutta) without reference or empirical result.** We appreciate the reviewer highlighting the need for verification. Our original statement regarding \"these methods struggle with high-dimensional systems, noisy data assimilation, and adaptive mesh refinement\" primarily referred to the inherent step size sensitivity and stability limitations of explicit fixed-step methods like RK4 in certain contexts. The extensive new comparison data (\"revised Section 4 and Appendiices, Figs 1-2, Tables 1-5\") provides explicit empirical evidence supporting this specific claim.\n\nIn the revised manuscript, we have corrected the statement to \"these methods face inherent challenges (inaccuracy, instability) under coarser temporal resolutions [3]\" in Section 1 and added \"new reference [3]\" below to clarify.\n\n\"[3] Butcher J C. Numerical methods for ordinary differential equations[M]. John Wiley & Sons, 2016.\"\n\n    \n**Response to Comment 4: Emphasis on technical novelty (Nyström + LS-SVM for ODE Learning).** We sincerely appreciate the reviewer's insightful comment regarding the novelty of our approach. We agree that the technical challenges and innovations deserve clearer emphasis in our manuscript. *The core novelty of this work lies in developing a novel computational framework specifically designed for ODE learning*. This framework uniquely integrates:\n\n- *A modified Least Squares Support Vector Machine (LS-SVM) formulation:* Adapted explicitly for the task of learning derivatives and satisfying differential constraints inherent in ODEs, distinct from standard LS-SVM for regression or classification.\n- *Embedded Nyström Approximation:* Strategically applied within this ODE-specific LS-SVM setting to achieve computational efficiency on potentially large-scale problems arising from dense temporal or state observations.\n\n*Key innovations focus on this novel framework and its successful application:*\n\n- *ODE-Specific Formulation:* We defined a tailored objective function within the LS-SVM paradigm, explicitly incorporating the ODE residual as a central constraint. This differs significantly from off-the-shelf LS-SVM usage and addresses the fundamental challenge of learning dynamics from data.\n- *Effective Computational Strategy:* While leveraging well-known concepts (Nyström, standard equidistant sampling for landmarks), the novelty stems from their seamless integration and application within this novel ODE-LS-SVM framework. This integrated design *enables significant computational gains* by circumventing large dense kernel matrices while focusing the approximation where needed. We demonstrate it *effectively solves the target ODE problem class.*\n- *Successful Application to ODE Learning:* To the best of our knowledge, this represents the *first successful implementation and empirical demonstration* of utilizing the Nyström-accelerated LS-SVM methodology specifically for learning and solving ODEs, proving its feasibility and advantages in robust handling.\n- *The proposed new method is scalable*—not only in terms of modelling flexibility but also regarding computational scalability for large-scale, long-duration applications.\n\nIn the revised manuscript, the authors have added text to clarify the novelty in Section 1.2. The added text is: \"This work introduces a novel computational framework specifically designed for ODE learning, combining two key components in an original way. First, it develops a modified Least Squares Support Vector Machine (LS-SVM) formulation explicitly tailored for learning derivatives and satisfying ODE constraints, distinct from standard regression/classification uses. Second, it incorporates an embedded Nyström approximation within this ODE-specific LS-SVM to achieve computational efficiency for large-scale problems with dense temporal/state observations. The primary innovations focus on three aspects. (1) **ODE-Specific Formulation**: A custom objective function within LS-SVM explicitly incorporates ODE residuals as constraints, fundamentally addressing the challenge of learning dynamics from data. (2) **Integrated Computational Strategy**: Though leveraging known concepts (Nyström, equidistant landmark sampling), the novelty lies in their seamless integration within the proposed framework. This design enables significant computational gains by avoiding large dense kernel matrices while maintaining approximation focus. (3) **First Successful Application**: To our best knowledge, this is the first empirically demonstrated implementation of Nyström-accelerated LS-SVM specifically for solving ODEs, proving its feasibility and robustness in handling dynamics\".\n\n**Response to Comment 5: Typo in line 225, should be \"second-order\".** Thank you for your thorough review. \"Line 225 has been corrected accordingly\". Additionally, the author has reviewed the entire text to prevent similar errors from occurring. Thank you once again for your valuable comments."}, {"review_id": "UnFMssn4Sz", "reviewer": "Reviewer_hduQ", "summary": "This paper introduces a Nyström-accelerated LS-SVM framework that reduces the $\\mathcal{O}(n^3)$ complexity of kernel-based ODE solvers to $\\mathcal{O}(m^3)$ using explicit primal-space formulations. The method enables efficient handling of linear/nonlinear ODEs through closed-form feature mapping and derivative computations. Experiments on fifteen benchmark problems show substantial speedups over LS-SVMs and PINNs with comparable or better accuracy.", "questions": "1. The use of uniform random sampling for Nyström landmarks may limit approximation quality in nonlinear ODEs. Have you considered evaluating more informed strategies (e.g., leverage score sampling), or include an ablation study to show how sampling affects stability and performance.\n2. The iterative optimization solver is critical in nonlinear ODEs, yet the paper provides no theoretical discussion on its convergence or stability behavior. A brief analysis of convergence would help build confidence in the solver’s reliability.\n3. The PINN used for comparison is based on a shallow architecture with identity activation, which may not reflect the performance of recent PINN variants. It would be helpful if the authors could justify the architectural choice and consider recent variants, especially those designed to reduce computation costs.\n4. The experiments focus primarily on first-/second-order ordinary differential equations, without evaluating higher-order systems or real-world dynamics. A discussion on the potential generalizability of the method to more complex systems would be valuable. Adding even one case study on a higher-order or real-world problem would significantly enhance the paper's applicability.\n5. Adding a schematic pipeline diagram and more visualization of the results would improve clarity. In addition, the notation could be made more consistent throughout the paper, for example, symbols like $\\varphi_k^p $ and $f_k$​ appear in both bold and non-bold forms across different equations and lines (e.g., Eq.(13), Eq.(15) , Eq.(16) ,Eq.(17), Eq.(18) and line 177), which may confuse readers.\n6. While the paper adopts an RBF kernel, it does not explain why this choice was made or how it compares with other kernels. It is suggested that author includes an empirical comparison with alternative kernels, along with a basic sensitivity analysis for RBF kernel parameters.\n7. Since the method involves iterative optimization, it would be important to clarify whether the results in Appendix A.7 and A.8 are based on single runs or averages over multiple trials. The definition and measurement procedure for the reported \"Time\" metric should also be clearly described.\nIf the authors address these concerns or provide convincing clarification, I would be open to revisiting my final evaluation.", "limitations": "yes", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. By introducing Nyström approximation in the primal space, the method effectively reduces the computational cost of LS-SVM-based ODE solvers from $\\mathcal{O}(n^3)$ to $\\mathcal{O}(m^3)$, improving scalability for large time steps and high-resolution settings.\n2. The paper derives explicit expressions for the feature mapping $\\varphi(t)$ and its derivatives, enabling model training without relying on the kernel trick.\n3. The paper is well-structured and clearly written, and the open-sourced code facilitates reproducibility and follow-up research.\n\nWeaknesses:\n1. The method employs simple uniform random sampling to select landmarks, which may lead to suboptimal approximation quality in highly nonlinear ODEs. The choice of sampling strategy could impact robustness and generalization.\n2. While the paper systematically integrates Nyström approximation into LS-SVM-based ODE solvers, it lacks theoretical analysis on subspace stability and solver convergence in nonlinear settings.\n3. The PINN baseline used in comparison adopts a shallow architecture with identity activation, without benchmarking against more recent improvements in the PINN literature.\n4. The experiments mainly focus on first/second-order ODEs, with no discussion of higher-order or real-world systems, limiting the evaluation scope.\n5. Mathematical notation could be further clarified, and visual explanations (e.g., diagrams) are insufficient. Specific concerns are outlined in the Questions section.", "was_revised": true, "final_justification": "I still think the work lacks theoretical analysis and the empirical comparison is weak. Authors state that the convergence behavior of their method follows standard Newton-method, then the theoretical contribution of this work is not significant. As for the experiments, the original baselines are quite weak. Authors state that they added a new recent PINN variant for comparison, but it is not convincing enough without detailed results.", "rebuttal": "Thank you very much for your time and effort in reviewing our manuscript. We sincerely appreciate your insightful comments and constructive suggestions, which have been invaluable in improving the quality and clarity of our work. We have carefully considered all the points raised and have revised the manuscript accordingly. Below, we provide a point-by-point response to each Weakness/Question.\n\nThis document is organized in the following manner: Each comment is numbered and noted in bold, and the response is written directly afterward, in normal font,with the quoted text directly from the revised manuscript or revised sections noted in \"...\".\n\n**Response to question 1 (Uniform sampling for Nyström landmarks).** W appreciate the reviewer's valuable suggestion. We would like to clarify that we performed evenly-spaced sampling to select the subset from $n$ discrete time points in our original manuscript. In the revised manuscript, we conducted extensive ablation studies on sampling strategies ( \"Section 4, Table 6\") for all 15 benchmark ODE problems plus a new fourth-order ODE. The added text to describe the ablation studies on sampling strategies is: \"While leverage score sampling improved mean accuracy by about 8% over uniform sampling, it increased computation time by about 120%. Uniform sampling exhibited high instability (higher variance across 10 runs). In contrast, Evenly-spaced sampling achieved comparable accuracy to leverage score sampling (mean error < 0.3%) with about 65% faster computation time and near-zero variance. This deterministic approach ensures reproducibility without sacrificing accuracy. Table 6 presents the mean and standard deviation of R2, RMSE and MAE for these three sampling methods across 10 runs\".\n\n\n**Response to question 2 (Convergence analysis of iterative solver).** We sincerely thank the reviewer for emphasizing the importance of theoretical convergence analysis. Our solver employs Newton's method to solve the nonlinear system derived from ODEs (e.g., Equation 21). The convergence behavior follows standard Newton-method theory: \"sufficient conditions for convergence are satisfied when regularization parameter $\\gamma$ and RBF kernel parameter $\\sigma^2$ ensure positive definiteness of the coefficient matrix in equation 21\". Additionally, extensive experiments confirm robust convergence. Across all tested systems (including stiff problems in Section 4), the solver consistently achieved residuals < 0.13% within 50-500 iterations without divergence (\"Table 3-4\"). This empirically confirms reliability for our problem class. We also would like to emphasize that our numerical results validate the solver’s robust and efficient convergence behaviour across the range of problems studied herein.In the revised manuscript, \"we have also supplemented the number of iterations of the proposed method in solving nonlinear ODEs in Table 1\". \n\n\n**Response to question 3 (PINN baseline comparison).** Thank you for your suggestion. To further demonstrate the accuracy and efficiency of the proposed method, the authors have supplemented a comparison of the performance of recent PINN variant (Fourier Feature PINNs) used to solve ODEs. The computational efficiency of the method proposed in this study is still greatly improved compared to them, with an increase of 10-1000 times, and its prediction accuracy is comparable, verifying the conclusions proposed in this paper. Further, we also tested a deep architecture with 8 layers and tanh activation function, the accuracy is not significantly improved compared to the PINNs employed in this paper, but the computational time increases significantly.\n\n**Response to question 4 (Limited experimental scope).** Thank you for your valuable suggestion. In the revised manuscript, \"we have added an example of a higher-order ODE ($\\frac{d^4y}{dt^4} =120t,t\\in [-1,1],y(-1)=1,y{}' (-1)=5,y(1)=3,y{}' (1)=5$) in this paper (Section 4)\". Its computational accuracy is comparable to that of PINN, but its computational efficiency has been significantly improved by a factor of about 1000.The results of this example further demonstrate the performance of the proposed method.\n\n**Response to question 5 (Visualization and notation clarity).** Thank you for your careful review and valuable suggestion. In the revised manuscript, \"we have reviewed and standardized all symbols throughout the text. Additionally, a flowchart illustrating the proposed method has been created to facilitate readers' comprehension and understanding\".\n\n**Response to question 6 (RBF kernel justification).** Thank you for your valuable suggestions. Within the framework of the proposed approach, the choice of kernel function directly affects the model’s ability to fit nonlinear relationships. This paper selects the Radial Basis Function (RBF) kernel primarily based on the following theoretical and empirical grounds: \n    \n- *Universality and performance advantages*— \"the RBF kernel can map samples into an infinite-dimensional feature space, possessing universal approximation capability and being insensitive to data scaling. Compared to the linear kernel (which is only suitable for linearly separable problems)and the polynomial kernel (which tends to overfit at higher degrees), the RBF kernel demonstrates more robust performance on complex nonlinear problems\". \n\n- *The performance of the RBF kernel depends on the kernel parameter*. \"We have conducted a systematic grid search for this hyperparameter. The results of the sensitivity analysis have be included in the manuscript\".\n\n \n**Response to question.7(Reproducibility metrics).** Thank you for your careful review. The results in Appendices A.7 and A.8 are based on a single run, as our method (evenly-spaced sampling) does not alter results under multiple runs when regularization and RBF kernel parameters are fixed.  \"Time\" metrics exclusively measure **wall-clock solution time** from algorithm initiation to end, computed via Python's time.perf_counter() on our test platform (\"Appendix A.8\"), and the authors have clarified this in the paper (\"Appendix A7 and A8\")."}, {"review_id": "24f2Ws60om", "reviewer": "Reviewer_uEyG", "summary": "The paper focuses on solving ordinary differential equations using kernel methods, with the focus being on the use of Least-Squares SVMs (LS-SVMs). Similar to the classic supervised learning setting, LS-SVMs can struggle on large-scale tasks owing to their cubic scaling. In the ODE setting, the cost of LS-SVMs, scales cubically with $n$, the number of timepoints, making the method prohibitive as the number of timepoints becomes large.\n\nTo address this challenge, the paper adapts the classic Nyström method to their setting, and show that it reduces the cost to $O(m^3)$, where $m$ is the number of landmarks selected. Numerical experiments on several benchmark ODEs, shows the proposed method exhibits good performance relative to the LS-SVM baseline and PINNs.", "questions": "I provided suggestions for improvement in my reply under strengths and weaknesses.", "limitations": "Yes.", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 1, "strengths_and_weaknesses": "**Strengths**\n\nThe main strength of the paper is that the proposed method is a lot faster than the LS-SVMs baseline without compromising performance too much. The method also outperforms the basic PINNs baseline on certain problems.\n\n**Weaknesses**\n\nThe paper has many weaknesses, which I summarize below, that make it inappropriate for publication at Neurips.\n\n*Lack of theoretical support*\n\n    The paper does not perform any theoretical analysis that would help validate the empirical \n    effectiveness of the proposed algorithm.\n    This is a significant weakness in my view, in light of the lack of novelty of the approach (see below), \n    and the fact that most papers that use the Nyström provide strong theoretical guarantees on its performance, i.e. [1], [2], [3], [4].\n\n*Lack of Novelty*\n\n    The Nyström method is a classic approach for reducing the cost of kernel methods, from SVMs to \n    kernel ridge regression. \n    The main contribution of the paper seems to just be applying this technique to the setting of solving ODEs \n    and verifying it works well empirically.\n    As mentioned in my previous point, the work introduces no new theory for this specific application that justifies the effectiveness of the method.  \n    While the use of Nyström might be new in the context of solving ODEs, by itself this is not enough to warrant \n    acceptance at Neurips. \n    Especially since prior works in Neurips such as [5], that use the Nyström method to accelerate SciML come with strong theoretical support.\n\n*Weak baselines* \n\n    The baselines used for empirically evaluating the method are quite weak to me. \n    \n    The baseline LS-SVMs method is from 2012. \n    The main bottleneck of this approach is the need to solve kernel linear systems in Newton's method. \n    Since this method was introduced, the sketching and randomized linear algebra literature have developed efficient \n    methods for solving large kernel systems [6,7,8,9], as well as Newton's method for solving non-linear equations [9]. \n    Rather than solving the system that arises in solving the LS-SVM method using a direct method \n    which naively costs $O(n^3)$, \n    a more reasonable comparison would be to use one of the randomized iterative methods I just referenced to solve the \n    Newton system, which cost at most $O(n^2)$. \n\n    As for the PINNs baseline,  the vanilla PINN is very rarely used in practice, as it is known to exhibit failure modes [8, 9], which has spawned the development of much more effective variants. \n    The empirics would be more convincing if you compared with a more sophisticated PINN architecture, \n    ideally one tailored to solving ODEs or a specific class of ODE.\n\n*Related work*\n\n    The Nyström method has a rich history and has played a significant role \n    in making large-scale kernel computations possible, and the present paper is missing quite a few references. \n    Works 1-5, as well as 8-9 below should all be cited as examples of how the Nyström method has accelerated \n    kernel learning.\n\n**Summary**\n\n    In its present form this submission is not a good fit for Neurips. \n    The submission develops no new theory and the numerical baselines are too weak to demonstrate whether or not the proposed method represents a real practical advance. \n\n    If the authors do not plan to develop more theoretical support for their method, \n    I would recommend instead that they improve the numerical baselines and submit to a \n    computational science/engineering journal as it would be a better fit.\n\n**References**\n\n[1] Rudi, A., Camoriano, R. and Rosasco, L., 2015. Less is more: Nyström computational regularization. Advances in Neural Information Processing Systems, 28.\n\n[2] Rudi, A., Carratino, L. and Rosasco, L., 2017. Falkon: An optimal large scale kernel method. Advances in Neural Information Processing Systems, 30.\n\n[3] Marteau-Ferey, U., Bach, F., & Rudi, A. (2019). Globally convergent newton methods for ill-conditioned generalized self-concordant losses. Advances in Neural Information Processing Systems, 32.\n\n[4] Della Vecchia, A., De Vito, E., Mourtada, J. and Rosasco, L., 2024. The Nyström method for convex loss functions. Journal of Machine Learning Research, 25(360), pp.1-60.\n\n[5] Meanti, G., Chatalic, A., Kostic, V., Novelli, P., Pontil, M., & Rosasco, L. (2023). Estimating Koopman operators with sketching to provably learn large scale dynamical systems. Advances in Neural Information Processing Systems, 36, 77242-77276. \n\n[6] Tu, Stephen, Shivaram Venkataraman, Ashia C. Wilson, Alex Gittens, Michael I. Jordan, and Benjamin Recht. \"Breaking locality accelerates block Gauss-Seidel.\" In International Conference on Machine Learning, pp. 3482-3491. PMLR, 2017.\n\n[7] Gazagnadou, Nidham, Mark Ibrahim, and Robert M. Gower. \"RidgeSketch: a fast sketching based solver for large scale ridge regression.\" SIAM Journal on Matrix Analysis and Applications 43, no. 3 (2022): 1440-1468.\n\n[8] Frangella, Zachary, Joel A. Tropp, and Madeleine Udell. \"Randomized nyström preconditioning.\" SIAM Journal on Matrix Analysis and Applications 44, no. 2 (2023): 718-752.\n\n[9] Díaz, Mateo, Ethan N. Epperly, Zachary Frangella, Joel A. Tropp, and Robert J. Webber. \"Robust, randomized preconditioning for kernel ridge regression.\" arXiv preprint arXiv:2304.12465 (2023).\n\n[10] Yuan, Rui, Alessandro Lazaric, and Robert M. Gower. \"Sketched Newton--Raphson.\" SIAM Journal on Optimization 32, no. 3 (2022): 1555-1583.\n\n[11] Krishnapriyan, A., Gholami, A., Zhe, S., Kirby, R. and Mahoney, M.W., 2021. Characterizing possible failure modes in physics-informed neural networks. Advances in Neural Information Processing Systems, 34, pp.26548-26560.\n\n[12] Rathore, Pratik, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. \"Challenges in training PINNs: A loss landscape perspective.\" arXiv preprint arXiv:2402.01868 (2024).", "was_revised": true, "final_justification": "After carefully reading the authors' latest rebuttal, I decided to increase my score to 4. \nMy reasoning primarily lies in the fact that the authors have now made clear the novelty of their approach. \nIn particular, it is clearer now that the paper is not just straightforwardly applying the Nyström method to the LS-SVM formulation from the prior work in 2012. \nI am also happy that the authors will add a convergence argument for their Newton solver in the final version.\nIn addition, they also addressed my concerns about the experiments by adding some more compelling baselines to show the advantage of their approach. \n\nI hope the authors will implement their proposed revisions, which include further references, clarification of novelty/relation to prior work, and additional experiments/theory, as this was crucial for convincing me to recommend acceptance.", "rebuttal": "Thank you very much for your time and effort in reviewing our manuscript. We sincerely appreciate your insightful comments and constructive suggestions, which have been invaluable in improving the quality and clarity of our work. We have carefully considered all the points raised and have revised the manuscript accordingly. Below, we provide a point-by-point response to each Weakness/Question.\n\nThis document is organized in the following manner: Each comment is numbered and noted in bold, and the response is written directly afterward, in normal font,with the quoted text directly from the revised manuscript or revised sections noted in \"...\".\n\n**Response to weakness 1 (Lack of theoretical support and novelty)** Thank you for your thorough review of this paper. We fully acknowledge the importance of the issues you raised regarding the lack of theoretical analysis and insufficient innovation. However, it is undeniable that the method we proposed is entirely characterised by theoretical innovation. To further emphasise this feature, we have significantly enhanced the theoretical depth and the explanation of contributions in the revised manuscript through the following work.\n\n### Theoretical exposition\n\n**NLS-SVMs for learning $p$-th Order Linear ODE**\n\n“Let us now consider the general $p$-th order ODE with initial value problem (IVP) in [7] of the following form:  \n$y^{(p)}(t)-\\sum_{k=1}^{p} f_{k}(t) y^{(p-k)}(t)=r(t) ,\\quad t \\in[t_1, t_n] , \\quad y(t_1)=v_{1},\\quad y^{(k-1)}(t_1)=v_{k}, \\quad k=2, \\ldots, p$\n\nThe approximate solution can be obtained by solving the following optimization problem: \n\n$$\n\\begin{aligned}\n\\underset{\\omega, b, e_i}{\\operatorname{minimize}}\\quad \n    & \\frac{1}{2} \\boldsymbol{\\omega}^{T} \\boldsymbol{\\omega}\n    +\\frac{\\gamma}{2} \\sum_{i=1}^{n} e_{i}^{2} \\\\\n\\text{s.t.}\\quad \n    & y^{(p)}(t_i) = \\boldsymbol{\\omega}^{T} \\boldsymbol{\\varphi}^{(p)}(t_{i})\n    =  \\omega^{T}\\left[\\sum_{k=1}^{p} f_{k}(t_{i})  \\varphi_{i}^{(p-k)}\\right] \n    + f_{p}(t_{i}) b + r(t_{i}) + e_{i}, \\quad i=2,\\ldots,N \\\\\n    & y(t_1) = \\boldsymbol{\\omega}^{T} \\boldsymbol{\\varphi}(t_{1}) + b = v_{1} \\\\\n    & y^{(\\ell-1)}(t_1) = \\boldsymbol{\\omega}^{T} \\boldsymbol{\\varphi}^{(\\ell-1)}(t_{1}) = v_{i}, \\quad \\ell=2,\\ldots,p\n\\end{aligned}\n$$\n\n***Theorem 1*** Given a positive definite kernel function $K$ : $R × R → R$ and a regularization parameter γ ∈ R+, the solution to optimization problem is given by the primal problem.(See Section3.1,Eq(18))\"\n\n\n**NLS-SVMs for Learning $p$-th Order Nonlinear ODE**\n\n“When the function is nonlinear, the $p$-th order ODE is nonlinear and written in the following form (initial value problems for illustration): \n\n$y^{(p)}=f(t, y^{(p-1)},\\cdots,y^{\\prime\\prime},y^{\\prime},y), \\quad y(t_1)=v_{0},\\cdots,y^{(p-1)}(t_1)= v_{p-1} \\quad t_1 \\leq t \\leq t_n$\n\nThen start by assuming the approximate solution to be of the form ${ \\hat{y}(t)}={\\omega } ^{T} \\varphi ({t} )+  b $. Additional unknowns $y_i$ are introduced to keep the constraints linear in $\\omega$. This yields the following nonlinear optimization problem.(See Section2.2, Eq8)\n\n***Theorem 2*** Given a positive definite kernel function $K$ : $R × R → R$ and a regularization parameter γ ∈ R+, the solution to optimization problem is given by the primal problem.(See Section3.2,Eq(21))\"\n\n\nThe mathematical proof for these two theoritical analyses is provided in Appendix A.1-A.6.\n\n**The aforementioned theoretical analysis has been updated in the revised manuscript.** The authors “present the theorem and the theoretical underpinnings for its proof”. This theory offers significant advantages in the context of solving Ordinary Differential Equations (ODEs). While classical Least Squares Support Vector Machine (LS-SVM) optimization algorithms have achieved substantial success in ODE applications, they inherently lack scalability. In contrast, **the method proposed in this paper extends it to large-scale applications.** This scalability constitutes the core contribution of the present research and is supported by a solid theoretical foundation.\n\nIt is worth noting that when learning a class of $p$-th order nonlinear ordinary differential equations that need to be solved, the changes in the ODE primarily affect the $\\frac{\\gamma}{2} \\Big( \n\\boldsymbol{\\omega}^{T} \\boldsymbol{\\varphi}^{(p)}(t_{i})- f(t_i,y_{t_i}^{(p-1)},...,y_{t_i}',y_{t_i}) \\Big)^{T}\\Big( \n\\boldsymbol{\\omega}^{T} \\boldsymbol{\\varphi}^{(p)}(t_{i})- f(t_i,y_{t_i}^{(p-1)},...,y_{t_i}',y_{t_i}) \\Big)$ function, as indicated by the Lagrangian loss function. Therefore, the Jacobian matrix needs to be updated based on the derivation of the four functions $\\frac{\\partial {Y_1} }{\\partial \\omega  }$ , $\\frac{\\partial Y_1}{\\partial y_i} $,$\\frac{\\partial {Y_{p+3}} }{\\partial \\omega  } $and $ \\frac{\\partial Y_{p+3}}{\\partial y_i}$. Therefore, in terms of algorithmic contribution, we propose **the first scalable Nyström framework** for ODEs. Our work is not merely an application of Nyström to ODEs, but rather **addresses the theoretical and algorithmic commonalities in the ODE solving process.**\n\nWe firmly believe that the improved theoretical exposition and the solutions to the specific challenges of ODEs have significantly enhanced the depth and originality of the paper. Thank you once again for prompting us to refine this work!\n\n**Response to weakness 2 (Weak baselines and nonlinear solver)**\n\nThank you for your suggestion. To further demonstrate the accuracy and efficiency of the proposed method, the authors have supplemented a comparison of the performance of recent PINN variant (Fourier Feature PINNs) used to solve ODEs. The computational efficiency of the method proposed in this study is still greatly improved compared to them, with an increase of 10-1000 times, and its prediction accuracy is comparable, verifying the conclusions proposed in this paper. Further, we also tested a deep architecture with 8 layers and tanh activation function, the accuracy is not significantly improved compared to the PINNs employed in this paper, but the computational time increases significantly.\n\nRegarding the nonlinear solver, the authors would like to clarify that the existing efficient method (e.g., Conjugate Gradient Method) cannot guarantee convergence. In the revised manuscript, the authors have conducted a theoritical analysis on convergence using Newton method.The detail is that \"sufficient conditions for convergence are satisfied when regularization parameter $\\gamma$ and RBF kernel parameter $\\sigma^2$ ensure positive definiteness of the coefficient matrix in equation 21\". Additionally, extensive experiments confirm robust convergence. Across all tested systems (including stiff problems in Section 4), the solver consistently achieved residuals < 0.13% within 50-500 iterations without divergence (Table 3-4). This empirically confirms reliability for our problem class. We also would like to emphasize that our numerical results validate the solver’s robust and efficient convergence behaviour across the range of problems studied herein.In the revised manuscript, we have also supplemented the number of iterations of the proposed method in solving nonlinear ODEs in Table 1.    \n\n**Response to weakness 3 (Related work)**\n\nThank you for your suggestion. The paper has been “supplemented with relevant literature on the Nyström method” for implementing large-scale kernel computations.\n\n[1] Rudi, A., Camoriano, R. and Rosasco, L., 2015. Less is more: Nyström computational regularization. Advances in Neural Information Processing Systems, 28.\n\n[2] Rudi, A., Carratino, L. and Rosasco, L., 2017. Falkon: An optimal large scale kernel method. Advances in Neural Information Processing Systems, 30.\n\n[3] Marteau-Ferey, U., Bach, F., & Rudi, A. (2019). Globally convergent newton methods for ill-conditioned generalized self-concordant losses. Advances in Neural Information Processing Systems, 32.\n\n[4] Della Vecchia, A., De Vito, E., Mourtada, J. and Rosasco, L., 2024. The Nyström method for convex loss functions. Journal of Machine Learning Research, 25(360), pp.1-60.\n\n[5] Meanti, G., Chatalic, A., Kostic, V., Novelli, P., Pontil, M., & Rosasco, L. (2023). Estimating Koopman operators with sketching to provably learn large scale dynamical systems. Advances in Neural Information Processing Systems, 36, 77242-77276.\n\n[6] Frangella, Zachary, Joel A. Tropp, and Madeleine Udell. \"Randomized nyström preconditioning.\" SIAM Journal on Matrix Analysis and Applications 44, no. 2 (2023): 718-752.\n\n[7] Díaz, Mateo, Ethan N. Epperly, Zachary Frangella, Joel A. Tropp, and Robert J. Webber. \"Robust, randomized preconditioning for kernel ridge regression.\" arXiv preprint arXiv:2304.12465 (2023)."}]

Author rebuttal:

Meta-review: This paper proposes a Nyström-accelerated LS-SVM framework for solving ODEs by reformulating the problem in the primal space and deriving explicit Nyström-based mappings and derivatives. This reduces the cost from $O(n^3)$ to $O(m^3)$. Experiments on multiple benchmarks show significant speed-ups  over baseline approaches with comparable accuracy.

Reviewers agree that the paper addresses an important problem and demonstrates strong empirical performance. The authors were very responsive, adding convergence analysis, ablation studies, and new baselines in their rebuttal. Two reviewers (AWxF, EPWB) recommend acceptance, one (uEyG) raised their score to borderline accept after rebuttal, and one (hduQ) remains skeptical, citing limited theoretical novelty and insufficiently strong baselines.

Overall, while the theoretical contribution remains debated, the method is practically useful, well-presented, and reproducible. Therefore, I recommend acceptance.
Decision: Accept (poster)