RandNet-Parareal: a time-parallel PDE solver using Random Neural Networks
Authors: Guglielmo Gattiglio
Lyudmila Grigoryeva
Massimiliano Tamborrino
Venue: NeurIPS 2024 poster (2024)
Keywords: Random Neural Networks
Extreme Learning Machines
Parareal
Parallel-in-time algorithms
Scalability
Ordinary and Partial Differential Equations
Primary area: machine_learning_for_physical_sciences

Abstract: Parallel-in-time (PinT) techniques have been proposed to solve systems of time-dependent differential equations by parallelizing the temporal domain. Among them, Parareal computes the solution sequentially using an inaccurate (fast) solver, and then ``corrects'' it using an accurate (slow) integrator that runs in  parallel across temporal subintervals. This work introduces RandNet-Parareal, a novel method to learn the discrepancy between the coarse and fine solutions using random neural networks (RandNets). RandNet-Parareal achieves speed gains up to x125 and x22 compared to the fine solver run serially and Parareal, respectively. Beyond theoretical guarantees of RandNets as universal approximators, these models are quick to train, allowing the PinT solution of partial differential equations on a spatial mesh of up to $10^5$ points with minimal overhead, dramatically increasing the scalability of existing PinT approaches. RandNet-Parareal's numerical performance is illustrated on systems of real-world significance, such as the viscous Burgers' equation, the Diffusion-Reaction equation, the two- and three-dimensional Brusselator, and the shallow water equation.

RandNet-Parareds a time-parallel PDE solver using Random Neural Networks
Guglielmo Gattiglio Department of Statistics University Warwick Coventry, CV4 7AL., UK
Guglielmo cawwiggooo@arwwwckkkkkkk........
Lyudmila Grigoryeva'
Faculty of Mathematics and Statistics University of St. Gallen
Rosenbergstrassee 20, CH-9000 St. Gallen, Switzerland
Lyudmila Grigoryeva@uniag....... ch
Massimiliano Tamborrino Department of Statistics University of Warwick Coventry, CV4 7AL, UK
Massimi liano. Tamborrino@warwicc. ac.ut
Abstract
Parallel-in-tin (PinT) techniques have been proposed to solve systems of time- dependent differential equations by parallelizing the temporal domain. Among them, Parareal computes the solution sequentially using an inaccurate (fast) solver, and then "corrects" it using accurate (slow) integrator that runs in parallel across temporal subintervals. This work introduces LandNet-PPrarea. learn the discrepancy between the coarse and fine solutions using random neural networks (RandNets). HandNet-PPararee achieves speed gains up to x125 and x22 compared the fine solver run serially and Parareal, respectively. Beyond theoreti- cal guarantees RandNets as universal approximators, these models are quick train, allowing the Pin'T solution of partial differential equations on spatial mesh up to 10t points with minimal overhead, dramatically increasing the scalability of existing PinT approaches. Sanddet-Parareal'' 's numerical performance is illus- trated on systems of real-worke significance. such as the viscous Burgers' equation, the CnftusionnRRecttii equation, the two- and three-dimensional Brusselator, and the shallow water equation.
novel method
Introduction
Parallel-in-time (PinT)) methods have been used to overcome the saturation of well-established spatial parallelism approaches for solving (prohibitively expensive) initial value problems (IVPs) for ordinary and partial differential equations (ODEs and PDEs), described by systems of d E N ODEs (and
"Honrrrry Associate Professor. Department of Statistics, University of Warwick, Coventry, CV4 7AL. UK. (Lyudmila. Grigoryyvaawwarvick.k ac. uk.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

similarly for PDEs)
du
dt = h))(t),t) on ante E twl, with u (to) = u. NeN,
(1)
where h : Fd x [x, - Fd is smooth multivariate function, [do ty] + Rd is the time dependent column vector solution, and u² E IR' is the initial value at to. PinT schemes are particularly important when the sequential application of an accurate numerical integrator & over [Lo, (v] is infeasible in a reasonable wallclock time. There are three general approaches for PinTT computation: parallel across-the- problee, parallel-across-the-ttr- and arrallel-acros--th--metet another classification is provided: multiple shooting, methods based on waveform relaxation and domain decomposition, multigrid approaches. and direct time- aaee-palllll methods. Parallel-across the- step methods, in which solutions at multiple time-grid points computed simultaneously. include Parareal (approximation of the derivative in the shooting method) (45), Parallel Full Approximation Scheme in Space and Time (PFASST) (multigrid method) [13 50 and Multigrid Reduction in Time (MGRIT) [14 [6] methods (see [9 for details). Among them, Parareal 445 has garnered popularity, with extensive theoretical analyses, improved versions, and empirical applications 557) This is due to its non-intrusive nature which allows seamless integration with arbitrary temporal and spatial discretizations, and to its successful performance across diverse fields, such as plasma physics [64] 66, 67], finance 4 [66. and weather modeling [59] 60). Limited theoretical results are available for MGRIT and PFASST. with a few extensions and empirical applications. Interestingly, combined analyses have shown equivalences between Parareal and MGRIT, and connections between MGRIT and PFASST. In Parareal, a coarse and fast solver L is run sequentially to obtain first approximation of the solution, which is then corrected by running fine (accurate) but slow integrator & in parallel across N temporal subintervals. This procedure is then iterated until convergence criterion is met after l S N iterations, leading to a speed-up compared to running f sequentially over the entire time interval. A recent advancement, GParareal [57], improves Parareal convergence rates (measured as k/N) by learning the discrepancy & 3 using Gaussian Processes (GPs). This method outperforms Parareal for low-dimensional ODEs and moderate number of computer cores N. However, the cubic cost (in the number of data points, roughly kN at iteration k)) of inverting the GP covariance matrix hinders its broader application. Subsequent research introduced nearest neighbors (nns) GParareal (nnGParareal) [21], enhancing GParareal's is scalability properties in both N and through data reduction. Significant computational gains were achieved by training the GP on small subset of nns, resulting in algorithm loglinear in the sample size. This allowed scaling its effectiveness up to systems with a few thousand ODEs, beyond which it loses its potential. Indeed, being based on the original GP framework, it uses costly hyperparameter optimization procedure that requires fitting one GP per ODE dimension.
In 17 55),
This study introduces RandNet-Parareal a new approach using random neural networks (RandNets) to learn the discrepancy F S. RandNets are family of Singe-hidden--ee feed-forward neural networks (NNs), where hidden layer weights are randomly sampled and fixed, and only the output (or readout) layer is subject to training. Compared to standard artificial NNs, RandNets are hence much simpler to train: the input data are fed through the network, the predictions observed, and the weights of the linear output (or readout) layer are obtained as minimizers of a penalized squared loss between the NN outputs and the training targets. Since this optimization problem admits a closed- form solution, no backpropagation is required, and the issues of vanishing and exploding gradients persisting for standard fully trainable NNs are therefore avoided. The literature on the topic is rich and somewhat fragmented, and different names are used for essentially the same model. RandNets are related to Random Feature Networks 66
49 62
63
65 and Reservoir Computing 24 26 25 27 28
61] 70 744. Some authors use the name refer to RandNets, while others use the term 7 for the same paradigm. RandNets show excellent empirical performance, and have been used in the context of mathematical finance 2
Random Fourier Features (RFFs) and kernel methods 41
Extreme Learning Machines (ELMs) 34
37
randomized or random NNs 55] 32 39
46
3
38
mathematical physics (52), electronic circuits [69], photonic
[7] and quantum systems 23 48
random deep splitting schemes [53], scientific computing [00
81) and have shown excellent
empirical performance in numerous further applications. Moreover, recent work [22 25
proves that RandNets are universal approximators within spaces of sufficiently regular functions, and provides explicit approximation error bounds, with these results generalized to large class of Bochner spaces in [52]. These contributions show that RandNets are a reliable machine learning paradigm with provable theoretical guarantees.
2

---

In this paper, we show that endowing Parareal with RandNets- based learning of 3 - G, the new proposed SandNet-Parareee algorithm. leads significantly improved scalability, convergence speed. and parallel performance with respect to nnGParareal, GParareal, and Parareal. This allows us to solve PDE systems on a fine mesh of up to 10 discretization points with negligible overhead, outperforming nnGaarreeal by two orders of magnitude and reducing its model cost by several orders.
Here, we compare the performance of Parareal, nnGParareal, and RandNet-Parareal on five in- creasingly complex systems, some of which are drawn from an extensive benchmark study of time- denendennddntt PDEs [75). These include the one-dimensional viscous Burgers' equation, the two- dimensional Diffusion-Reaccee formation 776,, the two- and daree-dimensiona Brusselator, known for its complex behavior, including oscillations, spatial patterns, and chaos, and the shallow water equations (SWEs). Derived from the compressible Navier-Stokes equations, the SWES are a system of hyperbolic PDEs exhibiting several types of real-world significance behaviors known to challenge numerical integrators, such as sharp shock formation dynamics, sensitive dependence on initial conditions, diverse boundary conditions, and spatial heterogeneity. Example applications include of tsunamis or flooding simulations.
equation, a challenging benchmark used to model biological pattern
We intentionally chose two hyperbolic equations (Burgers' and SWE) to challenge RandNet-PPaared on systems for which Parareal is known to struggle, with slow or non-convergent behavior [2] 3 g
18
72.. Previous works have developed hoe coarse solvers address Parareal's slow convergence
for Burgers' 7 40 68 [7]], and for SWE [1] 31 54 [73.. Here. we adopt a different strategy: by leveraging the generalization capabilities of RandNets within the Parareal algorithm, we enhance the performance of standard, off-the-shelff integration methods such as Runge- Kutta, obtaining speed gains up to x 15 and x22 compared to the accurate integrator I and Parareal, respectively. All experiments have been executed on Dell PowerEdge C6420 compute nodes each with 2 Intel Xeon Platinum 826 (Cascade Lake) 2.9 GHz 24-core processors, 48 cores and 192 GB DDR4-2933 RAM per node. To illustrate our proposed algorithm and facilitate code adoption, we provide step-by-step Jupyter notebook outlining RandNet-PParaeee Moreover, all simulation outcomes, including tables and figures, are fully reproducible and accompanied by the necessary Python code at https:/gitthb..com www/Paralll--iinnnnaccccaaaa Time- siine--Dfferemtilll Seafferentill-Equattionsss/// Reaiiaett/naneeeeeeeeeeeeee
It is well acknowledged that comparing PinT methods based on different working principles is extremely hard, with [55] representing a recent survey article with some comparisons. Quoting 55 "caution should be taken when directly comparing speedup numbers across methods and im- plementations. particular, some of the speedup and efficiency numbers are only theoretical in nature, and many of the parallel time methods do not address the storage communication overhead of the parallel time integrator". [1 is one of very few recent attempts to systematically compare different PinT classes. However, it is limited exclusively to the Dahlquist problem. Thus, it has become conventional to compare new techniques to the existing state-of-the-art methods within the same group of solvers. This is why, in this work, compare RandNet-PPrareal with the original Parareal and its recently improved versions, GParareal [57], and nnGParareal (21].
The rest of the paper is organized as follows. In Section 2 we describe the Parareal algorithm. Section briefly explains GParareal and nnGParareal, focusing on the latter. RandNet-Parrreal is introduced in Section E. while Sections and 6 present our numerical results, and a final discus- sion. A computational complexity analysis of RandNet-Parareal, a robustness evaluation of the proposed algorithm, complementary simulation studies, and other additional results are available in the Supplementary Material.
We denote by v E IR column vector with entries i E f1....,7,, and by |||||| and ||ollo its Euclidean and infinity norms, respectively. We use A 1 xm denote a real-valued n X m matrix, n, m E N, with elements Aij, jth column A(. j E {l.... m.. and ith row A(.-). i E {l.....n) .... We write AT, AI, and || for the A matrix transpose, Penrose pseudoinverse, and Frobenius norm, respectively. I, denotes the identity matrix of dimension n.
Notation.
2 The Parareal algorithm
The idea of Parareal is to solve the d-dimensional ODE (and similarly PDE) system 1 in parallel- in-time fashion, dividing the original IVP into N sub-IVPs
du, h (u, (t .)), te [t,, u, (ta) U.. for i == O..... N dt
3

---

where the number of time intervals N is also the number of available machines/cores//rroesssss u, (t U,) is the solution at time of the 7th IVP with initial condition u(t,) U, E i O,..., N - 1. If the initial conditions were known and satisfied the continuity conditions U; = U, (1) (for the coherent temporal evolution of the system across sub-intervals), then the sub-IVPs could be trivially solved in parallel on a dedicated machine. Unfortunately, this is not the case, as only the first initial condition Uo = u E 12 at time to appears available. To account for this, Parareal introduces another numerical integrator G, much faster but less accurate than $, to approximate the missing initial conditions U,, i = I,.... .,, 1, sequentially, trades off accuracy for computational feasibility, usually taking seconds/minutes instead of hours/days of .7
The algorithm works as follows. We use U to denote the Parareal approximation of (t, U
at iteration &o 2 0. At k = O, the initial conditions {U;) N-1 j= are initialized using a sequential application of the coarse solver g, obtaining U = g
I,..., N - 1, with U Uo- At k 2 1, the obtained initial conditions - are "propagated" through I in parallel on N cores to obtain $ ((*)))) = = -. ,, .N. Note that for every initial condition U we compute both & $(* i.e. a precise evaluation of u, (4: (4 and & (UŞ=T)) inaccurate evaluation of the same term. Hence, we can interpret & and & as functions mapping an initial condition to the next one, thereby evolving ( by one interval. We can then use their difference, (# - (1) to correct the inaccuracy of & on future evaluations. This gives rise the original Parareal predictor-correcnt rule U 9 (U
( - ()) with .....,,.. N 1, k 21 [88] where the sequential prediction 9 (U is corrected by adding the discrepancy I - I computed at the previous iteration k 1. However, this formulation can be changed to use data from the current iteration & [57], and generalized to account for different ways of computing the discrepancy, leading to (21)
U = ((* + f(U)
(2)
where f: _2 + Fd specifies how the correction function & - 16 is computed or approximated based on some observation U E Rd. Parareal uses
Frara(U)_ = (# - (()))
while other variants will be introduced in the subsequent sections. The Parareal solution (2) is considered converged for a given threshold > 0 and up to time 4. S IN. if solutions across consecutive iterations have stabilized. That is, for some pre-defined accuracy level > holds that
| - U. 1la <, o<is SISSLNN-L
(4)
Other stopping criteria are also possible 166 67]. Converged Parareal approximations U is L, are 21
57 58). Instead, unconverged solution values U i > L, are updated during future iterations by first running. & in parallel and then using the cnddiction-correction rule (2). The Parareal algorithm stops at some iteration KPara N when all initial conditions have converged, that is when (4) is satisfied with L = N - 1 and thus KPara = k. Note that during every Parareal iteration > 1, the "leftmost" fine solver evaluation & (U)) ))) is either run from the outcome of previous fine computation U & $1tZ)) or from a converged initial condition J - < This guarantees that, either way, the maximum number of iterations to convergence for any Parareal- -based algorithm is K Para = N, in which case it sequentially attains the fine solver solution, with the added computational cost of running I and evaluating J N times. A Parareal pseudocode is presented in Algorithm in Supplementary Material A
no longer iterated to avoid unnecessary overhead [12 20
3 GParareal and Nearest Neighbors GParareal
The performance of Parareal can be improved by a careful selection of f in ( combined with a better use of the available information present at iteration k. Let DD denote the dataset consisting of Nk: pairs of inputs U U IRd and their corresponding outputs (5 () E Rd, I,...,N, N,
j = =0, O,..., k- - 1, that is
Dj =
= N, j = o...., 1}.
(5)
F. and of can be two different solvers the same solver with different time steps.

---

While Parareal relies on one observation to construct the correction f in 3 GParareal and following works, including this one, use all the discrepancy terms I V and information in Dr to make their predictions. The idea of GParareal is to learn the map R2 + R. -RR.U + (5 () 1): via d independent scalar GPs Rd R. - - GPaaa (U% s= ,....,,d, d, one per ODE dimension, ................................ Grara (U) e Rd. and finally plugged into the predictor-correcn rule 2 In particular, each GP prediction (5)
(Ts
whose predictions are concatenated into faram (U-_)) = (F Grara (U, 1)........
Ja TUF is obtained as the GP posterior mean PDD (U;-)) E R. computed by conditioning
the corresponding GP prior on the dataset PR, i.c. fCPara (U*-)) = p ( We refer to Supplementary Material B and [57] for a thorough description of the algorithm, including all relevant quantities of interest, namely the d GP priors, the likelihood, the hyperparameters and their optimization procedure, and an explicit expression of the posterior means. Here, it is worth highlighting that the GPs are trained once per iteration to leverage the new incoming data, and then their predictions are used to sequentially update the initial conditions in ( Using all information stored in De instead of a single observation (as for Parareal) is the primary driver of faster convergence rates experienced by GParareal. Other benefits of this algorithm are increased stability to different initial conditions, the ability to incorporate legacy data (that is, the possibility of using datasets coming from previous runs of the algorithm with different starting conditions or settings, leading to faster convergence), lower sensitivity to poor choices of the coarse solver I, and the possibility of parallelizing the training of the d GPs over the N available cores. The main drawback of GParareal is the heavy computational burden incurred when inverting the GP covariance matrices, which is of order O (d((N))²) at iteration k. This negatively impacts the algorithm's 's wallclock time, which may be higher than Parareal despite a lower number of iterations needed to converge. This is why GParareal has been proposed mainly for low-dimensional ODE systems with a relatively small number of processors/intervals N (up to hundreds), limiting its use and parallel scalability (57].
The nnGParareal algorithm 21 has been proposed to tackle GParareal's scalability issue, sensibly reducing the computational time and memory footprint of GPs by using their nns version (nnGPs). In this framework, at iteration k, the d GPs are all trained on smaller dataset of size m, D,- composed out of the m nns (in Euclidean distance) of U - in DE. leading to the nnGParareal P(a)
correction FauGPara (U = (f (U ..., Panaaaa (UT with (8)
(U%_1) = 1. (U+_1), 8= 1,,,..,, d.
(5) E 5R,5= 1, d, denotes the nnGP posterior mean computed by conditioning the corresponding GP prior the reduced dataset 1,- of size m. Due to the decreased sample size, each nnGP covariance matrix can be inverted at a cost of O(m²) independent of k or N. However, contrary to GParareal which trains the GPs once per iteration, the nnGPs are re-trained every time a new prediction FanGPara ( (U+,) is made, which are at most N k at iteration ls (as at least k intervals have converged at iteration k), yielding a combined O(d(N - -k)m²) complexity. Several experiments on different ODE and PDE systems have shown that m [15,..., 20} offer accuracy comparable to the full GP 21 at a much lower cost. Although faster than GParareal, nnGParareal still exhibits some of the drawbacks inherited from the GP framework, such the cost of optimizing the hyperparameters through a numerical maximization of a non-convex likelihood, and the use of d scalar nnGPs The latter is particularly critical. On the one hand, despite the possibility of training the d nnGPs in parallel, the inversion of m X m matrix is so efficient that parallel overheads may outweigh the theoretical benefits. On the other hand, when solving PDEs, nnGParareal will incur additional costs due to insufficient hardware resources, as usually d >> N, forcing the d nnGPs to queue among the N available processors, which why the algorithm has been proposed for high- dimensional ODE and PDE systems with d < N. We refer to Supplementary Material B and [21] for more details on nnGParareal, and to Algorithm in Supplementary Material A for the pseudocode of the nnGP training. In the next section, we address the nnGParareal issues by introducing RandNets.
Here, 12 1.
4 Random neural networks Parareal (RandNets-Parareeal
In RandNet-Pararent propose to learn the map Rd U (F via RandNets, obtaining the Panarel-Parareal correction /RandNet.Pan Para. which we then use within the predictor- corrector rule (2). Prior to that, we define how RandNets work in general setting with input U Rd
5

---

and output or target Y E R.. Later in the text we will go back to the input of interest U. Let M denote the number of hidden neurons, and HW (U) be single-hidden-laye feed-forward neural network used to learn & 3, given by
H H$ = W + +) UER.,
(6)
where A E RR- is the matrix of random, non-trainable weights of the hidden layer, E RM a random non-trainable bias vector, and W E MR is the matrix of trainable output weights. Here, RMd + RM denotes an activation function obtained as the componentwise application of non-linear map a IR + R which we choose to be ReLU ozx) = max(z, 0) with æ E R, to satisfy the assumption of Proposition I below. The entries of A and Ç are randomly sampled from given distributions A and Pc. respectively, and kept fixed. After observing the dataset Ds. the output weights W are obtained as the minimum l norm least squares (or simply min-norm least squares) estimator or as the solution of the following penalized empirical minimization problem:
W Dr
lim arg min X-+0
Σ (U,Y)
-Y
2
WeRMxx
which is also called "ridgeless" (interpolation) estimator 30 and can be more compactly written as
pPP lim (x X + ALM) Y 2
(7)
Here, X E N is matrix with (Xa.))) = ((A(Ua.)) + C). 1 ,,..., NK, and U,Y Y. Nixd are the collection of inputs and outputs of Dr in matrix resp.eespy...velly de- fined as (Ua..))
jN + i + 1, O, .,...,N. j. 0, 1. Whenever NI: 2 M and the rank of X E RMM is M, 7 reduces to the standard least squares estimator W Dc (x X while if the rank of X X is Nk, the solution admits closed form
U. (Ya.)) =Y,,
WP (x x) X
We get inspired by 21, where only m nns are used in the training. In this setting, M >>>>> Nk = m, and in this overparametrized linear regression case, the ridgeless estimator interpolates the training data, which is a desirable feature since the problem is genuinely deterministic (29 49].
Several ingredients control the performance of RandNets, such as the dimension of the network M and the choice of distributions PA and Pc. In this work, we take the rows of the weight matrix A and the bias entries of Ç be independent and uniformly distributed. For this case, the approximation bounds are available [25] Proposition 3]. which we report below using our notation.
+ - (*) be an unknown function we wish to approximate with H.y defined in (6). Suppose H* can be represented as (U) = =sge e:lwUU g( ddw for some complex-valued function g on I& and all U E Rd with ||U|| S Q. where (...) is the inner product on Rdd Assume that fe+ max (1, [|] 22222 Ig(w))? 1 dw < 00. For > O, suppose the rows of A are i..dd random variables with uniform distribution on Bp d the Euclidean ball of radius p around O, and that the M components
Proposition (Approximation bound, (25), Proposition 3). Let H* :R + R.
Ç are i.i.d. uniform random variables on [- maam(a..p Qp, 1), max (Qp, 1)). Assume that A and Ç are independent and let σ: : RR- + R be given by oxx) max(ι, 0). Then, there exist a R M xd valued random variable W and explicit (see (33) in (25)) constant C* > 0 such that
E
H
S M. S (U) H' -U) W
A.
and for any 8 E (0, 1), the random neural network A. satisfies
IP( (( -U 1/2
2 1 1-6.
Our choice of PA and Pe satisfies the conditions of Proposition i if |U] < Q. If this is not met, we rescale the ODE/PDE system via a change of variables. We found these bounds empirically useful in informing good choice for the sampling distribution, which we follow. If no prior

---

information were available, the common approach would have been to take PA ~ Unif(- a, a) a) Pc Unif(-b, b) and optimize a, b E RC via expensive voos-validation procedure.
Unlike nnGParareal, GParareal, and the corresponding nnGPs and GPs, training RandNets is so fast that parallelization across the d dimensions is unnecessary. Hence, the predictions of the random network are computed jointly on all d coordinates, yielding the RandNet-Parann correction function
PaaddNet-Paral (U -U R..
(U²_1)
(8)
Here, the estimated weights WP....K are obtained using the reduced dataset 1,-1.x consisting of the WRaddNNr nns of - requiring the retraining of the RandNet for every prediction. Employing a multi-output model instead of independently training d scalar-output models addresses one of the pitfalls of GPs, allowing for better scalability when d > N. The fact that training the RandNets reduces closed-form ridgeless interpolation solution presents a substantial difference and im- provement with respect to (nn)GPs. Moreover, expensive hyperparameter optimization is avoided in RandNets, addressing the other major pitfall of GParareal and anGParareal The pseudocode for training RandNets is reported in Algorithm 3 in Supplementary Material A
In Supplementary Material C we derive the theoretical computational costs of nnGParareal and Paraareal, illustrating them as a function of dimension d and number of processors N in Figure 3 These theoretical findings confirm the significantly superior scalability of RandNet-Parrreal which we observe in the numerical experiments reported in Section 3
In Supplementary Material D we study the robustness of RandNet-Pararee to changes in the number of nns RandNet (and thus the input data size), the number of neurons M, and the randomly sampled network weights A, Ç. Intuitively, one might anticipate that a larger data sample would yield more accurate approximation of the correction 50 - 3, and that a higher number of neurons M would reduce the prediction error of RandNets (as in Proposition TT. One may also suspect the algorithm to be sensitive to the particular sampling seed. Remarkably, our empirical findings demonstrate that these factors have limited impact on the number of iterations needed by RandNet-Parareal to converge, which remains largely consistent (up to few iterations) across different values and ODE/PDE systems, for sensible choices of MeandNet and M. For the end user, this eliminates the need of ad-hoe tuning, making the proposed RandNet-Pareed Paraaeal a convenient out-of-the-hoo algorithm.
Numerical Experiments
In this section, we first compare the performance of Parareal, nnGParareal, and RandNet-Pararees on the viscous Burgers' equation (one spatial dimension and one variable, also considered in nnG- Parareal (2I)), to showcase Parareal and nnGParareal challenges as the number of space discretization and, correspondingly, the dimensions d, increases. Then, we consider the larger system defined on a two-dimensional spatial domain with two non-linearly coupled variables, and the SWES (two spatial dimensions and three variables), representing a suitable framework for modeling free-surface flow problems on two-dimensional domain. Two additional challenging systems, the 2D and 3D Brusselator PDEs, known for their complex behavior, including oscillations, spatial patterns, and chaos, are considered in Supplementary Material [E] The simulation setups used for obtaining the results in this section are provided in Supplementary Material G with the corresponding accuracies and runtimes for RandNet-Pararent Parareal, and nnGParareal reported in Supplementary Material FF
equation,
Let Ti and Tig be the time it takes to run F and I over one interval [t,, f441), respectively, and let N& and Ny denote the number of steps for the fine and coarse solvers over one interval, respectively. We can measure the parallel efficiency of an algorithm via its parallel speed-up Salg. defined as the ratio of the serial over the parallel runtime, i.e. Salg := NT. /Talg Salg captures the wallclock gains of parallel procedures and, unlike other quantities (such as the number of algorithm iterations needed to converge), also includes the model training cost.
5.1 Viscous Burgers' equation
Our initial example is a non-linear, one-dimensiona PDE (illustrated in Figure of Supplementary Material H exhibiting hyperbolic behavior (68), described by the equation
Uz VUax
VCai (x, E [-L, L] x [to, tw),
(9)
7

---

Table Empirical scalability and speed-up analysis for viscous Burgers' equation
d = 128, N
128
Algorithm
K
NTg
Ts
Tmodel
Talg
Salg
Fine
13h Sm
Parareal
90
Os
6m
Os
8h 54m
1.47
nnGPararcal
14
Os
6m
12m
Ih 39m
7.90
RandNctt-Paaaral
Os
Ih 2m
12.61
d = 1128, N
128
Algorithm
K
NTw
Ts
T Tmodel
Talg
Salg
Fine
18h 52m
Parareal
Os
9m
12h 57m
1.41
unGParareal
2s
9m
1h 25m
2h 17m
8.26
RandNet- Parareal
Randreet-Pararaa
4
28
9m
1s
38m
29.98
Speed-up Salg of Parareal, nnGParareal (MmnGF=18), and RandNee-PParares (MrandNet=4, M=100) for the ID viscous Burgers' equation. Ti and Ty are the interval runtimes of the fine and coarse solvers, respectively, K the number of iterations converge. moodel the overall time to evaluate across K iterations, including training and predicting, and Talg thealgorithm runtime.
with initial condition v(x,tg) = vo(x), x E I-L, 0, and Dirichlet boundary conditions 2( t) v(L. t), 2z () (L. t [u. ty]. We use the same setting and parameter values in (21). More specifically, we choose 1, diffusion coefficient
0.01, and discretize the spatial domain using finite difference [I5) and equally spaced points æi+1 zj + Ax, with Ax = 2L/d and 0......d
hence reformulate the PDE as dd-dimensional ODE system.
In our first numerical experiment, we choose N = 128, vo(x) = 0.5(cos( + 1), and ty = 5.9 as in [2D] and consider &/ = IRK, $0 RK8, Ng and N& = 4et, where RK1 stands for Runge- Kutta of order 1, and similarly for RK4 and RK8. The results, reported the top of [able show how RandNet-Parte Parareal converges in fewer iterations and has higher speed-up than Parareal and nnGParareal. The difference the model training costs striking, with the nnGP's being approximately 700 times higher than that of RandNets, reducing thus its potential speed-up.
1),4,, lo = O,
As real-world (ine-dimensional) problems would require higher spatial discretization, increase by one thousand to d = 1 1128, keeping N fixed. Unlike assuming matching hardware resources the system size (as implicitly done in (21], where d N). we deliberately do not increase to assess the algorithms' performances under constrained conditions. Instead, both time discretization numbers are increased to account for the finer spatial mesh 43 As observed from the bottom of Table I as d/N > 1, nnGParareal's issues become more pronounced, as the d scalar GPs cannot be all in parallel across the N processors, but need d/N = 10 runs instead, slowing down the algorithm. In contrast, RandNet- Paarrrell has a training cost comparable with the previous example, leading to an even higher speed-up, running in approximately 38 minutes compared to the almost 13 hours of Parareal.
N₂ 6 and Ng 293 (resulting thus longer Ty and Tig times)
5.2 Diffusion-RRectionn system
We now turn more challenging case study. The
equation 775 (illustrated
in Figure 8 in Supplementary Material H is a system of two non-linearly coupled variables, the activator u(l, y) and the inhibitor v(t, x, v). defined two-dimensional spatial domain
as
dzu DuO Dudgyu Ru: 0 D + Do R.
Here, D., D, are the diffusion coefficients for the activator and inhibitor, respectively, and R. Ry [u, v). Rs = R, (u, v) are their reaction functions defined by the Saizhugh-Nagune equation 42
R2 (u, v)
3
U,
RS (u, u)
where = -3 3 D,, = -e-3, and D, -3 We take (I. (-1, 1)² and E [0, 20). The initial condition u(0, x, y) is generated standard Gaussian noise. We apply no- no-flow Neumann boundary

---

condition D. 0, domain is discretized by the finite volume method [ST], resulting a N,, ddimensional dimmeiiiidimni ODE with N₂ and Ny the number of space discretizations along and y/ respectively. The time integration conducted with RK of variable order for & and y (see Table Supplementary Material G
=
for (z.v) E (-1,1)². The spatial
As in the previous example, we conduct two experiments for this system, with speed-ups and runtimes reported Figure [ In the first one, we increased d and N proportionately (with d/N E [11, 13) while maintaining all other quantities (i.e. 3. MnnGPP RandNe) fixed until N = 256. This scenario reflects situation where more resources are allocated solve larger problem sizes. In contrast, in the second experiment, N remains fixed at 512, with d increasing proportionately with Ng maintain algorithm stability. Moreover, F chosen be RK8, with Ns automatically selected by the used Python library scipy [77] This second setting simulates a scenario with constrained resources, where the user aims solve the system using finer spatial mesh. Table Supplementary Material shows that for N 2 256 and d/N > 1, PaGarrre fails to converge within 48-hour budget. Parareal converges always, albeit a considerably slower rate than LandNet-Pararees which is x3-5 faster than Parareal (and up x120 than the fine solver).
Number processors
Number of processors
256
128
256
120
X
Parareal
3.0
Parareal
100
+
RandNet- Parareal
nnGParareal
Runtimc (hours, loguo)
2.5
RandNet Parareal
nnGParareal
Fine solver
2.0
Fine solver
60
X
40
X
X
0.5
X
20
X
0.0
4
01
7
-0.
A
A
3.5
32
4.7
logu (d)
log (d)
Figure 1: Speed-ups (left) and runtimes (right) of Parareal, nnGParareal (MmnGP=20), and RandNet- system versus the
Parareal (MRandNe:=4, M=100) for the -io-dimeniional number d of dimensions (bottom x-axis) and N cores (top x-axis) capped at 512 to simulate limited resources.
5.3 Shallow water equation
Finally, we focus SWEs on two-dimensional domain, described by system of hyperbolic PDFs
0gh + Vhu = 0, öghu + V(u2h + 1. 122) - g.h
where u (u, represents the velocities in the horizontal
u(t, x. v) and vertical v = v(t, x, v) directions, h = (t,, .r, y) denotes the water depth, b b= b(x, 3) describes a (given) spatially varying bathymetry, and hu can be interpreted as the directional momentum components. The parameter gr describes the gravitational acceleration, while 0, denotes the partial derivative with respect to time, and V the gradient of function f. Following (75). we solve a radial dam break scenario where Gaussian-shapet water column (blue) inundates nearby plains (green) within a rectangular box subject to Neumann boundary conditions, causing the water to rebound off the sides of the box, as depicted in Figure 2 More details the simulation setup are given in Supplementary Material G.I
In this case, our algorithm also converges much faster than Parareal, with speed gain of x 1..3-3.6, while nnGParareal fails to converge within the 48-hour time budget as d >>>>>> N. Although the speed gain is lower than for the Diffusion-Reacet takes to 4-10 hours and 37 days less than the Parareal and sequential solver, respectively.
the improvements are remarkable. GandNeet-Paaaree
6 Discussion and limitations
This study improves the scalability properties, convergence rates, and parallel performance of Parareal and more recently proposed PinT solver for ODEs and PDEs, nnGParareal [21]. By replacing the

---

0.6
2-
0.0
Time 2.55
Time 4.26
Time 8.51
Time 12.77
Time 20.00
0.100
0.075
0.050
0.025
0.000
Figure 2: Numerical solution of the SWE for (a, v) E [-5,, 5] x [0,5] with N 264 and N₂ 133
for a range system times t. Only the water depth h (blue) is plotted.
Table 2: Sp analysis for the shallow water PDE as a d-dimensional ODE system,
235
d
Kpaza
KRandNec-Para
T>
Traa
TRandNet-Para
Spara
SkandNet
15453
52
14
22h 54m
5h 8m
Ih 25m
4.47
16.16
31104
50
13
3d 2h
15h 43m
4h 9m
4.68
17.69
60903
13d 15h
19h 30m
12h 34m
16.73
25.92
105336
8
38d 4h
1d 7h
23h 34m
29.37
38.90
K. is the number iterations converge, T. wallelock time and S. speed-up for the Parareal (Para) and RandNet-Pararees (MARandNet= =4, M=100). T is the sequential runtime of F. The results for nnGParareal (##)n)==2)) are not reported as fails to converge within a 48-hour time budget.
nnGP with random networks, we decreased the model costs (in learning the discrepancy between the fine and coarse solvers) by several orders of magnitude. The reasons behind this are multi-fold. Training of RandNets cheap due to the availability the closed-form solution for its output (readout) weights, and avoids any expensive hyperparameter optimization. Moreover, it is possible to simultaneously learn and predict the d-dimensional correction map instead of d scalar maps (in parallel if the number of processors N is comparable to d. or queuing if smaller). The latter "liberates" RandNet-Paas PPrrrrrl from requiring key/notable improvement with respect to anGParareal. We tested the proposed algorithm systems of real-world significance, such as the Diffusion- solving them fine spatial mesh of up 10 discretization points. These systems and requirements align with those outlined the benchmark PDE dataset (75 as necessary prerequisites for using such algorithms practical scenarios. The strength of GandNett-Paareaa which can be embedded within Parareal with virtually no overhead, irrespective of the implementation or solvers, leading to notable speed gains over Parareal (x8.6-21.2 for viscous Burgers', x3-5 for x X1.3-3.6 6.6 for SWE, and x3.4-4.4 for Brusselator). Moreover, training RandNets is easily conducted with established linear algebra routines, and requires no ad-hoc parameter tuning.
N, extending its application to high-dimensional settings, a
equation, the SWE, and the Brusselator.
the cheap cost of RandNets,
Despite its excellent performance, RandNet- HandNet-Pararcal has limitations common to all Pararcal algo- rithms, its rate of convergence relies the accuracy of the coarse solver G. Although neural net- works can help mitigate the impact of suboptimal choices of & (as observed for GPs in (ararrarr,,),), if the solver is mismatched for the system for example, an unstable solver for stiff ODE RandNet-Parareal, similar to Parareal and (arrrrra,,,,,,, is likely to exhibit non-convergent behavior. would then be of interest investigate MandNet-Paaarea performance when using customized solvers tailored to specific systems, such as those outlined in Section for the shallow water equation and the viscous Burgers' equation, which defer future research.
Acknowledgments and Disclosure F Funding
GG is funded by the Warwick Centre of Doctoral Training in Mathematics and Statisties. GG thanks the hospitality of the University of St. Gallen where part of the results in this paper were obtained.
10

---

References
[1] N. Abel, J. Chaudhry, R. D. Falgout, and J. Schroder. Mulligriidredductionnn--iii-n for the rotating shallow water equations. Technical report, Lawrence Livermore National Lab (LLNL), Livermore, CA (United States), 2020.
[2] G. Ariel, S. J. Kim, and R. Tsai. Parareal multiscale methods for highly oscillatory dynamical systems. SIAM Journal on Scientific Computing, 3866):33340-3354, 2016.
[3] G. Bal. On the convergence and the stability of the parareal algorithm to solve partial differential equations. Domain Decomposition Methods Science and Engineering, pages 425-433 Springer, 2005.
[4] G. Bal and Y. Maday, A "parareal" time discretization for non-linear PDE's with application to the pricing of an american put. In Recent Developments in Domain Decomposition Meth- ods, volume 23 of Lecture Notes in Computational Science and Engineering, pages 289-202. Springer, 2002.
[5] W. Cao, X. Wang, Z. Ming, and Gao. A review on neural networks with random weights. Neurocomputing, 275:278-287, 2018.
[6] L. Carratino, A. Rudi, and L.. Rosasco. Learning with sgd and random features. Advances in neural information processing systems, 31, 2018.
[7] F. Chen, J. S. Hesthaven, and X. Zhu. On the use of reduced basis methods to accelerate and stabilize the parareal method. Reduced Order Methods for modeling and computational reduction, pages 187-214, 2014.
[8] V. Csomor.
Pararealmi, 2023. URL hawpargg//ppp.ooo/rpree///rrrrre
[9] X. Dai and Y. Maday. Stable parareal in time method for first-and second-order hyperbolic systems. SIAM Journal on Scientific Computing, (5(1):A52-A78, 2013.
[10] S. Dong and Z. Li. Local extreme learning machines and domain decomposition for solving linear and nonlinear partial differential equations. Computer Methods in Applied Mechanics and Engineering, 387:1 129, 2021.
[11] V. Dwivedi and B. Srinivasan. Physics informed extreme learning machine (PIELM) rapid method for the numerical solution of partial differential equations. Neurocomputing. 391: 96-118, 2020.
[12] W. R. Elwasif, S. Foley, D. E. Bernholdt, L. A. Berry, D. Samaddar, D. E. Newman, and R. Sanchez. A dependency-driven formulation of parareal: parallel-in-time solution of PDEs as a many-task application. In Proceedings of the 2011 ACM international workshop on Many task computing grids and supercomputers, pages 15-24, 2011.
[13] M. Emmett and M. Minion. Toward an efficient parallel time method for partial differential equations. Communications in Applied Mathematics and Computational Science, 7(1): 101):005-332 2012.
[14] R. D. Falgout, S. Friedhoff, T. V. Kolev, S. P. MacLachlan, and J. B. Schroder. Parallel time integration with multigrid. SIAM Journal on Scientific Computing, 36(6):C635-C66) 2014.
[15] B. Fornberg. Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation, 51 (184):699-706, 1988.
[16] S. Friedhoff, R. D. Falgout, T. V. Kolev, S. MacLachlan, and J. B. Schroder, multigrid-in-timn algorithm for solving evolution equations in parallel. University of North Texas Libraries, UNT Digital Library, 12 2012.
[17] M. J. Gander. 50 years of time parallel time integration. In T. Carraro, M. Geiger, S. Körkel, and R. Rannacher, editors, Multiple Shooting and Time Domain Decomposition Methods, pages 69-113, Cham, 2015. Springer International Publishing. ISBN 978-3-319-23321-55
11

---

[18] M. J. Gander and S. Vandewalle. Analysis of the parareal time- iieeparallel time-integration method. SIAM Journal on Scientific Computing, 29(2):556-577 578, 2007.
[19] M. J. Gander, T. Lunet, D. Ruprecht, and R. Speck. A unified analysis framework for iterative parallel-in -intteee----mee algorithms. SIAM Journal on Scientific Computing. 45(5) 2023, doi: 10.1 00.1137/22M488716.
[20] Garrido, B. Lee. G. Fladmark, and M. Espedal. Convergent iterative schemes for time parallelization. Mathematics of Computation, 75(255):1403-4428 2006.
[21] G. Gattiglio, L. Grigoryeva, and M. Tamborrino. Nearest neighbors GParareal: Improving scalability of Gaussian processes for parallel-in-time solvers. arXiv:2405. 12182v1 2024.
[22] L. Gonon. Random feature neural networks learn Black- type PDEs without curse of dimensionality. Journal of Machine Learning Research, 24:1-51, 2023.
[23] L. Gonon and A. Jacquier. Universal approximation theorem and error bounds for quantum neural networks and quantum reservoirs. urXiv:2307.12904 2023.
[24] L. Gonon and J.-P. Ortega. Reservoir computing universality with stochastic inputs. IEEE Trunsactions on Neural Networks and Learning Systems, 31(1):100-112, 2020.
[25] L. Gonon, L. Grigoryeva, and J.-P. Ortega. Approximation bounds for random neural networks and reservoir systems. The Annals of Applied Probability, 33(1):28 69, 2023.
[26] L. Gonon, L. Grigoryeva, and J.-P. Ortega. Infiniteddimensionall reservoir computing, Neural Networks, 179:106486, 2024.
[27] L. Grigoryeva and J.-P. Ortega. Universal discrete-time reservoir computers with stochastic inputs and linear readouts using non-homogeneous state-affine systems. Journal of Machine Learning Research, 19(24): 2018.
[28] L. Grigoryeva and J.-P. Ortega. Differentiable reservoir computing. Journal of Machine Learning Research, 20(179):1-2,, 2019.
[29] Q. Han and X. Xu. The distribution of ridgeless least squares interpolators. arXiv:2307.02044, 2023.
[30] T. Hastie, A. Montanari, S. Rosset, and R. J. Tibshirani. Surprises in high-dimensional ridgeless least squares interpolation. The Annals of Statistics, 50(2):949 986, 2022.
[31] T. Haut and B. Wingate. An asymptotic pame method for highly oscillatory PDEs. SIAM Journal on Scientific Computing, 36(2):A693 2014.
[32] Y.-L. He, X.-Z. Wang, and J. Z. Huang. Fuzzy nonlinear regression analysis using a random weight network. Information Sciences, 364:222-240, 2016.
[33] C. Herrera, F. Krach, P. Ruyssen, and J. Teichmann. Optimal stopping via randomized neural networks. Frontiers of Mathematical Finance, 3(1):31-77, 2024.
[34] G.-B. Huang. An insight into extreme learning machines: random neurons, random features and kernels. Cognifive Computation, 6:376-390, 2014.
[35] G.-B. Huang, Q.-Y. Zhu, and C.-K. Siew. Extreme learning machine: a new learning scheme of feedforward neural networks. In IEEE International Joint Conference on Neural Networks (IEEE Cat. No. 04CH3754/), J, volume 2, pages 985-999) 2004.
[36] G.-B. Huang, L. Chen, and C.-K. Siew. Universal approximation using incremental constructive feedforward networks with random hidden nodes. IEEE Transactions on Neural Networks, 17 (4):879-892 2006.
[37] G.-B. Huang, Q.-Y. 7hu, and C.-K. Siew. Extreme learning machine: theory and applications. Neurocomputing, 70(1-3) (011-3):49-5011 2006.
12

---

[38] A. Jacquier and Z. Zuric. Random neural networks for rough volatility. arXi::2305.01033 2023.
[39] J. Ji, H. Jiang, B. Zhao, and P. Zhai. Crucial data selection based on random weight neural network. IEEE International Conference on Systems, Man, and Cybernetics, pages 1017- 1022, 2015.
[40] B. Jin, Q. Lin, and Z. Zhou. Learning coarse propagators parareal algorithm. arXiv:2311.15320 2023.
[41] Kar and H. Karnick. Random feature maps for dot product kernels. In Proceedings of the 15th International Conference on Artificial Intelligence and Statistics, volume 22 of Proceedings of Machine Learning Research, pages 583-591, 2012.
[42] G. A. Klaasen and W. C. Troy. Stationary wave solutions of system of seaction-diffusion equations derived from the Fitzhugh-Naguma equations, SIAM Journal on Applied Mathematics, 44 4401);96-110, 1984.
[43] R. J. Vequee Finite Difference Methods for Ordinary and Partial Differential Equations: Steady- and Deee-Deenddent Problems. SIAM, 2007.
[44] M.-B. Li, G.-B. Huang, P. Saratchandran, and N. Sundararajan. Fully complex extreme learning machine. Neurocomputing, 68:306-314, 2005.
[45] J.-L. Lions, Y. Maday, and G. Turinici. Résolution d'EDP par schéma temps pararéel. Comptes Rendus de l'Académie des Sciences-Series I-Mathematics, 332(():661-666 2001.
[46] J. Lu, J. Zhao, and F. Cao. Extended feed forward neural networks with random weights for face recognition. Neurocomputing, 136:96- 102, 2014.
[47] A. Lupo, L. Butschek, and S. Massar. Photonic extreme learning machine based on frequency multiplexing. Opt. Express, 29:28257-28276 2021.
[48] R. Martínez-Peña and J.-P. Ortega. Quantum reservoir computing in finite dimensions. Physical Review E Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, 107(3): 035306. 2023.
[49] S. Mei and A. Montanari. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 75:667-766, 2019.
[50] M. Minion. A hybrid parareal spectral deferred corrections method. Communications Applied Mathematics and Computational Science, 5(2):265-301, 2011.
[51] Moukalled, L. Mangani, and M. Darwish. The Finite Volume Method in Computational Fluid Dynamics. Springer Cham, 1st edition, 2016.
[52] A. Neufeld and P. Schmocker. Universal approximation property of random neural networks. a7X1v:2312.08410, 2023.
[53] A. Neufeld, P. Schmocker, and S. Wu. Full error analysis of the random deep splitting method for nonlinear parabolic PDEs and PIDES with infinite activity. 4rXiv:2405.05// 05290,999922 2024.
[54] A. S. Niclsen, G. Brunner, and J. S. Hesthaven. Communication-awee aware adaptive parareal with application nonlinear hyperbolic system of partial differential equations. Journal of Computational Physics, 37 :483-505, 2018.
[55] B. W. Ong and B. Schroder. Applications of time parallelization. Computing and Visualization in Science, 23(1):1-15, 2020.
[56] G. Pages, O. Pironneau, and G. Sall. The parareal algorithm for american options. SIAM Journal Financial Mathematics, 9(3):966-99) 2018.
13

---

[57] K. Pentland, M. Tamborrino, T. J. Sullivan, J. Buchanan, and L. C. Appel. GParareal: time-parallel ODE solver using Gaussian process emulation. Statistics and Computing, 33(1): 23, 2023.
[58] K. Pentland, M. Tamborrino, D. Samaddar, and L. C. Appel. Stochastic parareal: an application of probabilistic methods to time- aame-prallelization. SIAM Journal on Scientific Computing, 45: S82-S102, 2022.
[59] B. Philippi and T. Slawig. The parareal algorithm applied to the FESOM ocean circulation model. arXiv:2208. 07598, 2022.
[60] B. Philippi and T. Slawig, A micro-macro parareal implementation for the occan-circulation model FESOM 2. arXiv:2306. 17269, 2023.
[61] A. Rahimi and B. Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems (NeurIPS), volume 20. Curran Associates, Inc., 2007.
[62] A. Rahimi and B. Recht. Uniform approximation of functions with random bases. In 46th Anmual Allerton Conference on Communication, Control, and Computing, page 555-561. Curran Associates, Inc., 2008
[63] A. Rahimi and B. Recht. Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning. In Advances in Neural Information Processing Systems (NeurIPS), volume 21, page 1313-1320. Curran Associates, Inc., 2009.
[64] M. Reynolds-Barredo, Elwasif. Mechanisms for the convergence of time-parallelized, parareal turbulent plasma simulations. Journal of Computational Physics, 231(2)):7811-788, 2012.
Newman, R. Sánchez, D. Samaddar, L. A. Berry, and W. R.
[65] A. Rudi and L. Rosasco. Generalization properties of learning with random features. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS 17, 2017.
[66] D. Samaddar, D. E. Newman, and R. Sánchez. Parallelization in time of numerical simulations of fully-developed plasma turbulence using the parareal algorithm. Journal of Computational Physics, 229 6555--6573, 2010.
[67] D. Samaddar, D. P. Coster, X. Bonnin, L. A. Berry, W. R. Elwasif, and D. B. Batchelor. Application of the parareal algorithm simulations of ELMs in ITER plasma. Computer Physics Communications, 235:246-257, 2019.
[68] A. Schmitt, M. Schreiber, P. Peixoto, and M. Schäfer. A numerical study of a semi lemi-lagrangian parareal method applied to the viscous burgers equation. Computing and Visualization in Science, 19(1):45-55) 2018.
[69] C. Sheldon, A. Kolchinsky, and F. Caravelli. Computational capacity of Irc, memristive, and hybrid reservoirs. Physical Review E, 106:045310, Oct 2022.
[70] A. Sinha and J. C. Duchi. Learning kernels with random features. In Advances in Neural Information Processing Systems (NeurIPS), volume 29. Curran Associates, Inc., 2016.
[71] B. Song, J.-Y. Wang, and Y.-L. Jiang. Analysis of new krylov subspace enhanced parareal algorithm for sime-periodic problems. Numerical Algorithms, pages 1-22, 2023.
[72] G. A. Staff and E. M. Rønquist. Stability of the parareal algorithm. In Domain Decomposition Methods in Science and Engineering, pages 449 Springer, 2005.
[73] J. G. C. Steinstraesser, P. da Silva Peixoto, and M. Schreiber. Parallel-in-time integration of the shallow water equations on the rotating sphere using parareal and MGRIT. Journal of Computational Physics, 496: 12591, 2024.
[74] Y. Sun, A. Gilbert, and A. Tewari. But how does it work in theory? Linear SVM with random features. In Advances in Neural Information Processing Systems (NeurIPS), volume 31. Curran Associates, Inc., 2018
14

---

[75] M. Takamoto, T. Praditia, R. Leiteritz, D. MacKinlay, F. Alesiani, D. Pflüger, and M. Niepert. Pdebench: An extensive benchmark for scientific machine learning. In Advances in Neural Information Processing Systems (NeurIPS), volume 35, pages 1596-1611, Curran Associates, Inc., 2022.
[76] A. Turing. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society B, 237:37-72, 1952.
[77] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski,
P. Peterson, W. Weckesser, J. Bright, et al. Scipy 1.0: fundamental algorithms for scientific computing in python. Nature methods, 17(3):261-272, 2020.
[78] W. Wan, Z. Zhou, J. Zhao, and F. Cao. novel face recognition method: Using random weight networks and quasi-singular value decomposition. Neurocomputingg 151:1180-1188 2015.
[79] Y. Wang and S. Dong. An extreme learning machine-based method for computational pdes in higher dimensions. Computer Methods Applied Mechanics and Engineering, 418:116578, 2024.
[80] C. K. Williams and C. E. Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
[81] Y. Yang, M. Hou, and J. Luo. A novel improved extreme learning machine algorithm in solving ordinary differential equations by Legendre neural network methods. Advances in Difference Equations, 429, 2018.
[82] L. Zhang and P. N. Suganthan. A survey of randomized algorithms for training neural networks. Information Sciences, 364:146-15: 2016.
15

---

A Pseudocodes
This section provides pseudocodes for the implementation of Parareal (Algorithm D and the training procedure for learning the discrepancy F - via nnGPs in nnGParareal (Algorithm2) 2) and RandNets
Sanddee--Pararell (Algorithm 3)
Algorithm I: Parareal (generic)
Input: Initial condition u at time to, number of intervals N
Ji=1 with K the number of iterations to Output: Converged initial conditions {U; K N-1 convergence
Initialization
Rescale the ODE/PDE system such that each coordinate takes values in [-1, 1]
L+1
U u for +- to N 1 do U 1 (U)-1)
end
for k 1 to N do
Compute $ (UU i 1, I, in parallel for i + L 1+110 N do U (U f(U))))
/* Update the initial conditions */
end
Convergence checks for + + to 1 do if U U 1 al < c then
L 1
/* Update converged interval counter */
else end
break
end
L. == N then
break
/* All intervals have converged +/
end
end
16

---

Algorithm 2: nnGP training procedure within nnGParareal
Input: Input U 1-1: dataset DR2 number of nearest neighbors MnnGP, number of random restarts for loss maximization Rsstart Output: Prediction (un of ( (U)))
Initialization
/* Find the MMnGP nns to U i-1 and compute the reduced dataset (10) D1-1. + ((U)))m)) (l-an) y (l-nn)), I == l,... MGGGP C Di
Ut-1"
/* Both loops can be massively parallelized */
for < to N do
/* Training */
for j to Nsstart do
/* Random restarts to avoid local minima when maximizing 12 */
Sample 0f at random
Maximize 12
numerically using 09 as initial value; obtain 0]
end
Find 0*
such that
log p( (U,) logp(Y(,))U,03) j I,..., Retart
/* Predicting */
Compute u 1.k (U 1 with using
end
Set Jun (U 1 + (WD, 1. (1) PD, (U++1))
Algorithm 3: RandNets training procedure within RandNet-Parareal
Input: Input U dataset DR, number of neurons M, number of nearest neighbors MRandNet
Output: Prediction JRamdNer (U+=1) of (5 - {))())
Initialization
Ensure each ODE/PDE coordinate takes values in [-1, to U and compute the reduced dataset 10 +/ Di-1,K + {(U (I-an) (1-mn)
/* Find the Net
Ut_ U ), 1= I,... MRandNa C Dr:
Sample Aw.j Uniform -1, 1), w = 1,..., ..MMj = .....,
Sample S., Uniform( 1, 1), 1, 1.....M. M
Let X E RmxM
Training
T < a (AUT + +))
Using broadcasting on Ç */
if Trank(X* X1 TTX) X) = M < m then
WD.-1A -X
Least-squares estimator +/
else
WD (X XT end
/* Ridgeless interpolator */
Predicting
JRandNer (UL1)) -WWPiLL.) (AU +-1 { +5)
17

---

B Additional details on the nnGParareal correction function
In this section, we provide more details on the nearest neighbors (nns) Gaussian process modeling. the mathematical expressions of the anGParareal correction function SunCPara and the reduced dataset 1, These are not explicitly presented in the main text as they require additional notation, which we believe does not enrich the explanation. While the deseription of GPs presented here is for nnGPaaareeal (and the corresponding nnGPs), it immediately generalizes to GParareal by replacing the reduced dataset D1-1. with the full dataset Dz. The interested reader can find more details the original papers, [2 and [57].
Let the set of inputs U E IR and outputs (3 - (U) R i= N. = o...., k collected by iteration k, be denoted by Uy, and Vn, respectively, Now, define D2-1,8 as the restriction of Dr to the m nns of U in Us, namely
D D1-1 {(U ((ul) U._ Y Y( 1= I......, m) C DR.
where y-(kua) = (5 (U ) E Vi, and U (b-m)) is the (th nn of U 1 in Dry i.e. the /th ordered Uh
U
statistics the set formed out of Euclidean distances | - U between U and any U U
(((-m)) U4 such that, for any U' Uss U' + Ur 1,...,4
That is, there exists U1, U,...., U1 Uḥ
have
- ... - s |0%_1 U.l s |1 U'I
(10)
Finally, let Ü, Ý Rmdd be the matrices of input nns and outputs collected in D; respectively.
In nnGParareal, following the Bayesian framework, GP prior is placed over the correction function F ② for each of the coordinates as
F- GP Kar). 1,.....,d
where uGP Rdd + R the prior mean function, taken be zero for all 1,..., d, and KGp: : Rd Rd IR the exponential prior variance kernel function
(U, (U, = oー -U1////)''
with of and of denoting the input and output length scales, respectively. Differently from the prior mean, the prior variance is the same across the d components. Then, each nnGParareal prediction
((s) JunGPara ( E R, = 1.....d. d, is obtained from the GP posterior mean (D..I.K (U R. computed on the reduced dataset 1.- given by
Flal InCPara (U 1 I. (U
(a)
:= ((U,
oeg Im) -
(1a);
(11)
where ICÚÜ, U. RT is vector of covariances between every input collected in Ü and U defined as (KCU (((U,UL-))) = (T)) ,U?,),,
.,....,......... m. and K(Ü,ū) E Rmxm the (U(r.)) ( T), r,q = l,..., m. Here, of
covariance matrix, with (AC(,, = Ker((U (U()))) denotes a regularization term, also known as nugget, jitter, or regularization strength, which is added improve the numerical stability when computing the inverse matrix, see (21) for further details. The hyperparameters 0 = (o o oPeg) entering into the posterior mean and prediction (IT) control the performance of the GP, and optimized by numerically maximizing the marginal log-likelihood:
(.,)) log (,,),,,,U... (,, (K, Ú) TT. - log det deC(( Ú)),
(12)
where (C.,,)) depends on 0 through the kernel KGP, and det((A) (A) denotes the determinant of square matrix A. For a thorough treatment Gaussian processes, including derivation of the likelihood and of the posterior distribution (which is Gaussian with mean in I1 see (80).
C Computational complexity analysis
Consider the d-dimensional initial value problem T for some (O/P)DE, Let the number of subintervals (data points) at each kth iteration of the PinT algorithm. For any kth iteration of
18

---

the scheme, a total of Nk data points, each d-dimensional, are available. Here, we provide the computational cost of PaandNet-Paarrall and compare that of mnGParareal, the state-of-the Parareal algorithm proposed in [21]. Both RandNet-Parareal and nnGParareal use only the reduced data set of m nns to a given point to construct its image-prediction via 2 Note that the m nns (in Euclidean distance) to some point U E Rd among Nl: available points are found at cost which is at most linear in the sample size, that is O(mNk) (for moderate dimensions d, one can get an improved cost O(m log(NN) k)), logarithmic in the sample size) [21]. Since our goal is to compare the computational complexities anGParareal and RandNet-Pararge as a function of d, we consider the worst-case complexity of the nns search.
Given an input - E R i = , at iteration k. the computational model cost of a prediction U produced by all d models of ManGP-nnGPS iteration k via the predictor-corrector rule (
with nnGParareal correction (TT) and MmGGP nns is given in 21 as
TanGP (k) CDaCP (Tatart V1)x
(
ManGPd
+ manGrd +
M MBGP + -anGPPd ManGP-Kk
-== C==K(U,U)
D:= DB menn
B F.D
BD.Y
nearest neighhors
=CmaGP.Nki +CuuGp-NK(nreeMtaat Naten N
V 1) MUnGP + d(m + 2mmnGP) + Nk))
with CHnGP being some constant that in general does depend on k, MnnGPP and d. Also, Rrog and Vistart correspond to the number of random restarts and the number of explored values of the regularization penalty in the kernel regression (associated to the hyperparameter optimization (see [21 Section 4.5l). respectively. Furthermore, V is the maximum operator, and the factor (Natart? d/ N V 1) > 1 follows from the fact that d independent nnGPs and hyperparameter optimization are parallelized over the N cores.
In Parareel,, the correction term JRandNet-Paraa is modeled by the random weights neural network and evaluated as (8). Again, only RandNet nns (in Euclidean distance) to U are used to construct the prediction, leading to the following computational model cost at iteration k:
TRandNot (k) S CRand Nk:
MidtRaanNee + M MRRRaanNet + Mr
(AU4()
E==X-XT
Et
MmRanaNee d +
M2d
+ MdmTRandNet
MandNNa Nk)
E1
W
WT.X
nearest neighbors
k M2RndadNee + ddM² + mRRadNNee) + MRandNet Nk)),
where M is the number of hidden neurons, r is the rank of the covariance of activated neurons Σ (mind that the pseudoinverse of Σ would contribute cubically m only if Σ of full rank numerically, which is not observed empirically) and CRandNet is constant independent on N, k, M, d, MRandNet- The factor 1 I/N in the first inequality corresponds to parallelization over N processors.
We note the following differences in costs between these two algorithms according to realistic situations:
d >>>> N in most relevant applications, especially for PDEs. Hence, (rega d/V1) > 1, limiting the benefits from parallelization for nnGParareal. In the considered experiments, we had access to a maximum of approximately N = 500 processors, while we considered up d 10 It is easy see that THaGP is quadratic in dimension d, while TRaaNNet is only linear. This difference is mainly due to the factor (Netart V 1) 2 1 in ThaGP as opposed 1/N in RandNett
Although M > MPnGP. M 100 is sufficient for consistent performance across range of systems, as shown in our numerical experiments.
nnGParareal incurs additional cost due to hyperparameter optimization (21,, necessary for tuning the kernel input and output scales and the regularization strength for each of the d dimensions, which is performed by maximizing the loglikelihood. First, to explore the parameter space and allow for multiple starting points given the nonconvex optimization
19

---

Model computatioal cost
Model computational cost with fine solver 7129
Ra nalneet Pararcal
RandNet-Parareal
l0go hours
ninPPararea
rnGiaraneals
logo IT1
de2392
deli000 Wo532
deilate
de2dasg No512
Ceb M6512
M-256
No512
Mo512
No512
No512
No512
Panel
Panel
Figure 3: Theoretical model cost (panel A) and thcoretical total cost (panel B), as functions of the dimension (and the corresponding N). The results reported in terms of log1o (hours).
problems, both freg and Rsttart should be set large. Second, each loglikelihood maximization conducted per dimension of the system requires large number of iterations performed sequentially. Hence, GDaGP >> CRandNet. with CRoGP depending, in general, on k, MnGGP d, opposed to RandNet- Indeed, RandNet requires no tuning significant advantage with respect to GPs and nnGPs, making it more user-friendly), neither for the distribution of the random weights nor for the regularization parameter due the use of the ridgeless estimator, Empirically, we observed CanGP/ CCRandNet to be up to 1000 (this can be seen in Figure T..
We emphasize that since the ReLu function is chosen activation in RandNet, the matrix X activated neurons is sparse with sparsity degree 7y Hence, the computational complexity RandNet could be further improved, the computational complexity of sparse operations is proportional to the number of nonzero elements in the matrix. We intentionally left these arguments out of the complexity analysis, since we do implementation.
sparse operations in code
The upper bound of NandNet could potentially be improved further, as additional paral- lelization may occur during standard matrix operations, depending the specific computing environment.
Figure 3 illlstrates the theoretical model costs THmGP and T TaadNNet (Panel A) and theoretical total costs obtained by adding the coarse and fine solver costs (Panel B), as functions of the dimension (and the corresponding N). The results are reported in terms of log (hours). To calibrate the constants in both complexity bounds, used the total empirical computational cost in Figure [ together with its breakdown described in Table g Panel shows that RandNet-Pararee Parareal displays significant improvement scalability with respect the state-of-the-art Parareal algorithm nnG-
Parareal, while Panel В demonstrates that whenever the cost of the fine solver is added, our results in full coherence with the empirical results.
D Robustness study
In this section, we study the robustness of Pararett to changes in the number of nns WRandNet. the number of neurons M, and the randomly sampled values of neural network weights A, Ç. Our empirical findings (for two of the three considered PDEs) demonstrate that the iterations KRandNet-Pand Para convergence for FandNet-Paranee remain largely consistent despite variations in these factors. This ensures robust performance across a broad spectrum of parameter values, reducing users' need for extensive tuning. For computational tractability, we limit the robustness analysis to relatively small systems, such as Burgers' equation with d = 128, and the Diffusion-Reacii equation with d 722, conducting 100 weight samplings for each system. For every set of weights, we iterate RandNet-Paane Parareal across MRandNet values ranging from 2 to 20, and M values ranging from 20 to 500 increments of 10. The proportions of iterations needed converge across 100 runs for different values of MRaudNet and M for the Burger's s and diffusion -eaaiionn equations are reported in Figures [n respectively. Although we observe some minor differences between two systems, the main trend is clear: as long as reasonable values of MRRanNe and M are chosen, the iterations to convergence for RandNet-Paran
vary at most by a few units when changing the values of m Rend M particular sampling seed of weights. Nevertheless, larger M might improve
20

---

the performance, since, in this case, RandNets operate in the interpolation regime, as discussed in Section
1.0
1.0
0.8
0.8
Krunches
Proportion
0.6
0.6
Kaurdeistar 10
Krardeestar 12
0.4
0.4
Kiandeerraa Kussnaa
0.
0.2
0.0
14
15
16
20
0.0
50
80
110
100
200
520
260
320
390
500
410
640
470
500
M. number neurens
Figure 4: Histogram of the iterations convergence K RandNet-Parss of LandNee-PPrareal for 128 for Burgers' equation. We sample the network weights A, Ç 100 times. For each set of weights, we run RandNet-Parareal for MandNNet (2,3,.. 3....... 20} and M {20, 30, 40,.., 500]. The left and right panels show the aggregated histograms of K RandNet-Para versus RandNet and M, respectively.
1.0
1.0
KrurdMistee 10
Impportton
Proporttoo
tvo.
Kennden 12
Krundhwshar
KrundMistar Kenthaw
15
1
80
140
170
200
220
260
200
370
350
380
410
440
470
500.
Miandia
M. number neurens
Figure 5: Histogram of the iterations to convergence K RRandNet-Para of RandNet-Parareal for = 722 for Diffusion-Reaction weights, we run RandNet-Pararee for Randne E (2, 3,..., 20} and M E {20, 30, ,0, The left and right panels show the aggregated histograms of K RandNet-Pars respectively.
equation. We sample the network weights A, Ç 100 times. For each set of
RandNet and M,
E
Additional numerical experiments: 2D and 3D Brusselator PDE
Here, we carry out an additional scalability study for the 2 and spatial dimensional Brusselator PDE. This model is two-component reaction system that exhibits complex behavior, including oscillations, spatial patterns, and chaos. is described by
0,u Do 2 (1 b)u+ ou?
and
b bu cu..
In chemistry, the components refer to the concentration of two substances, whereas the constants Du, Dj are the respective diffusivity each component, indicating the rate which the substances spread out in space. Moreover, the parameters and are related to reaction rates. In our experiments, we used D, 0.1, D, Dp. a = and 3. We take [0, 35), (u, v) E (-1, ( for the 2D Brusselator, and (u, v) E 1)3 x (-1, 1)3 for the three spatial dimension case. We initialize the u values at time = by setting them equal to a, and the values by taking them normally distributed over the spatial grid. Further details regarding the number of spatial
21

---

Table 3: Simulation setup for the 2D and 3D Brusselator
Domain
Nu No
d
8f
Sat
f
Fst
N
v, f) (-1,1)2 (-1,1)²
|0,35
32
2048
RK1
0.034
RK4
1e-
512
(u, v, f) (-1,1) (-1,1) [0.
8192
RK1
0.033
RK4
1c-7
512
(u, v, E (-1,1)*
(-1,1)* x [0,35
20
16000
RKI
0.052
RK4
le"7
512
(u, v. () E (-1,1) (-1,1)) [0.35]
25
31250
RK1
0.057
RK4
1e-7
512
N₂ and N, are the number of spatial discretization points for and along each spatial dimension, yielding a d = 2N or d = - dineennsioo ODE, depending on the considered system. & and f denote the coarse and fine solvers, respectively, while the At subscript refers the timestep.
The number of nns used for D.-1.k in nnGParareal and RandNet-Parareat are MaG = 20 and ManndNet 4, respectively, N is the total number of intervals.
discretizations, the number of intervals N and the order of the solvers & and is given in Table 3 Figure 6 hiighlighss the strong scaling advantages of GandNet-Pararees compared to nnGParareal, setting N 512 and restricting the runtime budget to maximum of 48 hours, as done in the other
test cases.
Spatial dim dimensions
Spatial dimensions
250
1
Parareal
2.0
Runtime thours. logi))
PanaNetsParareal
nnGParareal Fine solver
Spieed up
X
X
100
X
X
X
A
0.5
50
-1.0
X
15
4.2
logus (d)
logio (d)
Figure 6: Scalability study for the 2 and 3 spatial dimensional Brusselator PDE. We used N 512 and 48 hours runtime budget. nnGParareal for (d) 3.9 is estimated, the algorithm does converge within the 48 hours runtime budget.
F Accuracy and runtimes across models and algorithms
In Table beloww we report the accuracies and runtimes (shown in parentheses) for GandNet-PPaareet Parareal, and nnGParareal. The accuracy is measured with maximum absolute error (mean across intervals) with respect to the true solution obtained by running &0 sequentially. Interestingly, accuracies are far below the pre-defined accuracy level e, with RandNet-Parareal achieving the lowest one in all but one experiment, with much smaller runtimes across all case studies.
G Simulation setups
This section summarizes the simulation setups used for producing the results discussed in Section in the main text. The tables below report the space and time domain of the considered PDEs, the number of spatial discretization points N₂ (and Ny. in case of -vo-dimensional spatial systems), the numerical solvers used for and G, their corresponding numbers of time steps per interval, the number of intervals N, and the number of nns used for nnGParareal (ManGP) and RandNet- Parareal (aande) particular, Table refers to the viscous Burgers' equation, Table 6 to the Diffusion Reacttio equation, and Table 7 to the shallow water equations (SWEs).
22

---

Table 4: Accuracy and computational cost of the three considered algorithms
PDE
RandNet-Pararcl
Pararcal
nnG Pararcall
Burgers' d 128
1.06e (Ih 2m)
1.85e (8h 54m)
1.32e (Ih 39m)
Catuuioon-Reactonn 7.2c.
3.56e (23m)
1.83e (1h 40m)
5.71e (Ih 11m)
Diffusion-Reactonn = 3.3e
8.56e (33m)
2.45c (7h 52m)
not converged
Cnfrusion-Reattion
2.5e*
8.00c (Ih 57m)
7.43e (9h 50m)
not converged
SWE d = 3.1e*
6.75e (4h 9m)
5.15e - (15h 43m)
not converged
SWE = 6.
8.54e (12h 34m)
2.84e (19h 30m)
converged 3.38c (7h 31m)
Brusselator 2D 2c
2.09e (2m)
3.16c -8 (7m)
Aceuracy and computational cost comparison GeandNet-Parareal Parareal, and nnGParareal for different PDEs, with runtimes reported in parentheses. The accuracy is measured maximum absolute error (mean across intervals) with respect to # run sequentially.
Table 5: Simulation setup for the viscous Burgers' equation
Domain
N₂
d
Ng
F
N
N
MenGP
MRandNet
(x, () [-1, 1] (z, 1
0.9.
128
128
RKI
4
RK8
4e
128
18
5.9
1128
1128
RKI
293
RK8
6e
128
N₂ is the number of space discretizations, the same d here. g and I denote the chosen coarse and fine solvers, with corresponding time discretization steps per interval Ng and Nş, respectively, Here N is the number of intervals, while ManGP and ManddNet the numbers of nns used to create 1 for nnGrrreel and RandNet-Parareal, respectively.
G.1 Simulation setup for the SWEs
Here, we give more details on the radial dam break simulation of Section 5.3
Our domain consists of a rectangular box defined as (x., 3) [-5, 5] [0, 5), which we evolve temporally over E [0, 20]. Following [75]. as an initial condition, we place chhusian--shapd column of water centered at
0.25
0
(z. v) = 2.5, 1.5). with covariance matrix Σ
We use Neumann boundary
0
0.25
conditions, and evolve the system using N = 235 intervals over four increasingly finer spatial meshes, as described in T We used the ParareaML. [8] Python package to implement the SWEs and corresponding numerical solvers.
23

---

Table Simulation setup for the Diffusion-Reaction equation
Domain
N₂
N₂
d
g
N₂
N
N
(z, E [-1, 1] 00,
19
19
722
RKI
RK4
NA
(z, E -1,
00, 20
28
28
1568
RK1
1
RK4
NA
128
(2, E -1,1] [0, 20
3362
RKI
RK4
NA
256
(x.y, E [1,1] x 20
77
77
11858
RK4
1
RK8
NA
512
(2.y, E 1] 20,
113
113
25538
RK4
2
RK8
NA
512
(2,y,t) E -1, 00, 20
164
164
53792
RK4
4
RK8
NA
512
(2, E 11,,1
0, 20
235
235
110450
RK4
8
RK8
NA
512
N₂ and N, are the number of spatial discretization points for æ and y, respectively, yielding d 2N,,,, S..yy-dimensiona ODE. & and & denote the coarse and fine solvers, respectively. The number of nns used for D;-1.k in nnGParareal and RandNet- LandNet-PPrareal are MnaGP = 20 and ManddNet = 3.
respectively. Nig is the time discretization steps of v per interval. Ng NA since H's step size chosen by scipy Runge- Kunge-Kutta method [77].
Table 7: Simulation setup for the SWES
Domain
N,.
N
Ng
f
N
(I.8,f) E
5,
x
5
[0, 20
15453
RK1
7
RK4
1
235
(x, E
5, 5]
X
0, 5
[0, 20)
144
72
31104
RK1
8
RK4
2e
235
(x, E
5,
X
[0, 20
101
60903
RKI
RK4
4e
235
(x. E 5,,
5]
X
0,
5
[0, 20
264
133
105336
RKI
24
RK4
5e
235
N₂ and Ny are the number of spatial discretization points for x and y, respectively, leading to an ODE of dimension d 3.N, 9 and F denote the chosen numerical coarse and fine solvers, respectively, with Ni and N being their corresponding time discretization steps per interval. In all cases, the number of nns used to create D;-1,k MpnGP 20 for nnGParareal, and
MandiNt = for RandNet-Pararce
24

---

Illustration of some PDE solutions
1.00
0.75
0.50
0.50
0.5
0.00
0.25
-0.75
1.00
Figure 7: Numerical solution of viscous Burgers' equation over (xr,
[-1,1]
[0,5.9) with
1128 and initial conditions and additional settings as described in Section 5.1
Time 0.00
Time 0.78
Time 1.56
Time 3.12
Time 5.86
1.0
0.0
-0.5
-1.0
Time 8.79
Time 11.72
Time 14.65
Time 17.58
Time 20.00
1.0
0.0
-0.1
-1.0
-1
-1
-1
-1
Figure 8: Numerical solution of the Diffusion- Diffusion-Reaction equation over (x, v) E [], with N₂- Ng = 235 for range of system times t. Only the activator uit, x, y) is plotted. The initial conditions and additional settings are described in Section 5.2
I Additional simulation results for the Diffusion-Reaction equation
Here, we complement the results of the speed-ups and wallclock times reported in Figure I in the main text, with a detailed breakdown of the number of iterations to convergence, the runtimes of the coarse and fine solvers, the overall cost of training the model (up to convergence), and the total runtime, reported in Table
25

---

Table 8: Speed-up up analysis for the SrffusioonRReaction equation
722,
Algorithm
K
NT
Ts
Tmodel
Talg
Sale
Fine
2h 2m
Parareal
53
Os
2m
Os
Ih 40m
1.22
nnGParareal
16
Os
2m
40m
11m
1.72
RandNet- Pararraa
12
2m
23m
5.36
d 1568, N
128
Algorithm
K
NT8
Tmodel
Talx
Salg
Fine
21m
Parareal
93
Os
2m
3h 3m
1.42
PaGParaeea
20
Os
2m
10h 32m
11h 12m
0.39
RandNet- Parrrraa
12
Os
2m
4s
25m
10.26
3362,
256
Algorithm
K
NTw
Ts
Tmodel
Tagg
Say
Fine
10h 58m
1
Parareal
195
2m
7h 52m
1.40
RandNet-Paare Paraareal
12
3m
33m
19.87
11858, N
512
Algorithm
K
NT&
Ts
Tmodel
Talg
Saig
Fine
2d 16h
Parareal
58
1s
7m
Os 1m
6h 59m
9.23
RandNet-Pararee
6
2s
8m
49m
78.44
25538,
512
Algorithm
K
NTg
Ty
T Tmodel
Talz
Saig
Finc
7d 16h
1
Parareal
27
8s
22m
Os
9h 50m
19.25
RandNet- Parrrraal
5
23m
57m
97.40
d 53792, N
512
Algorithm
K
NT
Ts
Tmodel
Talg
Salg
Fine
21d 7h
Parareal
19
13m
26.60
RandNet-Parareal
4
42s
60m
4m
4h 6m
124.87
110450, N 512
Algorithm
K
NT4
Ts
T Tmodel
Talg
Saig
Fine*
56d 2h
Parareal
14
3m
2h 38m
1s 7m
1d 14h
35,84
RandNet-Pararee
4
3m
2h 37m
10h 48m
124.52
Simulation study on the empirical scalability and speed-up of Parareal, nnGParareal (with MonGP = 20), and RandNet-Parareal (with MRRndNet = 4 and M = 100) for the Diffusion- Recationn equation. T3 and T& refer to the runtimes per interval of the fine and coarse solvers, respectively, while NTg is the runtime of the coarse solver over N intervals. Tmodel corresponds the overall time evaluate f. including training and predicting. until convergence iteration K. Talg the total algorithm runtime, while Salg is the parallel speed-up. "Fine** indicates that the total runtime has been estimated extrapolating data from the other algorithms. Missing anGParareal rows for d 2 3362 are due to convergence failure within a 48-hour time budget.
26

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: All claims, including but not limited to time gains, model training times, numbers of iterations to convergence, speed-ups, and scalability are backed up by empirical results reported in Tables [ 2 and Figure [ in the main text, and Figure and Table 8 in Supplementary Material [|
comparison between theoretical and empirical results is also provided. These claims are stated in the abstract, introduction, and in the final section.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the 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 can be expected to generalize other settings.
It is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: relevant limitations are discussed in Section 6 second paragraph. The robust- ness of the proposed algorithm to several factors, such as the number of neural networks M, the number of nearest neighbors m and the sampled random weights A, introduced at the end of Section 4 and investigated in details in Supplementary Material [ The scaling performance of the algorithm with respect to the number of cores and model dimensions d is extensively discussed in Section 5 Finally, rescaling of the system is proposed if the data do meet the condition |U
Q in Theorem
Guidelines:
The answer NA means that the paper has no 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 a separate "Limitations" section in their paper.
The paper should point 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 few datasets 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 is low 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 of the proposed algorithms and how they seale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
27

---

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 to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Proposition | clearly stated with all the required assumptions. The proof is not given, cite this result from [25], adjusting their notation to match our, clearly mention. The derivation of the computational complexity analysis is provided with all relevant details Supplementary Material
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 to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We have taken care of ensuring the reproducibility of all results through precise use of notation, and by detailing pseudocodes for the algorithms in Supplementary Material A Additionally, we comprehensively describe the simulation setups both in the main text and in Supplementary Material G Moreover, a link to GitHub repository with a step-by-step Jupyter notebook outlining RandNet-Paree Parareal, and the necessary code to reproduced the results has been provided in Section I in the main text.
Guidelines:
The answer NA means that the paper does not include experiments.
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 to 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 is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is 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.
28

---

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 new algorithm, the paper should make it clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.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 to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: 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?
Answer: [Yes]
Justification: All experimental results are fully reproducible, with code provided via GitHub repository, with a link shared in Section the main text. Each simulation and its corresponding analysis are clearly labeled, and step-by-step Jupyter notebook is provided aid the reader in becoming familiar with the API's usage. The repository follows the best practices of the most common ML repositories.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips cc/
Sublic/guidese www.ddee//oeeuumiieeeeeeeeeeeeee
(CodeSSumissiioooooliooo. for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not 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 NeurIPS code and data submission guidelines (https
//nips.cc/publl httpsssccpbbblccuuiiiiiie wwwisiicguieee/eeeeeeeee http://///.iieee//eeuuubiiiii ssi ssionpolicy for more details.
The authors should provide instructions on data access and preparation, including how to 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 much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: 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?
Answer: [Yes]
Justification: All the experimental setups, training details, and simulation parameters are described in the text, mainly Sections and Moreover, they are also summarized in Supplementary Material
29

---

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, or as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
Justification: We acknowledge that error quantification for the speed-up might be of interest in some situations. However, given the runtime of our experiments, this would be too computationally expensive obtain. Nevertheless, we reported two robustness studies for two different, smaller systems among the ones considered (Figures 4 and (5), where the performance the algorithm averaged across multiple runs. There, display the more informative empirical distribution instead of just the error bars.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
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).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report I-sigma error bars, but one should state it. The authors should preferably report 2-sigma error bar than state that they have a 96%
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of 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
For each experiment, does the paper provide sufficient information the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: All results report the execution runtime. Details the hardware used are provided Section I in the main text.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't 't make into the paper).
30

---

9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics /tt.h//....rpp/..oooo.....oooooa httpp/ub.....tti.cccc.e..eeeeee.
Answer: [Yes]
Justification: We have reviewed the Code of Ethics and found particular area of concern regarding our research.
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.. there a special consid- eration due to laws regulations in 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: The potential positive impacts are implicitly mentioned in Section ! and in Section (when referring having chosen systems of real-world significance, with the necessary prerequisites for using the proposed algorithm in practical scenarios). By enabling faster convergence times with minimal overhead, LandNet-Pararees can be applied wide range of applications, such as plasma physics simulation, weather forecasting (both mentioned in the introduction), leading positive societal impacts.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no socictal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g.. disinformation. generating fake profiles, surveillance), fairness considerations (e.g.. deployment 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 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 as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (c.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: We relied on publicly available models, simulating the relevant data as de- scribed in the main text and in Supplementary Material.
31

---

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 to 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 challenging, and many papers do not require this, but encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
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?
Question:
Answer: [Yes]
Justification: All creators have been properly credited both in terms of published scientific papers, and publicly available code and libraries (e.g. for some specific Python libraries).
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
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 (c.g.. website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, phcorswithcode thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license 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 is available online, the authors are encouraged to reach out 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]
The code, simulations and associated analyses are publicly released with permissive licence.
Justification:
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of 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 anonymize your assets (if applicable). You either create an 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 screenshots, if applicable, as well details about compensation (if any)?
32

---

Answer:
[NA]
Justification:
The paper involves neither crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, 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
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 on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: The paper involves neither crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does 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 to 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.
33


Reviews:
[{"review_id": "w79I6TSOC6", "reviewer": "Reviewer_FREq", "summary": "The paper introduces their method RandNet-Parareal, which is a method to solve differential equations. Their method can be categorized as a Parallel-in-time technique that aims at parallelizing solvers in the temporal domain. They extend the Parareal algorithm by using RandNets which learn the difference between a coarse and a fine solver. They show experimentally a speed-up of their method compared to previous methods and provide theoretical guarantees for their method.", "questions": "1. Can you put your paper into a broader perspective? This paper only compares methods that are based on the Parareal algorithm. Are there also alternatives and how do they compare to your method?  \n2. What is the “Fine” algorithm in Table 1 and in Fig. 1? I could not find its definition.\n3. How can the Runtime in Figure 1 be negative? Is log(Runtime) shown in this plot?\n4. Is there a way to measure the solution quality of the solver? Is the solution quality fully defined by the converged initial condition $\\epsilon$ and therefore the same for all methods? \nIf the solution quality can be quantified, how does the solution quality of your method compare to other methods?\n5. Are the subscripts $k$ and $i$ dropped in Equation 6? If so this is a bit confusing since these are used in the previous Section and then dropped without a notice.", "limitations": "The authors addressed some limitations of their work.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The paper evaluates their method on three different PDEs.\nThe authors show that their method yields a significant speed-up.\nThe authors provide theoretical guarantees of their method.\nThe authors claim that the algorithm can be used as a convenient out-of-the-box algorithm.", "weaknesses": "Paper Structure:\nThe authors provide experimental details within the introduction (lines 85-90) that should belong to the experimental section.\nThe paper is not structured as most papers in the conference. I think Section 2 should be put into a Section called “Problem Description”. Section 3 should be included in Section 2 or into a “Related Work” section. The paper misses a related work section.\n\nTo me, it is questionable if this paper is of interest to the audience of the conference. The paper replaces a subpart of a numerical algorithm with a two-layer Neural Network where only the last layer is trained. Therefore, while the contribution is sound the impact of the work to other fields of machine learning is very limited.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate that you recognize how the proposed method yields a significant speedup compared to competing ones. We believe to have successfully replied and tackled all your questions (Q) and weaknesses comments (W). We hope that you reconsider your score based on the new provided information.\n\n\n**Q1: Lack of background on existing parallel-in-time (PinT) schemes**  \n**A**: This is a very important point, which we will address by adding the relevant background in the revised introduction. There are three general approaches for PinT computation: parallel across-the-problem, parallel-across-the-step, and parallel-across-the-method. [8,11] provide another classification: multiple shooting, methods based on waveform relaxation and domain decomposition, multigrid approaches, and direct time-parallel methods. Parallel-across-the-step methods, in which solutions at multiple time-grid points are computed simultaneously, include the Parareal (approximation of the derivative in the shooting method), Parallel Full Approximation Scheme in Space and Time (PFASST) (multigrid method) [5,10], and Multigrid Reduction in Time (MGRIT) [6,7] methods (see [9] for details). \n\nAmong them, Parareal has received more attention in the literature, with extensive theoretical analyses, improved versions, and empirical applications [8,11]. Limited theoretical results are available for MGRIT and PFASST, with a few extensions and empirical applications. Interestingly, combined analyses have shown equivalences between Parareal and MGRIT, and connections between MGRIT and PFASST.\n\nIt is well acknowledged that comparing PinT methods based on different working principles is extremely hard, with [11] representing a recent survey article which nonetheless does not offer a systematic comparison either. Quoting [11], \"caution should be taken when directly comparing speedup numbers across methods and implementations. In particular, some of the speedup and efficiency numbers are only theoretical in nature, and many of the parallel time methods do not address the storage or communication overhead of the parallel time integrator\". [9] is one of very few recent attempts to systematically compare different PinT classes. However, it is limited exclusively to the Dahlquist problem. Thus, it has become conventional to compare new techniques to the existing state-of-the-art methods within the same group of solvers. This is what we do in our paper, comparing RandNet-Parareal with the original Parareal and the recently improved version nnGParareal [1]. A broader comparison with alternative PinT approaches would be insightful, but it is beyond the scope of this work.\n\n **Q2: What is the “Fine” algorithm in Table 1 and in Fig. 1?**  \n**A**: The fine solver $\\mathcal{F}$ is the accurate solver defined on Page 1, line 28 in our submission. $\\mathcal{F}$ is typically chosen to be a higher accuracy method (e.g. Runge-Kutta 8) compared to the coarse solver (e.g. Runge-Kutta 2).\n\t\n**Q3: How can the Runtime in Figure 1 be negative? Is log(Runtime) shown in this plot?**  \n**A**: Yes, it is the $\\log_{10}$ runtime, as specified in the y-axis on Fig. 1 [paper].\n\t\n**Q4: Is there a way to measure the solution quality of the solver? How does the solution quality of the method compare to other methods?**  \n**A**: The fine solver $\\mathcal{F}$ is chosen by the user, and it determines the accuracy of the solution. All Parareal-type PinT schemes target such solution, with closeness controlled by $\\epsilon$, i.e., all converged solutions will be $\\epsilon$-close to the one of $\\mathcal{F}$. In the table below, we report the maximum absolute error committed by each of the RandNet-Parareal, Parareal, and nnGParareal with respect to the fine solver (run sequentially), averaged over intervals, together with the runtime in parentheses. \n\n|PDE|RandNet-Parareal|Parareal|nnGParareal|\n | - | - | - | - |\n |Burgers' $d=128$| $1.06e^{-8}$ (1h 2m)|$1.85e^{-8}$ (8h 54m)|$1.32e^{-7}$ (1h 39m)|\n |Diffusion-Reaction $d=7.2e^2$|$3.56e^{-8}$ (23m)|$1.83e^{-8}$ (1h 40m)|$5.71e^{-7}$ (1h 11m)|\n |Diffusion-Reaction $d=3.3e^3$|$8.56e^{-10}$ (33m)|$2.45e^{-8}$ (7h 52m)|not converged|\n |Diffusion-Reaction $d=2.5e^4$|$8.09e^{-11}$ (1h 57m) |$7.43e^{-9}$ (9h 50m)|not converged|\n |SWE $d=3.1e^4$|$6.75e^{-8}$ (4h 9m)|$5.15e^{-8}$ (15h 43m)|not converged|\n |SWE $d=6.1e^4$|$8.54e^{-9}$ (12h 34m)|$2.84e^{-8}$ (19h 30m)|not converged|\n \n**Q5: Are the subscripts $k$ and $i$ dropped in Eq. 6?**  \n**A**: Yes, they are indeed dropped, since we first present our method for a generic setting, and then use it for the specific inputs, adding the indices again on lines 237-238.\n\n**W1: Paper structure**  \nWe did not find any specific guidelines in terms of paper structure, and we have seen numerous papers with similar structure to ours. Note that what you call \"Related work\" is in fact Section 3, where we introduce two recent works improving on Parareal. We will nevertheless take your comments on board when preparing the camera-ready version.\n\n**W2: Lack of interest for the audience of the conference**  \nWe disagree with you, as this is not the first contribution on PinT methods appearing on NeurIPS. You write that *\"while the contribution is sound, the impact of the work to other fields of machine learning is very limited\"*. However, this is not a NeurIPS requirement, as long as the paper matches the subjects and topicality of the Primary Area. Nevertheless, we believe that our new randomized PinT numerical solver has far-reaching implications for the scientific computing field and others, where efficient numerical (O/P)DEs solvers are in need. We kindly refer you to our detailed answer to W1 of Ref awZW. There, we discuss the importance of PinT schemes for scientific machine learning (ML)/computing and science in general and the use of scientific PinT computing in other fields of machine learning (e.g. diffusion models, neural ODEs, optimization, normalizing flows, to name a few).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PqRIj8Rq8L", "reviewer": "Reviewer_Eu5i", "summary": "The authors introduce a numerical algorithm that computes the solution to a large system of ordinary differential equations (ODE) \"parallel in time\". The main idea of the solver is based on the existing \"Parareal\", which introduces parallelism by running a sequential, fast, and inaccurate ODE solver and then corrects it in parallel with a more accurate correction scheme. The novelty in this manuscript is that this correction scheme is done using random neural networks, i.e., networks where the hidden weights and biases are chosen at random and then fixed. The benefit of this is that training time for the network is reduced significantly, as only the last layer must be approximated, which is possible using a linear solve. The authors demonstrate this efficiency and the parallel in time property on several nonlinear, time-dependent partial differential equations, discretized to obtain the required large ODE system.", "questions": "1) The authors stress that training a GP is fundamentally more expensive than training the random network. However, I fail to see this: if we do not compute a full kernel matrix of N by N (N being the number of data points), but just randomly choose M center points and then solve a linear system of N by M, the complexity is exactly the same as in the random network case. The only difference is that instead of M neurons we have M kernel functions (k(x_i, *)), essentially forming a \"radial basis neural network\". This has been discussed at length in the GP literature, and many methods for this exist (e.g. Nystroem kernel approximations, inducing points methods, etc.). Why do the authors only compare to the full kernel case, which obviously is subpar in performance?\n\n2) Why do the authors compare to clearly sub-optimal numerical solvers in the PDE setting? I may have missed the reasoning why even the simplest PDE settings take so much time to solve.\n\n3) Is Parareal the only parallel-in-time solver available? The paper focuses mostly on different Parareal settings and approaches, but not on different parallel-in-time solution methods.\n\n4) How does the high ambient-space dimension d of the state U affect the nearest neighbor search?", "limitations": "The authors are missing a discussion that \"nearest neighbor\" algorithms do not work well in high-dimensional settings (as given here, with $d \\gg 1000$ in many cases). They also do not discuss the suboptimal scaling of the linear solution for the random neural networks (caused by the linear solve necessary for the outer weights).", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Using random neural networks in a parallel-in-time setting is (to my knowledge) a new approach, and helps to mitigate issues with earlier versions of algorithms (e.g. the GP version of ParaReal, as presented). To that end, the work is a combination of well-known techniques to solve a problem in a new way. The authors discuss and demonstrate how their approach differs from earlier work, mostly from ParaReal with GP, fine-scale solvers, and a classical ParaReal approach. The paper is very well written, with clear descriptions of the algorithm and numerical results, and even includes scaling results for multi-core systems. The authors also include a robustness study of the approach in the appendix.", "weaknesses": "1) There is no theoretical analysis beyond restating existing results from the literature. For example, there is no complexity analysis for the random network setting, which should have been relatively straightforward (given that it involves only a single linear solve, for which complexity results are available).\n2) The numerical scenarios chosen in the paper are not particularly challenging on their own. It seems to me that the parallel in time setting itself is what makes them challenging, not the high number of degrees of freedom (e.g., 10^5 is not very large), nor the particularly high spatial dimension (one and two), or the complexity of the PDE (mildly nonlinear). Of course, for parallel-in-time algorithms these are challenging problems nonetheless, but the examples do not demonstrate why parallel in time in general makes more sense than just solving the system with a better algorithm (e.g. higher order, other basis functions, etc.). For example: it is not reasonable that solving Burgers equation on a 1D line takes 13 hours on the \"fine scale\" (table 1). This looks like an unnecessary difficulty for a problem that should be simple to solve, especially with only 128 degrees of freedom (variable d).\n3) It seems to me that the main novelty of this paper is to replace the (numerically suboptimal) Gaussian process with a neural network where the internal weights are randomly chosen. It is not clear to me why this particular change is any better than improving the numerics for the GP setting, or using any other possible solution (e.g., just a linear mapping as approximation, or running a fine-scale solver on multiple scales, or using polynomial regression, ...). It is probably not possible to perform these additional experiments during the rebuttal time, but the manuscript does not even state any other possibilities that were ruled out.\n4) Similarly, the state of the art focuses very much on ParaReal, and not on any other parallel in time solver. While this may be acceptable if ParaReal was the most advanced PinT solver available, it is not clear (because it is not stated) if that is the case in general (or at least for the chosen examples).\n\nA few minor issues:\n1) I would avoid using \"nns\" for nearest neighbors (l.35), because it is confusing with \"NNs\" (neural networks).\n2) l.87: it is strange to say \"24 core processors, 48 cores\". It is clear that 24*2=48; I hope that is what the authors mean.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for the critical assessment of our work. We identified these main criticisms: (a) lack of a computational complexity study; (b) comparison to a GP with full kernel and lack of alternatives to the GP solution; \\(c\\) use of a suboptimal fine solver _F_. We addressed (a) conducting a detailed complexity study. We believe that (b)-\\(c\\) are based on a misconception and provide our explanation together with answers to questions (Q), weaknesses (W), and limitations (L). We hope it is evident that the points raised do not warrant a rejection of our work. \n\n**W1: Lack of complexity analysis**\n- The complexity of RandNet-Parareal is reported in Tbl. A [pdf] and is linear in $d$, while nnGParareal scales as $d^2$.\n- The theoretical computational costs as a function of $d$ and $N$ are plotted in Fig.A-B [pdf], confirming the empirical superior scalability of RandNet-Parareal given in Fig.1 [paper].\n- The speed-up of each approach can be immediately obtained as the ratio of the sequential runtimes of _F_ and the parallel runtimes, as in [1].\n\n**Q1: The proposed approach is only compared to a GP based on full-rank matrix**  \nThis statement is not correct. The competing approach does not use the GP with the full kernel matrix of $N\\times N$ ($N$ is the sample size), but only its reduced $m \\times m$ version with $m \\ll N$ nearest neighbors (nnGP [1]). This is the only existing PinT scheme tackling the full GP kernel issue.\n\n**Q2 & W2: Choice of subobtimal numerical solvers and high running times**  \nAlthough valuable for sequential solvers, these comments are not pertinent for PinT schemes. We will add the following comment to a camera-ready version for a better understanding of the customary empirical design for PinT schemes.\n- *Goal of PinT schemes:* There is some confusion regarding this point. The essence of any Parareal method is to obtain in a PinT manner the $\\epsilon$-close solution to that produced sequentially by a generic _F_. Each Parareal variant is agnostic to the choice of _F_, and \"canonical\" off-the-shelf numerical solvers are standard in the PinT literature (here, Runge-Kutta, as in [1]). WLOG, more appropriate (accurate/faster) _F_ can be used. This does not contrast our framework: RandNet-Parareal allows the seamless replacement of _F_ and inherits its properties.\n- *Contribution and choice of slow _F_*: As we propose a new PinT scheme, we compare it with existing PinT approaches rather than ad-hoc sequential solvers, focusing on speed-up rather than absolute performance of _F_. Burgers' PDE is particularly slow as we used the same high accuracy setup of [1].\n- *Additional challenging example*: Despite the short rebuttal time, we ran 2 more challenging examples (2D and 3D Brusselator PDE) with results in Fig.C [pdf]. Should we add them to Suppl. Material?\n\n**W3: Choice of RandNet instead of better numerics for the GP/alternative solutions**  \nA recent work [1] replaced the numerically suboptimal GP from [3] with nearest neighbor GPs. While improving GP numerics, it still suffers from quadratic complexity in $d$ (see Q1) and is sensitive to hyperparameter tuning (Tbl. A [pdf]). We propose abandoning the GP framework and using RandNets.\n\nAlternative choices, e.g. linear/polynomial functions, do not offer the desirable learning quality:\n - *Universal approximation (UA) properties:* Linear functions do not form a UA class, so do not guarantee high-accuracy learning. Polynomials are UA, but (1) their order depends on the problem; (2) they are known to yield ill-conditioning in the regression setting. Instead, RandNets (as neural network functions) are UA [4] and do not suffer from (1)-(2).\n- *Overcoming the curse of dimensionality (CoD):* Under specific assumptions, one can determine if RandNets suffer from the CoD [4] (see also Q3 of Ref nFhC). This opens multiple avenues for provable speed and uncertainty quantification improvements with RandNet-Parareal.\n\nDespite limited rebuttal time, we tested $p$-order polynomials on Diffusion-Reaction PDE, finding they lack sufficient learning accuracy ($p=1$ is linear). The table displays iterations to convergence $K$ and parallel speed-up (in parentheses). For all models $K=N=64$ (so they converge serially), while for RandNet-Parareal, $K=12$ with a speed-up of 5.36 times (Tbl. 6 [paper]).\n|PDE|$p=1$|$p=2$|$p=3$|$p=5$|$p=7$|\n|-|-|-|-|-|-|\n|$d=722, N=64$|64 (0.98)|64 (0.95)|64 (0.92)|64 (0.88)|64 (0.83)|\n\n**W4 & Q3: Is Parareal the state-of-art among PinT schemes?**  \nYes. Due to space constraints, we kindly refer you to our detailed answer to Q1 of Ref FREq who posed a similar question. We will add the needed background to the revised introduction.\n\n**W5: nns (nearest neighbors) vs NNs (neural networks)**  \nThis may indeed be confusing. In the literature, nearest neighbor GPs are typically denoted as nnGPs, while neural networks as NNs. We welcome alternative suggestions.\n\n**W6: 24 core processors, 48 cores**  \nIn the updated version, we will change this to read as \" 24 core processors, _yielding a total of_ 48 cores.\"\n\n**L1 & Q4: Issues of nearest neighbor algorithms in high dimensions**\n- *Suboptimal scaling of the linear solution for RandNets:* As shown in Tbl. A [pdf], it contributes as $M^3$ ($M$ the number of neurons) and has limited effect on performance for sensible $M$ values for our examples of (O/P)DE, spatial dimension, degrees of freedom.\n- *Complexity in high dimensions:* Naive linear search of the nearest neighbor is linear in $N$ and $d$, $O(Nd)$, which we took into account in complexity in WC1. Approximate Near. Neighbors would allow for further improvements [2].\n- *Performance in high dimensions:* Fig.D [pdf] shows that nearest neighbors (closest points in L2 distance) of target points are recovered for 2 high-dimensional systems $d>1e^4$, proving the validity of our approach. This nn structure is a feature of Parareal -- the dataset is _not made of random_ observations but of initial conditions converging in $k$.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "zemTQPw3bG", "reviewer": "Reviewer_nFhC", "summary": "This paper proposes a method to accelerate the simulation of partial differential equations (PDEs) by converting them into systems of ordinary differential equations (ODEs). It then utilizes a framework that merges the random neural network and the parareal approach, termed RandNet-Parareal. For validation, three complex systems of PDEs are considered: the Burgers equation, the diffusion-reaction equation, and the shallow water equation. Results show gains achieved in terms of computational cost compared to other methods.", "questions": "1. The proposed method's performance and compared methods are not quantified through a metric. It would be interesting to see the trade-off between accuracy and computational cost.\n\n2. The method is not compared with other neural PDE solvers, specifically random neural network-based PDE solvers aiming to accelerate PDE simulations. Assessing the gain achieved compared to the methods proposed in the literature would be interesting.\n\n3. As the method depends on traditional mesh-based numerical solvers, can the authors provide insights on how the method would perform in complicated geometries, particularly in higher dimensions? \n\n4. Can the authors comment on the performance of their method when using a deep network instead of a shallow network?", "limitations": "The authors have mentioned the limitations of the work related to coarse numerical solvers and stiff systems. There are no social or ethical issues concerning the paper.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1. The proposed method accelerates neural PDE simulations by combining random neural networks with parareal concepts.\n\n2. The considered PDEs are prototypical and complicated and exhibit the efficacy of the proposed approach.", "weaknesses": "1. The authors have not mentioned the trade-off between accuracy and computational cost, which is essential for high-fidelity simulations.\n\n2. State-of-the-art comparisons are lacking, restricting the evaluation of the proposed method with random neural network-based PDE solvers.\n\n3. The method is based on conventional numerical solvers, and hence, extending its applicability to complex geometries in higher dimensions may suffer from the curse of dimensionality.\n\n4. It is often the case that random neural networks do not scale well to deep networks, and not much gain is achieved when using a deep network. The paper does not discuss extending the method to deep neural networks and what gains could be achieved.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank you for acknowledging that our method succeeds in accelerating the solutions of complex PDEs. We appreciate you recognize the efficacy of our approach and see no limitations beyond those we mentioned in the paper.\n\nIn the paragraphs below, we replied in detail to your questions (Q) and weaknesses (W) grouped by topic. We hope to convince you that our paper merits a higher score, as it pushes the boundaries of numerical solvers in the field of scientific computing by proposing a new direction of randomized numerical parallel-in-time (PinT) (O/P)DE solvers.\n\n**W1 & Q1: Discussion of the trade-off between accuracy (w.r.t. some metric) and computational cost, essential for high-fidelity simulations**  \n**A**: Thank you for the question similar to Q3 of Ref FREq and W1 of Ref Eu5i. The fine solver _F_, chosen by the user, defines the accuracy of the competing methods. All Parareal PinT schemes target the solution provided by _F_. For any given $\\epsilon>0$, all the solutions for all converged algorithms are $\\epsilon$-close to the solution of _F_ (see eq.(4) in the paper, where another accuracy metric criterion could be used). \n\n_Empirical evidence_: In the table in our response to Q3 of Ref FREq (not copied here due to the space limit), we report the accuracies and runtimes for RandNet-Parareal, Parareal, and nnGParareal. Max. abs. error (mean across intervals) w.r.t. _F_ (run sequentially) is reported with the runtimes in parentheses. $\\epsilon=5e^{-7}$ is chosen for the examples in the paper. One can clearly see that RandNet-Parareal has the best accuracy-cost trade-off.\n\n_Theoretical results_: To further address your comments and W1 of Ref Eu5i on the computational costs, we derived the theoretical complexity of RandNets and compared it to nnGParareal (see pdf for the details). Figs. A-B plot the theoretical _model_ cost and _total_ cost (in $\\log_{10}$(hours)), respectively, across dimension $d$ (and cores/subintervals $N$). To calibrate the constants in both complexity bounds, we used the total empirical cost in Fig. 1 and its breakdown in Tbl. 6 [paper]. Fig. A shows significantly superior scalability of RandNet-Parareal w.r.t. nnGParareal, and Fig. B shows that once the cost of _F_ is added, our results are fully coherent with the empirical results in the paper and the table in our response to Q3 of Ref FREq.\n\n**W2 & Q2: Comparison with other state-of-the-art and neural PDE solvers, specifically RandNets-based ones, and assessment of gains achieved.**  \n**A**: We kindly refer you to our extensive responses to similar questions Q1 of Ref FREq and W2 & Q2 of Ref Eu5i, which we cannot repeat here due to the space limits. We detailed why we focus exclusively on PinT (i.e., Parareal, nnGParareal) and not sequential (neural PDE) benchmarks. Importantly, RandNet-Parareal is the first PinT scheme in general, and Parareal algorithm in particular, that exploits random NNs, so no other NN-based PinT benchmark could be included in our comparative study. As the RandNet-Parareal construction is valid for any fine solver, any (neural) (O/P)DE solver can be readily incorporated as _F_ WLOG.\n\n**W3 & Q3: The method is based on mesh-based numerical solvers. Comment on applicability to complex geometries in higher dimensions (curse of dimensionality (CoD)).**  \n**A**: We split the answer to this interesting question as follows:\n\n- _CoD in the learning_: Our goal is to approximate correction functions in high dimensions. NN functions are universal approximators (UA) (for example, they are dense in the set of continuous functions, defined on compacta, in the sup norm). Random NNs with random inner weights are proven to be UA of functions of certain regularity (related to Radon-wavelet integral representation in [4]; see more general setup in [12]). See also our reply to your W4 & Q4 below. Notably, under specific regularity assumptions on the correction function, it is possible to establish whether RandNets suffer or not from CoD (i.e., $M$ needs to grow exponentially in $d$ to maintain the same approximation accuracy) [4]. This opens avenues for provable (potentially not prone to CoD) speed and uncertainty quantification improvements with RandNet-Parareal, some of which we pursue in our ongoing work.\n- _CoD in the cost_: Albeit numerous attempts, the state-of-the-art Parareal algorithms still suffer from unfavourable on-dimension cost dependence. We argue in our answer to W1 of Ref Eu5i and our answer to your W1 & Q1 that our main competitor, nnGParareal, recently proposed in [1] to reduce the GP training cost, also has this issue (see pdf), while RandNet-Parareal needs no hyperparameter tuning and has linear-in-$d$ complexity. We note that the matrix $X$ of neurons activated with the ReLu activation can be sparse. This implies that RandNets' cost is further improvable as the complexity of sparse operations is proportional to the number of nonzero matrix entries. \n- _Complex geometries in higher dimensions_: _F_ can be a symplectic, variational, energy-preserving, or any other integrator preserving the system's geometric properties. Using RandNet-Parareal with such _F_ is straightforward.\n\n**W4 & Q4: Discussion of poor scaling of RandNets to deep architectures and possible gains of using a deep network instead of a shallow one**  \n**A**: The poor scaling of RandNets to deep architectures is not always observed. Statistical properties of RandNets (or random feature NNs) have been studied in a sizeable body of works [18,19]. Asymptotic characterizations for the test error of _shallow_ RandNets have been derived [15,16,17]. Similar studies in the deep case [13,14] show the equivalence of deep RandNets to deep linear Gaussian models and, importantly, provide instances where testing performance _improves_ as a function of depth. Nevertheless, deep RandNets require higher training cost, so using _shallow_ RandNets may still be preferable in our method. We explore this in our ongoing work.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "fOLL8vhfLW", "reviewer": "Reviewer_awZW", "summary": "This paper proposes a new method for sequentially predicting and correcting numerical simulations by introducing random neural networks. The RandNets are single-layer feed-forward neural networks and only the output layer is for training. The numerical experiments have shown the improved training efficiency and scalability of the proposed method compared to the baseline models.", "questions": "In lines 7-8, is it a character \"x\" or a symbol \"\\times\" for `x125`? Probably a symbol is preferred.", "limitations": "N/A.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- This paper is well-written and well-organized. The theoretical guarantee of RandNets-Parareal is provided.\n\n- Different types of PDEs have been used to evaluate the model performance.", "weaknesses": "- The studied problem is domain-specific, which might not be of general interest to the scientific machine-learning community. \n\n- In Sec 5 Numerical Experiments, this paper provides the comparison of speed-ups/runtimes. It would be better to also show the varying solution accuracy with different setups. \n\n- Sections 2 and 3 might be shortened and some of these parts can be moved to the Appendix. The Robustness study (Appendix C) can be moved to the main text since it highlights the benefits of the proposed method.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank you for recognizing that our method shows improved training efficiency and scalability compared to the baseline methods. We also appreciate you acknowledge that our paper is well-written and well-organized, contains theoretical developments, and exemplifies RandNet-Parareal on different types of PDEs.\n\nIn the paragraphs below, we replied in detail to all your commented weaknesses (W). We hope we can convince you that the paper merits a higher score, especially in view of our explanations of the importance of advancing the state-of-art on (parallel-in-time, PinT) numerical solvers and the implications of our work for the scientific computing community. We also advise you to look at the additional theoretical results (computational complexity, Tbl. A [pdf]) and empirical evidence (two more challenging examples, 2D and 3D Brusselator PDE, Fig. C [pdf]) we added in response to other referees' comments and questions during the time of the rebuttal, which we will incorporate in the camera-ready version of the paper. \n\n**W1: The studied problem is domain-specific, which might not be of general interest to the scientific machine-learning community.**  \n**A**: We divide our answer into the following two parts:\n- _Importance of PinT schemes for scientific machine learning (ML)/computing and science in general_: Time and space efficient solving of (O/P)DEs remains one of the most important research directions of science and engineering. Broader availability of more affordable multicore high-performance computing resources makes parallelization more attractive. Much interest in the literature is concentrated on new numerical schemes, which allow for sub-linear in-time parallelization surpassing slow sequential solvers, especially crucial for long time horizons. Innovation in the area of PinT schemes is actively supported by research funding agencies (e.g. the ExCALIBUR cross-cutting research project \"Exposing parallelism: Parallel in Time\" involving the Met Office and UKAEA). Our proposed PinT algorithm leverages approximation and generalization properties of a flexible ML random NNs class, advancing the state of research on numerical Parareal solvers. Our work matches the NeurIPS topic of Primary Area on \"Machine learning for physical sciences\", opens the door to future research on randomized PinT numerical methods, and has implications for scientific computing, physics, and beyond [20,21,24]. Potential extensions of RandNet-Parareal include solvers for SDEs, coupled ODEs, and PDEs with noisy boundary conditions. We pursue some of these directions in our ongoing work.\n\n- _Scientific PinT computing for ML_: efficient ODE solvers started to gain much attention in the context of Neural ODEs, which are used as continuous-time models of certain NNs [22]. Other notable examples of ML models where PinT solvers could be game-changers are diffusion and generative models [23], for which the developed dedicated solvers could be further parallelized in time. \n \n**W2: In Sec. 5 Numerical Experiments, this paper provides the comparison of speed-ups/runtimes. It would be better to also show the varying solution accuracy with different setups.**  \n**A**:  Thank you for the question, similar to Q3 of Ref FREq and W1 of Ref Eu5i. The fine solver $\\mathcal{F}$, chosen by the user, defines the accuracy of the competing methods. All Parareal PinT schemes target the solution provided by $\\mathcal{F}$. For any given $\\epsilon>0$, all the solutions for all converged algorithms are $\\epsilon$-close to the solution of $\\mathcal{F}$ (eq. (4) [paper]). $\\epsilon=5e^{-7}$ is set for all the examples considered.\n \n_Empirical evidence_: In the table below (as in our response to Q3 of Ref FREq), we report the accuracies and runtimes (shown in parentheses) for RandNet-Parareal, Parareal, and nnGParareal. The accuracy is measured as max. abs. error (mean across intervals) w.r.t. $\\mathcal{F}$ (run sequentially). RandNet-Parareal has the best accuracy-cost trade-off, and one can easily see that under the same time budget for all approaches, RandNet-Parareal achieves higher accuracies than benchmarks.\n\n|PDE|RandNet-Parareal|Parareal|nnGParareal|\n |- | -| -|- |\n |Burgers' $d=128$|$1.06e^{-8}$ (1h 2m)|$1.85e^{-8}$ (8h 54m)|$1.32e^{-7}$ (1h 39m)|\n |Diffusion-Reaction $d=7.2e^2$| $3.56e^{-8}$ (23m)| $1.83e^{-8}$ (1h 40m) |$5.71e^{-7}$ (1h 11m)|\n |Diffusion-Reaction $d=3.3e^3$|$8.56e^{-10}$ (33m)| $2.45e^{-8}$ (7h 52m)|not converged|\n |Diffusion-Reaction $d=2.5e^4$|$8.09e^{-11}$ (1h 57m)| $7.43e^{-9}$ (9h 50m)|not converged|\n |SWE $d=3.1e^4$| $6.75e^{-8}$ (4h 9m)| $5.15e^{-8}$ (15h 43m)|not converged|\n |SWE $d=6.1e^4$| $8.54e^{-9}$ (12h 34m)| $2.84e^{-8}$ (19h 30m)|not converged|\n\nTo address comments from other Referees, we derived the theoretical complexity of RandNets and compared it to nnGParareal (see Table A [pdf]). Figs. A-B [pdf] plot the theoretical _model_ and _total_ costs (in $\\log_{10}$(hours)), respectively, across dimension $d$ (and cores/subintervals $N$). To calibrate the constants in complexity bounds, we used the total empirical cost in Fig. 1 [paper] and its breakdown in Tbl. 6 [paper]. Fig. A shows significantly superior scalability of RandNet-Parareal w.r.t. nnGParareal. Fig. B shows that with the cost of $\\mathcal{F}$ added, our results are fully consistent with those presented in the paper and in the table in our response to Q3 of Ref FREq.\n\n**W3: Sections 2 and 3 might be shortened and some of these parts can be moved to the Appendix. The Robustness study (Appendix C) can be moved to the main text since it highlights the benefits of the proposed method.**  \n**A**: Thank you for this reasonable and nice suggestion. We hope to have a chance to introduce this change in the camera-ready version of the paper.\n\n**Minor Q. In lines 7-8, is it a character ``x'' or a symbol ``$\\times$'' for x125? Probably a symbol is preferred.**  \n**A**: Thank you for mentioning this. We will correct it accordingly.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We would like to thank the reviewers for their constructive feedback and suggestions. We are happy to read that the paper "is well-written and well-organized" (awZW), "very well-written" (Eu5i), "contains clear descriptions and robustness study" (Eu5i) and that the proposed method is a "new approach" (Eu5i) with "provided theoretical guarantees" (FREq) which "yields a significant speed-up" compared to previous methods (FREq) and "helps to mitigate issues with earlier versions of algorithms" (Eu5i). We are particularly encouraged by the observation that "considered PDEs are ... complicated and exhibit the efficacy of the proposed approach" (nFhC). We are glad that all reviewers agree on our contribution being the first one that embeds NNs within parallel-in-time (PinT) schemes. We hope this work will prompt further research in this area, ideally combining the recent advancements in neural PDEs solvers with PinT schemes.

We are grateful to the reviewers for suggesting several improvements, which will notably enhance the revised version of the paper. Two common objections are (i) the lack of background on alternative parallel-in-time (PinT) schemes and (ii) a discussion on the computational complexity of the proposed approach (with respect to both competing methods and accuracy). We addressed point (i) in the rebuttal comments, while for (ii) we derived the theoretical complexity of our method and compared it to the recently introduced nnGParareal (see Table A in pdf). 

Additionally, during the rebuttal period, following suggestions of Eu5i, we tested our method on two new challenging examples, the 2D and 3D Brusselator PDEs, which are known to exhibit complex behavior, including oscillations, spatial patterns, and chaos. We also compare RandNets with polynomial functions, showing the latter's suboptimal performance.

The attached pdf (cited as [pdf]) includes additional figures and the complexity table. In the camera-ready version, we will place our contribution within existing PiT schemes, add a discussion on computational complexity, and include results obtained during the rebuttal time in Supplementary Material if the reviewers deem it advisable.

We believe that our new randomized PinT numerical solver has far-reaching implications for the scientific computing field specifically and for science/engineering more broadly, where time- and space-efficient solving of (O/P)DEs remains one of the most important research directions. By combining randomized neural networks with existing PinT methodology, we provide a framework that couples theoretical and robustness guarantees with a drastic reduction of computational requirements, as demonstrated through extensive simulation. Finally, our method can offer speedups for numerous instances where solving ODEs in the machine learning realm is needed (for example, in the context of diffusion models, neural ODEs, and many others). 

We are deeply grateful for the helpful feedback received from the reviewers, which will make the final version of our manuscript (especially its introduction) more comprehensive and refined. 

**References used in the rebuttals:**

[1] G. Gattiglio, L. Grigoryeva, M. Tamborrino. Nearest neighbors GParareal: Improving scalability of Gaussian processes for parallel-in-time solvers. 2024.  
[2] P. Indyk, R. Motwani. Approximate nearest neighbors: towards removing the curse of dimensionality. 1998.  
[3] K. Pentland, M. Tamborrino, T. Sullivan, J. Buchanan, L. Appel. GParareal: a time-parallel ODE solver using Gaussian process emulation. 2023.  
[4] L. Gonon, L. Grigoryeva, J.-P. Ortega. Approximation bounds for random neural networks and reservoir systems. 2023.  
[5] M. Emmett, M. Minion. Toward an efficient parallel in time method for partial differential equations. 2012.  
[6] R. Falgout, S. Friedhoff, T. Kolev, S. MacLachlan, J. Schroder. Parallel time integration with multigrid. 2014.  
[7] S. Friedhoff, R. Falgout, T. Kolev, S. MacLachlan,  J. Schroder. A multigrid-in-time algorithm for solving evolution equations in parallel. 2012.  
[8] M. Gander. 50 years of time parallel time integration. 2015.  
[9] M. Gander, T. Lunet, D. Ruprecht, R. Speck. A unified analysis framework for iterative parallel-in-time algorithms. 2023.  
[10] M. Minion. A hybrid parareal spectral deferred corrections method. 2010.  
[11] B. Ong, J. Schroder. Applications of time parallelization. 2020.  
[12] A. Neufeld, P. Schmocker. Universal approximation property of random neural networks. 2023.  
[13] D. Schröder, H. Cui, D. Dmitriev, B. Loureiro. Deterministic equivalent and error universality of deep random features learning. 2023.  
[14] D. Bosch, A. Panahi, B. Hassibi. Precise asymptotic analysis of deep random feature models. 2023.  
[15] S. Mei, A. Montanari. The generalization error of random features regression: Precise asymptotics and the double descent curve. 2022.  
[16] H. Hu, Y. Lu. Universality laws for high-dimensional learning with random features. 2022.  
[17] F. Gerace, B. Loureiro, F. Krzakala, M. Mezard, L. Zdeborova. Generalisation error in learning with random features and the hidden manifold model. 2020.  
[18] J. Lee, Y. Bahri, R. Novak, S. Schoenholz, J. Pennington, J. Sohl-Dickstein. Deep neural networks as Gaussian processes. 2018.  
[19] A. De G. Matthews, J. Hron, M. Rowland, R. Turner, Z. Ghahramani. Gaussian process behaviour in wide deep neural networks. 2018.  
[20] F. Hamon, M. Schreiber, M. Minion. Parallel-in-time multi-level integration of the shallow-water equations on the rotating sphere. 2020.  
[21] D. Samaddar, D. Coster, X. Bonnin, L. Berry, W. Elwasif, D. Batchelor. Application of the parareal algorithm to simulations of ELMs in ITER plasma. 2019.  
[22] R. Chen, Y. Rubanova, J. Bettencourt, D. Duvenaud. Neural ordinary differential equations. 2018.  
[23] G. Papamakarios, E. Nalisnick, D. Rezende, S. Mohamed, B. Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. 2021.

Meta-review: This paper continues a line of work on *Parareal methods for accelerating PDE solvers on parallel machines. The results demonstrate significant speedup on a variety of canonical problems. The innovation is somewhat modest, employing well-studied random neural network as a drop-in replacement for Gaussian processes used in previous work.  Nevertheless, I think these results will be of interest to some in the community.
Decision: Accept (poster)