RandNet-Parareal: a time-parallel PDE solver using Random Neural Networks

Authors: guglielmo  gattiglio
lyudmila  grigoryeva
massimiliano  tamborrino
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-04 16:22:43

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.

Main paper: 
 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 ∈ N ODEs (and similarly for PDEs)EQUATIONwhere h : R d × [t 0 , t N ] → R d is a smooth multivariate function, u : [t 0 , t N ] → R d is the time dependent column vector solution, and u 0 ∈ R d is the initial value at t 0 . PinT schemes are particularly important when the sequential application of an accurate numerical integrator F over [t 0 , t N ] is infeasible in a reasonable wallclock time. There are three general approaches for PinT computation: parallel across-the-problem, parallel-across-the-step, and parallel-across-the-method. In [17, 55] , another classification is provided: multiple shooting, methods based on waveform relaxation and domain decomposition, multigrid approaches, and direct time-parallel methods. Parallel-across-thestep methods, in which solutions at multiple time-grid points are 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, 16] methods (see [19] for details). Among them, Parareal [45] has garnered popularity, with extensive theoretical analyses, improved versions, and empirical applications [17, 55] . 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, 56] , 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 G is run sequentially to obtain a first approximation of the solution, which is then corrected by running a fine (accurate) but slow integrator F in parallel across N temporal subintervals. This procedure is then iterated until a convergence criterion is met after k ≤ 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 F -G using Gaussian Processes (GPs). This method outperforms Parareal for low-dimensional ODEs and a 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 scalability properties in both N and d through data reduction. Significant computational gains were achieved by training the GP on a small subset of nns, resulting in an 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 a costly hyperparameter optimization procedure that requires fitting one GP per ODE dimension.This study introduces RandNet-Parareal, a new approach using random neural networks (RandNets) to learn the discrepancy F -G . RandNets are a family of single-hidden-layer 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 closedform 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 [6, 49, 62, 63, 65] and Reservoir Computing [24, 26, 25, 27, 28] , Random Fourier Features (RFFs) and kernel methods [41, 61, 70, 74] . Some authors use the name Extreme Learning Machines (ELMs) [34] [35] [36] [37] 44] to refer to RandNets, while others use the term randomized or random NNs [5, 32, 39, 46, 78, 82] for the same paradigm. RandNets show excellent empirical performance, and have been used in the context of mathematical finance [22, 33, 38] , mathematical physics [52] , electronic circuits [69] , photonic [47] and quantum systems [23, 48] , random deep splitting schemes [53] , scientific computing [10, 11, 79, 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 a large class of Bochner spaces in [52] . These contributions show that RandNets are a reliable machine learning paradigm with provable theoretical guarantees.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 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". [19] 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, we compare RandNet-Parareal 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 3 briefly explains GParareal and nnGParareal, focusing on the latter. RandNet-Parareal is introduced in Section 4, while Sections 5 and 6 present our numerical results, and a final discussion. 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.Notation. We denote by v ∈ R n a column vector with entries v i , i ∈ {1, . . . , n}, and by ∥v∥ and ∥v∥ ∞ its Euclidean and infinity norms, respectively. We use A ∈ R n×m to denote a real-valued n × m matrix, n, m ∈ N, with elements A ij , jth column A (•,j) , j ∈ {1, . . . m}, and ith row A (i,•) , i ∈ {1, . . . , n}. We write A ⊤ , A † , and ∥A∥ F for the A matrix transpose, Moore-Penrose pseudoinverse, and Frobenius norm, respectively. I n denotes the identity matrix of dimension n. 
the parareal algorithm: 
 The idea of Parareal is to solve the d-dimensional ODE (and similarly PDE) system (1) in a parallelin-time fashion, dividing the original IVP into N sub-IVPsdu i dt = h (u i (t | U i ) , t) , t ∈ [t i , t i+1 ] , u i (t i ) = U i , for i = 0, . . . , N -1,where the number of time intervals N is also the number of available machines/cores/processors,u i (t | U i ) is the solution at time t of the i th IVP with initial condition u(t i ) = U i ∈ R d , i = 0, . . . , N -1.If the initial conditions were known and satisfied the continuity conditions U i = u i-1 (t i |U i-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 U 0 = u 0 ∈ R d at time t 0 appears available. To account for this, Parareal introduces another numerical integrator G , much faster but less accurate than F , to approximate the missing initial conditions U i , i = 1, . . . , N -1, sequentially. G trades off accuracy for computational feasibility, usually taking seconds/minutes instead of hours/days of F2 .The algorithm works as follows. We use U k i to denote the Parareal approximation of u i (t i ) = U i at iteration k ≥ 0. At k = 0, the initial conditions {U 0 i } N -1 i=1 are initialized using a sequential application of the coarse solver G , obtainingU 0 i = G (U 0 i-1 ), i = 1, . . . , N -1, with U 0 0 = U 0 . At k ≥ 1, the obtained initial conditions U k-1 i-1 are "propagated" through F in parallel on N cores to obtain F (U k-1 i-1 ), i = 1, . . . , N . Note that for every initial condition U k-1 i-1 , we compute both F (U k-1 i-1 ), i.e. a precise evaluation of u i-1 (t i |U k-1 i-1), and G (U k-1 i-1 ), an inaccurate evaluation of the same term. Hence, we can interpret F and G as functions mapping an initial condition to the next one, thereby evolving (1) by one interval. We can then use their difference, (F -G )(U k-1 i-1 ), to correct the inaccuracy of G on future evaluations. This gives rise to the original Parareal predictor-corrector ruleU k i = G (U k i-1 ) + (F -G )(U k-1 i-1 ), with i = 1, . . . , N -1, k ≥ 1 [18], where the sequential prediction G (U k i-1 ) is corrected by adding the discrepancy F -G computed at the previous iteration k -1. However, this formulation can be changed to use data from the current iteration k [57] , and generalized to account for different ways of computing the discrepancy, leading to [21] EQUATIONwhere f :EQUATIONwhile 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 t L ≤ t N , if solutions across consecutive iterations have stabilized. That is, for some pre-defined accuracy level ϵ > 0, it holds thatEQUATIONOther stopping criteria are also possible [66, 67] . Converged Parareal approximations U k i , i ≤ L, are no longer iterated to avoid unnecessary overhead [12, 20, 21, 57, 58] . Instead, unconverged solution values U k i , i > L, are updated during future iterations by first running F in parallel and then using the prediction-correction rule (2) . The Parareal algorithm stops at some iteration K Para ≤ N when all initial conditions have converged, that is when (4) is satisfied with L = N -1 and thus K Para = k. Note that during every Parareal iteration k > 1, the "leftmost" fine solver evaluation F (U k L ) is either run from the outcome of a previous fine computationU k L = F (U k-1 L-1 ), or from a converged initial condition ∥U k L -U k-1 L ∥ ∞ < ϵ.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 G and evaluating f N times. A Parareal pseudocode is presented in Algorithm 1 in Supplementary Material A.Question: Does the paper discuss the limitations of the work performed by the authors? 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 out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a 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 scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a 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 important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
gparareal and nearest neighbors gparareal: 
 The performance of Parareal can be improved by a careful selection of f in (2), combined with a better use of the available information present at iteration k. Let D k denote the dataset consisting of N k pairs of inputs U j i-1 ∈ R d and their corresponding outputsEQUATIONWhile Parareal relies on one observation to construct the correction f in (3), GParareal and following works, including this one, use all the discrepancy terms F -G and information in D k to make their predictions. The idea of GParareal is to learn the mapEQUATIONGPara (U k i-1 )) ⊤ ∈ R d , and finally plugged into the predictor-corrector rule (2) . In particular, each GP prediction f(s) GPara (U k i-1) is obtained as the GP posterior mean µ (s)D k (U k i-1) ∈ R, computed by conditioning the corresponding GP prior on the dataset D k , i.e. f (s)GPara (U k i-1 ) = µ (s) D k (U k i-1). 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 (2) . Using all information stored in D k 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 G , 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 k) 3 ) at iteration k. This negatively impacts the algorithm'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 a smaller dataset of size m, D i-1,k , composed out of the m nns (in Euclidean distance) ofEQUATIONnnGPara (U k i-1 ), . . . , f(d) nnGPara (U k i-1 )) ⊤ , with f (s) nnGPara (U k i-1 ) = µ (s) D i-1,k (U k i-1 ), s = 1, . . . , d.Here, µD i-1,k ∈ R, s = 1, . . . , d, denotes the nnGP posterior mean computed by conditioning the corresponding GP prior on the reduced dataset D i-1,k of size m. Due to the decreased sample size, each nnGP covariance matrix can be inverted at a cost of O(m 3 ) 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 f nnGPara (U k i-1 ) is made, which are at most N -k at iteration k (as at least k intervals have converged at iteration k), yielding a combined O(d(N -k)m 3 ) 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 as 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 a m × 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 is why the algorithm has been proposed for highdimensional ODE and PDE systems with d ≤ N . We refer to Supplementary Material B and [21] for more details on nnGParareal, and to Algorithm 2 in Supplementary Material A for the pseudocode of the nnGP training. In the next section, we address the nnGParareal issues by introducing RandNets.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]Justification: Proposition 1 is clearly stated with all the required assumptions. The proof is not given, as we cite this result from [25] , adjusting their notation to match our, as we clearly mention. The derivation of the computational complexity analysis is provided with all relevant details in Supplementary Material C.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 crossreferenced. • 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 a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
random neural networks parareal (randnets-parareal): 
 In RandNet-Parareal, we propose to learn the map R d → R d , U → (F -G )(U ) via RandNets, obtaining the RandNet-Parareal correction f RandNet-Para , which we then use within the predictorcorrector rule (2) . Prior to that, we define how RandNets work in a general setting with input U ∈ R d and output or target Y ∈ R d . Later in the text we will go back to the input of interest U k i . Let M denote the number of hidden neurons, and H A,ζ W (U ) be a single-hidden-layer feed-forward neural network used to learn F -G , given byEQUATIONwhere A ∈ R M ×d is the matrix of random, non-trainable weights of the hidden layer, ζ ∈ R M is a random non-trainable bias vector, and W ∈ R M ×d is the matrix of trainable output weights.Here, σ : R M → R M denotes an activation function obtained as the componentwise application of a non-linear map σ : R → R which we choose to be ReLU σ(x) = max(x, 0) with x ∈ R, to satisfy the assumption of Proposition 1 below. The entries of A and ζ are randomly sampled from given distributions P A and P ζ , respectively, and kept fixed. After observing the dataset D k , the output weights W are obtained as the minimum ℓ 2 norm least squares (or simply min-norm least squares) estimator or as the solution of the following penalized empirical minimization problem:W D k = lim λ→0 arg min W ∈R M ×d    (U ,Y )∈D k H A,ζ W (U ) -Y 2 + λ ∥W ∥ 2 F    ,which is also called a "ridgeless" (interpolation) estimator [30] , and can be more compactly written asEQUATIONHere, X ∈ R N k×M is a matrix with (X (l,•) ) ⊤ := σ(A(U (l,•) ) ⊤ + ζ), l = 1, . . . ,) ⊤ = U j i , (Y (l,•) ) ⊤ = Y j i , l = jN + i + 1, i = 0, . . . , N -1, j = 0, . . . , k -1. Whenever N k ≥ M and the rank of X ⊤ X ∈ R M ×M is M , (7) reduces to the standard least squares estimator W D k = X ⊤ X -1 X ⊤ Y , while if the rank of X ⊤ X is N k, the solution admits a closed form W D k = X ⊤ X † X ⊤ Y.We get inspired by [21] , where only m nns are used in the training. In this setting, M ≫ N k = 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 P A and P ζ . In this work, we take the rows of the weight matrix A and the bias entries of ζ to be independent and uniformly distributed. For this case, the approximation bounds are available [25, Proposition 3], which we report below using our notation. Proposition 1 (Approximation bound, [25] , Proposition 3). Let H * : R d → R, U -→ H * (U ) be an unknown function we wish to approximate with H A,ζ W defined in (6) . Suppose H * can be represented as H * (U ) = R d e i⟨w,U ⟩ g(w)dw for some complex-valued function g on R d and allU ∈ R d with ∥U ∥ ≤ Q, where ⟨•, •⟩ is the inner product on R d . Assume that R d max 1, ∥w∥ 2d+6 |g(w)| 2 dw < ∞.For ρ > 0, suppose the rows of A are i.i.d. random variables with uniform distribution on B ρ ⊂ R d , the Euclidean ball of radius ρ around 0, and that the M components of ζ are i.i.d. uniform random variables on [-max(Qρ, 1), max(Qρ, 1)]. Assume that A and ζ are independent and let σ : R → R be given by σ(x) = max(x, 0). Then, there exist a R M ×d -valued random variable W and an explicit (see (33) in [25] ) constant C * > 0 such thatE ∥H A,ζ W (U ) -H * (U )∥ 2 ≤ C * M ,and for any δ ∈ (0, 1), the random neural networkH A,ζ W satisfies P R d ∥H A,ζ W (U ) -H * (U )∥ 2 µ U (dU ) 1/2 ≤ √ C * δ √ M ≥ 1 -δ.Our choice of P A and P ζ satisfies the conditions of Proposition 1 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 a good choice for the sampling distribution, which we follow. If no prior information were available, the common approach would have been to take P A ∼ Unif(-a, a) M ×d , P ζ ∼ Unif(-b, b) M , and optimize a, b ∈ R + via expensive cross-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-Parareal correction functionEQUATIONHere, the estimated weights W D i-1,k are obtained using the reduced dataset D i-1,k consisting of the m RandNet nns of U k i-1 , 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 to a closed-form ridgeless interpolation solution presents a substantial difference and improvement with respect to (nn)GPs. Moreover, expensive hyperparameter optimization is avoided in RandNets, addressing the other major pitfall of GParareal and nnGParareal. 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 RandNet-Parareal, 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-Parareal which we observe in the numerical experiments reported in Section 5.In Supplementary Material D, we study the robustness of RandNet-Parareal to changes in the number of nns m 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 a more accurate approximation of the correction F -G , and that a higher number of neurons M would reduce the prediction error of RandNets (as in Proposition 1). One may also suspect the algorithm to be sensitive to the particular sampling seed. Remarkably, our empirical findings demonstrate that these factors have a limited impact on the number of iterations needed by RandNet-Parareal to converge, which remains largely consistent (up to a few iterations) across different values and ODE/PDE systems, for sensible choices of m RandNet and M . For the end user, this eliminates the need of ad-hoc tuning, making the proposed RandNet-Parareal a convenient out-of-the-box algorithm.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: We have taken care of ensuring the reproducibility of all results through a 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 a GitHub repository with a step-by-step Jupyter notebook outlining RandNet-Parareal, and the necessary code to reproduced the results has been provided in Section 1 in the main text.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be 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 a 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 a specific model and empirical evaluation, it may be necessary to either make it 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 a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. 
numerical experiments: 
 In this section, we first compare the performance of Parareal, nnGParareal, and RandNet-Parareal on the viscous Burgers' equation (one spatial dimension and one variable, also considered in nnG-Parareal [21] ), to showcase Parareal and nnGParareal challenges as the number of space discretization and, correspondingly, the dimensions d, increases. Then, we consider the Diffusion-Reaction equation, a 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 a 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-Parareal, Parareal, and nnGParareal reported in Supplementary Material F.Let T F and T G be the time it takes to run F and G over one interval [t i , t i+1 ], respectively, and let N F and N G 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 S alg , defined as the ratio of the serial over the parallel runtime, i.e. S alg := N T F /T alg . S alg 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.Our initial example is a non-linear, one-dimensional PDE (illustrated in Figure 7 of Supplementary Material H) exhibiting hyperbolic behavior [68] , described by the equation EQUATIONwith initial condition v(x, t 0 ) = v 0 (x), x ∈ [-L, L], L > 0, and Dirichlet boundary conditions v(-L, t) = v(L, t), v x (-L, t) = v x (L, t), t ∈ [t 0 , t N ].We use the same setting and parameter values as in [21] . More specifically, we choose L = 1, diffusion coefficient ν = 0.01, and discretize the spatial domain using finite difference [15] and equally spaced points x j+1 = x j + ∆x, with ∆x = 2L/d and j = 0, . . . , d. We hence reformulate the PDE as a d-dimensional ODE system.In our first numerical experiment, we choose N = d = 128, v 0 (x) = 0.5(cos( 9 2 πx) + 1), t 0 = 0, and t N = 5.9 as in [21] , and consider G = RK1, F = RK8, N G = 4 and N F = 4e 4 , where RK1 stands for Runge-Kutta of order 1, and similarly for RK4 and RK8. The results, reported at the top of Table 1 , show how RandNet-Parareal converges in fewer iterations and has a higher speed-up than Parareal and nnGParareal. The difference in the model training costs is striking, with the nnGP's being approximately 700 times higher than that of RandNets, reducing thus its potential speed-up.As real-world (one-dimensional) problems would require a higher spatial discretization, we increase d by one thousand to d = 1128, keeping N fixed. Unlike assuming matching hardware resources to the system size (as implicitly done in [21] , where d = N ), we deliberately do not increase N to assess the algorithms' performances under constrained conditions. Instead, both time discretization numbers are increased to N F = 6e 5 and N G = 293 (resulting thus in longer T F and T G times) to account for the finer spatial mesh [43] . As observed from the bottom of Table 1 , as d/N > 1, nnGParareal's issues become more pronounced, as the d scalar GPs cannot be run all in parallel across the N processors, but need d/N = 10 runs instead, slowing down the algorithm. In contrast, RandNet-Parareal 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.We now turn to a more challenging case study. The Diffusion-Reaction equation [75] (illustrated in Figure 8 in Supplementary Material H) is a system of two non-linearly coupled variables, the activator u = u(t, x, y) and the inhibitor v = v(t, x, y), defined on a two-dimensional spatial domain as∂ t u = D u ∂ xx u + D u ∂ yy u + R u , ∂ t v = D v ∂ xx v + D v ∂ yy v + R v .Here, D u , D v are the diffusion coefficients for the activator and inhibitor, respectively, andR u = R u (u, v), R v = R v (u, v)are their reaction functions defined by the Fitzhugh-Nagumo equation [42] R u (u, v) = u -u 3 -c -v, R v (u, v) = u -v,where c = 5e -3 , D u = 1e -3 , and D v = 5e -3 . We take (x, y) ∈ (-1, 1) 2 and t ∈ [0, 20]. The initial condition u(0, x, y) is generated as standard Gaussian noise. We apply a no-flow Neumann boundarycondition D u ∂ x u = 0, D v ∂ x v = 0, D u ∂ y u = 0, D v ∂ y v = 0 for (x, y) ∈ (-1, 1) 2 .The spatial domain is discretized by the finite volume method [51] , resulting in a d = 2N x N y -dimensional ODE with N x and N y the number of space discretizations along x and y, respectively. The time integration is conducted with RK of variable order for G and F (see Table 6 in Supplementary Material G). As in the previous example, we conduct two experiments for this system, with speed-ups and runtimes reported in Figure 1 . In the first one, we increased d and N proportionately (with d/N ∈ [11, 13] ) while maintaining all other quantities (i.e. G , F , m nnGP , m RandNet ) fixed until N = 256. This scenario reflects a situation where more resources are allocated to solve larger problem sizes. In contrast, in the second experiment, N remains fixed at 512, with d increasing proportionately with N G to maintain algorithm stability. Moreover, F is chosen to be RK8, with N F automatically selected by the used Python library scipy [77] . This second setting simulates a scenario with constrained resources, where the user aims to solve the system using a finer spatial mesh. Finally, we focus on SWEs on a two-dimensional domain, described by a system of hyperbolic PDEs∂ t h + ∇hu = 0, ∂ t hu + ∇(u 2 h + 1 2 g r h 2 ) = -g r h∇b,where u = (u, v) represents the velocities in the horizontal u = u(t, x, y) and vertical v = v(t, x, y) directions, h = h(t, x, y) denotes the water depth, b = b(x, y) describes a (given) spatially varying bathymetry, and hu can be interpreted as the directional momentum components. The parameter g r describes the gravitational acceleration, while ∂ t f denotes the partial derivative with respect to time, and ∇f the gradient of a function f . Following [75] , we solve a radial dam break scenario where a Gaussian-shaped 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 on the simulation setup are given in Supplementary Material G.1.In this case, our algorithm also converges much faster than Parareal, with a speed gain of x1.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-Reaction, the improvements are remarkable. RandNet-Parareal takes up to 4-10 hours and 37 days less than the Parareal and sequential solver, respectively. 
discussion and limitations: 
 This study improves the scalability properties, convergence rates, and parallel performance of Parareal and a more recently proposed PinT solver for ODEs and PDEs, nnGParareal [21] . By replacing the Despite its excellent performance, RandNet-Parareal has limitations common to all Parareal algorithms, as its rate of convergence relies on the accuracy of the coarse solver G . Although neural networks can help mitigate the impact of suboptimal choices of G (as observed for GPs in (nn)GParareal), if the solver is mismatched for the system -for example, an unstable solver for a stiff ODE -RandNet-Parareal, similar to Parareal and (nn)GParareal, is likely to exhibit non-convergent behavior. It would then be of interest to investigate RandNet-Parareal's performance when using customized solvers tailored to specific systems, such as those outlined in Section 1 for the shallow water equation and the viscous Burgers' equation, which we defer to future research.This section provides pseudocodes for the implementation of Parareal (Algorithm 1), and the training procedure for learning the discrepancy F -G via nnGPs in nnGParareal (Algorithm 2), and RandNets in RandNet-Parareal (Algorithm 3).Algorithm 1: Parareal (generic)Input: Initial condition u 0 at time t 0 , number of intervals N Output: Converged initial conditions {U K i } N -1 i=1 , with K the number of iterations to convergenceRescale the ODE/PDE system such that each coordinate takes values in [-1, 1] L ← 1 U 0 0 = u 0 for i ← 1 to N -1 do U 0 i ← G (U 0 i-1 ) end for k ← 1 to N do Compute F (U k-1 i-1 ), i = 1, . . . , N in parallel for i ← L + 1 to N -1 do U k i ← G (U k i-1 ) + f (U k i-1 ) /* Update the initial conditions */ end Convergence checks for i ← L + 1 to N -1 do if ∥U k i -U k-1 i ∥ ∞ < ϵ then L ← L + 1 /*In this section, we provide more details on the nearest neighbors (nns) Gaussian process modeling, the mathematical expressions of the nnGParareal correction function f nnGPara , and the reduced dataset D i-1,k . These are not explicitly presented in the main text as they require additional notation, which we believe does not enrich the explanation. While the description of GPs presented here is for nnGParareal (and the corresponding nnGPs), it immediately generalizes to GParareal by replacing the reduced dataset D i-1,k with the full dataset D k . The interested reader can find more details in the original papers, [21] and [57] .Let the set of inputs U j i-1 ∈ R d and outputs (F -G )(U j i-1 ) ∈ R d , i = 1, . . . , N , j = 0, . . . , k -1, collected by iteration k, be denoted by U k and Y k , respectively. Now, define D i-1,k as the restriction of D k to the m nns ofU k i-1 in U k , namely D i-1,k := {(U (l-nn) U k i-1 , Y (l-nn) U k i-1 ), l = 1, . . . , m} ⊂ D k ,whereY (l-nn) U k i-1 = (F -G )(U (l-nn) U k i-1) ∈ Y k , andU (l-nn) U k i-1is the lth nn of U k i-1 in D k , i.e. the lth ordered statistics of the set formed out of Euclidean distances∥U j i-1 -U ′ ∥ between U j i-1 and any U ′ ∈ U k . That is, there exists U 1 , . . . , U l = U (l-nn) U k i-1 ∈ U k such that, for any U ′ ∈ U k , U ′ ̸ = U r , r = 1, . . . , l,we haveEQUATIONFinally, let U , Y ∈ R m×d be the matrices of input nns and outputs collected in D i-1,k , respectively.In nnGParareal, following the Bayesian framework, a GP prior is placed over the correction function F -G for each of the d coordinates asEQUATIONGP : R d → R is the prior mean function, taken to be zero for all s = 1, . . . , d, and K GP : R d × R d → R is the exponential prior variance kernel functionK GP (U , U ′ ) = σ 2 o exp(-∥U -U ′ ∥ 2 /σ 2 in ), with σ 2in and σ 2 o 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 f (s) nnGPara (U k i-1 ) ∈ R, s = 1, . . . , d, is obtained from the GP posterior mean µ (s)EQUATIONwhere K( U , U k i-1 ) ∈ R m is a vector of covariances between every input collected in U andU k i-1 defined as (K( U , U k i-1 )) r = K GP (( U (r,•) ) ⊤ , U k i-1), r = 1, . . . , m, and K( U , U ) ∈ R m×m is the covariance matrix, with (K( U , U )) q,r = K GP (( U (q,•) ) ⊤ , ( U (r,•) ) ⊤ ), r, q = 1, . . . , m. Here, σ 2 reg denotes a regularization term, also known as nugget, jitter, or regularization strength, which is added to improve the numerical stability when computing the inverse matrix, see [21] for further details. The hyperparameters θ := (σ 2 in , σ 2 o , σ 2 reg ) entering into the posterior mean and prediction (11) control the performance of the GP, and are optimized by numerically maximizing the marginal log-likelihood:EQUATIONwhere K(•, •) depends on θ through the kernel K GP , and det(A) denotes the determinant of a square matrix A. For a thorough treatment of Gaussian processes, including derivation of the likelihood and of the posterior distribution (which is Gaussian with mean as in (11) , see [80] .Consider the d-dimensional initial value problem (1) for some (O/P)DE. Let N be the number of subintervals (data points) at each kth iteration of the PinT algorithm. For any kth iteration of the scheme, a total of N k data points, each d-dimensional, are available. Here, we provide the computational cost of RandNet-Parareal, and compare it to that of nnGParareal, the state-of-the-art 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 ∈ R d among N k available points are found at a cost which is at most linear in the sample size, that is O(mN k) (for moderate dimensions d, one can get an improved cost O(m log(N k)), logarithmic in the sample size) [21] . Since our goal is to compare the computational complexities of nnGParareal and RandNet-Parareal as a function of d, we consider the worst-case complexity of the nns search.Given an input U k i-1 ∈ R d , i = 1, . . . , N at iteration k, the computational model cost of a prediction U k i-1 produced by all d models of m nnGP -nnGPs at iteration k via the predictor-corrector rule (2) with nnGParareal correction (11) and m nnGP nns is given in [21] asT nnGP (k) ≤ C nnGP N k(n start n reg d N ∨ 1)× ( m nnGP d B:=K(U,U k-1 i-1 ) ⊤ + m 2 nnGP d C:=K(U,U ) + m 3 nnGP D:=(B+σ 2 reg Im nnGP ) -1 + m 2 nnGP B • D + m nnGP d BD•Y + m nnGP N k nearest neighbors ) = C nnGP N k(n reg n start d N ∨ 1)(m 3 nnGP + m 2 nnGP + d(m 2 nnGP + 2m nnGP ) + m nnGP N k),with C nnGP being some constant that in general does depend on k, m nnGP , and d. Also, n reg and n start 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.5]), respectively. Furthermore, ∨ is the maximum operator, and the factor (n start n reg d/N ∨ 1) ≥ 1 follows from the fact that d independent nnGPs and hyperparameter optimization are parallelized over the N cores.In RandNet-Parareal, the correction term f RandNet-Para is modeled by the random weights neural network and evaluated as (8) . Again, only m RandNet nns (in Euclidean distance) to U k i-1 are used to construct the prediction, leading to the following computational model cost at iteration k:T RandNet (k) ≤ C RandNet N k 1 N (M dm RandNet X:=σ(A•U +ζ) + M 2 m RandNet Σ:=X•X ⊤ + M r 2 Σ † + M m RandNet d Σ † •X + M 2 d W :=Σ † X•Y + M dm RandNet W ⊤ •X + m RandNet N k nearest neighbors ) = C RandNet k(M r 2 + M 2 m RandNet + d(M 2 + 3M m RandNet ) + m RandNet N k),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 in m only if Σ is of full rank numerically, which is not observed empirically) and C RandNet is a constant independent on N , k, M , d, m RandNet . The factor 1/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, (n start n reg d/N ∨1) ≫ 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 to d ≈ 10 5 . It is easy to see that T nnGP is quadratic in dimension d, while T RandNet is only linear. This difference is mainly due to the factor (n start n reg d/N ∨ 1) ≥ 1 in T nnGP as opposed to 1/N in T RandNet . • Although M > m nnGP , M = 100 is sufficient for consistent performance across a 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 optimizationd=722 N=64 d=1568 N=128 d=3362 N=256 d=11856 N=512 d=25538 N=512 d=53792 N=512 d=110450 N=512Panel A problems, both n reg and n start should be set large. Second, each loglikelihood maximization conducted per dimension of the system requires a large number of iterations performed sequentially. Hence, C nnGP ≫ C RandNet , with C nnGP depending, in general, on k, m nnGP , d, as opposed to C RandNet . Indeed, RandNet requires no tuning (a 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 to the use of the ridgeless estimator. Empirically, we observed C nnGP /C RandNet to be up to 1000 (this can be seen in Figure 1 ).• We emphasize that since the ReLu function is chosen as activation in RandNet, the matrix X of activated neurons is sparse with sparsity degree γ. Hence, the computational complexity T RandNet could be further improved, as 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 not use sparse operations in our code implementation.• The upper bound of T RandNet could potentially be improved further, as additional parallelization may occur during standard matrix operations, depending on the specific computing environment.Figure 3 illustrates the theoretical model costs T nnGP and T RandNet (Panel A) and theoretical total costs obtained by adding the coarse and fine solver costs (Panel B), as functions of the dimension d (and the corresponding N ). The results are reported in terms of log 10 (hours). To calibrate the constants in both complexity bounds, we used the total empirical computational cost in Figure 1 , together with its breakdown described in Table 8 . Panel A shows that RandNet-Parareal displays significant improvement in scalability with respect to the state-of-the-art Parareal algorithm nnG-Parareal, while Panel B demonstrates that whenever the cost of the fine solver is added, our results are in full coherence with the empirical results.In this section, we study the robustness of RandNet-Parareal to changes in the number of nns m RandNet , 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 K RandNet-Para to convergence for RandNet-Parareal 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-Reaction equation with d = 722, conducting 100 weight samplings for each system. For every set of weights, we iterate RandNet-Parareal across m RandNet values ranging from 2 to 20, and M values ranging from 20 to 500 in increments of 10. The proportions of iterations needed to converge across 100 runs for different values of m RandNet and M for the Burger's and diffusion-Reaction equations are reported in Figures 4 and 5 , respectively. Although we observe some minor differences between the two systems, the main trend is clear: as long as reasonable values of m RandNet and M are chosen, the iterations to convergence for RandNet-Parareal vary at most by a few units when changing the values of m RandNet , M or a particular sampling seed of weights. Nevertheless, larger M might improve the performance, since, in this case, RandNets operate in the interpolation regime, as discussed in Section 4. Here, we carry out an additional scalability study for the 2 and 3 spatial dimensional Brusselator PDE. This model is a two-component reaction system that exhibits complex behavior, including oscillations, spatial patterns, and chaos. It is described by∂ t u = D 0 ∇ 2 u + a -(1 + b)u + vu 2 ,and∂ t v = D 1 ∇ 2 v + bu -vu 2 .In chemistry, the components u, v refer to the concentration of two substances, whereas the constants D 0 , D 1 are the respective diffusivity of each component, indicating the rate at which the substances spread out in space. Moreover, the parameters a and b are related to reaction rates. In our experiments, we used D 0 = 0.1, D 1 = 0.1D 0 , a = 1, and b = 3. We take t ∈ [0, 35], (u, v) ∈ (-1, 1) 2 ×(-1, 1) 2 for the 2D Brusselator, and (u, v) ∈ (-1, 1) 3 × (-1, 1) 3 for the three spatial dimension case. We initialize the u values at time t = 0 by setting them equal to a, and the v values by taking them normally distributed over the spatial grid. Further details regarding the number of spatial Nu = Nv d G G∆t F F∆t N (u, v, t) ∈ (-1, 1) 2 × (-1, 1) 2 × [0, 35] 32 2048 RK1 0.034 RK4 1e -7 512 (u, v, t) ∈ (-1, 1) 2 × (-1, 1) 2 × [0, 35] 64 8192 RK1 0.033 RK4 1e -7 512 (u, v, t) ∈ (-1, 1) 3 × (-1, 1) 3 × [0, 35] 20 16000 RK1 0.052 RK4 1e -7 512 (u, v, t) ∈ (-1, 1) 3 × (-1, 1) 3 × [0, 35]25 31250 RK1 0.057 RK4 1e -7 512 N u and N v are the number of spatial discretization points for u and v along each spatial dimension, yielding a d = 2N 2 x -or d = 3N 3 x -dimensional ODE, depending on the considered system. G and F denote the coarse and fine solvers, respectively, while the ∆t subscript refers to the timestep. The number of nns used for D i-1,k in nnGParareal and RandNet-Parareal are m nnGP = 20 and m RandNet = 4, respectively. N is the total number of intervals.discretizations, the number of intervals N and the order of the solvers F and G is given in Table 3 . Figure 6 highlights the strong scaling advantages of RandNet-Parareal compared to nnGParareal, setting N = 512 and restricting the runtime budget to a maximum of 48 hours, as done in the other test cases. Figure 6 : Scalability study for the 2 and 3 spatial dimensional Brusselator PDE. We used N = 512 and 48 hours runtime budget. nnGParareal for log 10 (d) = 3.9 is estimated, as the algorithm does not converge within the 48 hours runtime budget.In Table 4 below, we report the accuracies and runtimes (shown in parentheses) for RandNet-Parareal, Parareal, and nnGParareal. The accuracy is measured with maximum absolute error (mean across intervals) with respect to the true solution obtained by running F sequentially. Interestingly, all accuracies are far below the pre-defined accuracy level ϵ, with RandNet-Parareal achieving the lowest one in all but one experiment, with much smaller runtimes across all case studies.This section summarizes the simulation setups used for producing the results discussed in Section 5 in the main text. The tables below report the space and time domain of the considered PDEs, the number of spatial discretization points N x (and N y , in case of two-dimensional spatial systems), the numerical solvers used for G and F , their corresponding numbers of time steps per interval, the number of intervals N , and the number of nns used for nnGParareal (m nnGP ) and RandNet-Parareal (m RandNet ). In particular, Table 5 refers to the viscous Burgers' equation, Table 6 to the Diffusion-Reaction equation, and Table 7 to the shallow water equations (SWEs). 7.43e -9 (9h 50m) not converged SWE d = 3.1e 4 6.75e -8 (4h 9m) 5.15e -8 (15h 43m) not converged SWE d = 6.1e 4 8.54e -9 (12h 34m) 2.84e -8 (19h 30m) not converged Brusselator 2D d = 2e 3 2.09e -8 (2m) 3.16e -8 (7m) 3.38e -7 (7h 31m)Accuracy and computational cost comparison of RandNet-Parareal, Parareal, and nnGParareal for different PDEs, with runtimes reported in parentheses. The accuracy is measured as maximum absolute error (mean across intervals) with respect to F run sequentially.Table 5 : Simulation setup for the viscous Burgers' equationDomain Nx d G N G F N F N mnnGP m RandNet (x, t) ∈ [-1, 1] × [0, 5.9] 128 128 RK1 4 RK8 4e 4 128 18 3 (x, t) ∈ [-1, 1] × [0, 5.9] 1128 1128 RK1 293 RK8 6e 5 128 18 3N x is the number of space discretizations, the same as d here. G and F denote the chosen coarse and fine solvers, with corresponding time discretization steps per interval N G and N F , respectively.Here N is the number of intervals, while m nnGP and m RandNet are the numbers of nns used to create D i-1,k for nnGParareal and RandNet-Parareal, respectively.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, y) ∈ [-5, 5] × [0, 5], which we evolve temporally over t ∈ [0, 20].Following [75] , as an initial condition, we place a Gaussian-shaped column of water centered at (x, y) = (-2.5, 1.5), with covariance matrix Σ = 0.25 0 0 0.25 . We use Neumann boundary conditions, and evolve the system using N = 235 intervals over four increasingly finer spatial meshes, as described in Table 7 . We used the ParareaML [8] Python package to implement the SWEs and corresponding numerical solvers. 1, 2 and Figure 1 in the main text, and Figure 6 and Table 8 in Supplementary Material I. A 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 be 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 to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Justification: We have reviewed the Code of Ethics and found no 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 a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
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 1, and in Section 6 (when referring to 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, RandNet-Parareal can be applied to a wide range of applications, such as plasma physics simulation, weather forecasting (both mentioned in the introduction), leading to 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 societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or 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 a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is 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 of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: We relied on publicly available models, simulating the relevant data as described in the main text and in Supplementary Material.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] 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 use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: The code, simulations and associated analyses are publicly released with permissive licence. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of 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 is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS 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.F and G can be two different solvers or the same solver with different time steps. 
 
 Appendix:  


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: Paper Structure:
The authors provide experimental details within the introduction (lines 85-90) that should belong to the experimental section.
The 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.

To 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.
questions: 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?  
2. What is the “Fine” algorithm in Table 1 and in Fig. 1? I could not find its definition.
3. How can the Runtime in Figure 1 be negative? Is log(Runtime) shown in this plot?
4. 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? 
If the solution quality can be quantified, how does the solution quality of your method compare to other methods?
5. 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: limitations: The authors addressed some limitations of their work.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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).
2) 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).
3) 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.
4) 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).

A few minor issues:
1) I would avoid using "nns" for nearest neighbors (l.35), because it is confusing with "NNs" (neural networks).
2) 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.
questions: 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?

2) 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.

3) 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.

4) How does the high ambient-space dimension d of the state U affect the nearest neighbor search?
limitations: 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).
paper_summary: paper_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.
strength_weakness: strength_weakness: 1. The authors have not mentioned the trade-off between accuracy and computational cost, which is essential for high-fidelity simulations.

2. State-of-the-art comparisons are lacking, restricting the evaluation of the proposed method with random neural network-based PDE solvers.

3. 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.

4. 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.
questions: 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.

2. 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.

3. 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? 

4. Can the authors comment on the performance of their method when using a deep network instead of a shallow network?
limitations: 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: - The studied problem is domain-specific, which might not be of general interest to the scientific machine-learning community. 

- 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. 

- 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.
questions: questions: In lines 7-8, is it a character "x" or a symbol "\times" for `x125`? Probably a symbol is preferred.
limitations: limitations: N/A.

Comments: title: Raising my score
Last day for discussion
Follow-up answer to Reviewer FREq's comments and raised score
Raised Score
Follow-up to Reviewer Eu5i's comment
Follow-up to Reviewer awZW's comment
Response for rebuttal
Follow-up answer to Reviewer FREq's comments
Answer to Rebuttal
comment: Thank you so much for providing detailed answers. Considering all the answers, I have raised my score from 5 to 6. It would be nice to include trade off results in the revised paper.
Reviewer nFhC, today is the last day for discussion. I hope you can take a moment to respond to the authors' rebuttal.
We are glad that our rebuttal comments were helpful in addressing your concerns. We are also grateful for reflecting this in your score, and for your further suggestions.

We notice that there are still some lingering misconceptions about our work, which we believe are important to point out. 

Our work does not "present an improvement to a numerical ODE solver". Instead, it improves existing _Parallel-in-Time_ (PinT) methods for solving ODEs and PDEs, as written in the paper abstract, summarized by the other reviewers, and discussed at length during the rebuttal. We do, ultimately, solve an (O/P)DE with a numerical scheme, but our proposed scheme  
(1) solves DEs by parallelizing in time, belonging thus to a different class than that of the "simple" numerical ODE solvers;  
(2) models the correction term using particularly suited **random neural networks**, allowing us to achieve significant speed-up compared to existing PinT techniques;  
(3) can incorporate existing (O/P)DE solvers in its architecture;  
(4) is not meant to improve a specific ODE solver, but instead to make PinT techniques more effective and, hence, more appealing to practitioners in the field.

Moreover, our new approach (even if considered to be an improvement to an ODE solver, an oversimplification which we disagree with) is far from "solely relevant to obtaining training data for PDE machine learning methods". Numerical methods for solving DEs have profoundly impacted scientific progress since the last century, extending far beyond the generation of training data generation. In fact, once again, **our method can incorporate existing PDE machine learning solvers instead of training data for them.** Independently on that, landmark achievements such as weather forecasting, space exploration, and advancements in nuclear energy are just a few examples of the broad and transformative influence of numerical DE solvers.

Commenting on the relevance of designing efficient PinT solvers beyond AI, the following non-exhaustive list contains some other areas of natural sciences, social sciences, and engineering where the availability of efficient (O/P)DE solvers could be a game-changer leading to extremely important advances:

- *Atmospheric Dynamics:* Accurate solutions can mitigate climate change, enhance disaster preparedness, and optimize agriculture.
- *Biomedical Processes:* Real-time simulations of biological processes, such as tumor growth and blood flow, can aid in personalized treatment planning and targeted interventions.
- *Earthquake Engineering:* Refined modeling of seismic wave propagation and structural responses enhances earthquake predictions and leads to more resilient building designs.
- *Financial Modeling:* Improved pricing accuracy and real-time risk assessment can stabilize financial markets and reduce systemic risk.

We believe our work is another step toward providing faster numerical solutions for DEs in any field where this may be required, whether in the natural sciences, engineering, machine learning, or artificial intelligence.

We thank you for raising the score and hope that our further comments help convey the main scope of our contributions.
I appreciate the authors' efforts in addressing the concerns raised in the rebuttal. Based on this, I am increasing my score from 4 to 5.

The paper presents an improvement to a numerical ODE solver, which is solely relevant to obtaining training data for PDE machine learning methods.

To enhance the paper's relevance to the AI community, I suggest two improvements:

1.  Provide a GPU implementation of the proposed method.
2.  Include a comparison between this method and other ODE solvers commonly used in simulating the reverse process of diffusion models.

These additions would help demonstrate the practical applications and potential advantages of the method in AI-related contexts.
We are glad you appreciated our additional theoretical and numerical results. We thank you again for your constructive criticisms, we believe the revised paper will greatly benefit from them. We are also extremely grateful that you have upgraded your score notably.
Thanks for your comment on our rebuttal and W2.

We would appreciate if you could read what we wrote as a comment to Reviewer FREq about the impact of our paper on AI, which we report below for your convenience. We hope tis may reduce some of your reservations on the significance of our work. Properly communicating the relevance and impact of our method is very important for us. 

Our work exemplifies an instance when ML methods allow the advancement of various fields of science and engineering where solving (O/P)DEs is needed, see also our response to Ref FREq to his further question. We mentioned these applications in our initial rebuttal as our submission belongs to the Primary Area of *ML for Physical Sciences*.

However, AI is also one of the fields that could greatly benefit from progresses regarding efficient ODE solvers, and thus our method. This is why we gladly provide examples of how our approach can also assist ML and AI, adding these points in the camera-ready version of the paper.

ODEs are a crucial building block of some relevant techniques, such as *Diffusion models (DMs)*, *Neural ODEs*, *Optimal control and reinforcement learning*, and *Optimization for ML models*, to name a few. Below, we provide an example of how our work could lead to substantial advancements for DMs.

- _Solvers for DMs_: In the context of DMs, the continuous time reversal process can be described by a probability flow ODE, defined by the score function, usually approximated with deep (convolutional) NNs (note that the ODE formulation offers significant advantages over SDEs in high dimensions). Faster ODE solvers can improve sampling speed, yielding faster image synthesis. The main existing directions in the literature involve using classical sequential solvers [30], developing faster and dedicated ones such as DDIM [24], DDPM [29], DPM-Solver [25], Heun [26], and parallelization of the autoregressive sampling process of DMs (by using Picard-Lindelöf iterations for ODE/SDE) [27,28]. The latter emerged mainly due to the need for other solutions to avoid the bottlenecks of sequential solvers. To the best of our knowledge, modulo this particular parallelized sampling, PinT schemes have not yet been used in this context. **This is the gap this paper could fill, as the existing successful sequential dedicated solvers could be embedded into our proposed RandNet-Parareal**. We expect this to be straightforward, since the goal is to collect samples via solving the corresponding (diffusion) ODEs at $[0,T]$, where RandNet-Parareal immediately finds its purpose.

We refer to our latest comment to Reviewer FREq for an additional discussion on the GPU implementation which could further improve the efficiency of our method.

We hope these clarifications provide you with evidence of the expected impact/implications of our method in the context of generative AI models and beyond.

**References used above:**  
[24] J. Song, C. Meng, S. Ermon. Denoising diffusion implicit models. 2020.  
[25] C. Lu, Y. Zhou, F. Bao, J. Chen, C. Li, J. Zhu. DPM-solver: A fast ODE solver for diffusion probabilistic model sampling in around 10 steps. 2022.  
[26] T. Karras, M. Aittala, T. Aila, S. Laine. Elucidating the design space of diffusion-based generative models. 2022.  
[27] A. Shih, S. Belkhale, S. Ermon, D. Sadigh, N. Anari. Parallel sampling of diffusion models. 2023.  
[28] Z. Tang, J. Tang, H. Luo, F. Wang, T.-H. Chang. Accelerating parallel sampling of diffusion models. 2024.  
[29] J. Ho, A. Jain, P. Abbeel. Denoising diffusion probabilistic models. 2020.  
[30] Y. Song, J. Sohl-Dickstein, D.P. Kingma, A. Kumar, S. Ermon, B. Poole. Score-based generative modeling through stochastic differential equations. 2021.
Thanks for your rebuttal, especially the response to **W2**.  I still have some reservations about the significance of the work. I will be maintaining my score.
I highly appreciate the additional results, both experimental and theoretical. I very much like the general idea of the paper but was unsure about the soundness in these two aspects, both of which were addressed. I will raise my score to 7.
* "Should we add [the new experiments] to Suppl. Material?" Yes, please.
**AQ1-3**: Thank you for your feedback. We will carefully take it into account when preparing the final version. 

We welcome any comment on our response to **Q4** on the solution quality of the solver.  

**AQ5: Dropping indices without prior mention.** 
Please note that we specify that on Line 195 in [paper], which reads "Prior to that [i.e., defining RandNet for Parareal], we define how RandNets work in a general setting with input $\boldsymbol{U}$". Indices are introduced later, upon combining RandNets with Parareal, in line 237. Should you find this still unclear to the reader, we could change it to "Prior to that, we define how RandNets work in a general setting with input $\boldsymbol{U}$, *before going back to the input of interest $\boldsymbol{U}_i^k$ within the Parareal framework*".

**AW2 and Further Q: On the potential impact of RandNet-Parareal on ML/AI, and, specifically, on Diffusion Models (DMs).**
Thank you for engaging in further conversation about this matter. Properly communicating the relevance and impact of our method is very important for us. 

Our work exemplifies an instance when ML methods allow the advancement of various fields of science and engineering where solving (O/P)DEs is needed, see also our response to W2 above, and W1 of Ref awZW. We mentioned these applications as our submission belongs to the Primary Area of *ML for Physical Sciences*.

However, AI is also one of the fields that could greatly benefit from progresses regarding efficient ODE solvers, and thus our method. This is why we gladly provide examples of how our approach can also assist ML and AI, adding these points in the camera-ready version of the paper.

ODEs are a crucial building block of some relevant techniques, such as *DMs*, *Neural ODEs*, *Optimal control and reinforcement learning*, and *Optimization for ML models*, to name a few.

- _Solvers for DMs_: In the context of DMs, mentioned in your answer, the continuous time reversal process can be described by a probability flow ODE, defined by the score function, usually approximated with deep (convolutional) NNs (note that the ODE formulation offers significant advantages over SDEs in high dimensions). Faster ODE solvers can improve sampling speed, yielding faster image synthesis. The main existing directions in the literature involve using classical sequential solvers [30], developing faster and dedicated ones such as DDIM [24], DDPM [29], DPM-Solver [25], Heun [26], and parallelization of the autoregressive sampling process of DMs (by using Picard-Lindelöf iterations for ODE/SDE) [27,28]. The latter emerged mainly due to the need for other solutions to avoid the bottlenecks of sequential solvers. To the best of our knowledge, modulo this particular parallelized sampling, PinT schemes have not yet been used in this context. **This is the gap this paper could fill, as the existing successful sequential dedicated solvers could be embedded into our proposed RandNet-Parareal**. We expect this to be straightforward, since the goal is to collect samples via solving the corresponding (diffusion) ODEs at $[0,T]$, where RandNet-Parareal immediately finds its purpose.

- _GPU Implementation_: There are two increasingly sophisticated approaches for implementing RandNets-Parareal on GPUs: (i) parallel implementation of the solvers (model parallelism); (ii) parallel implementation of Parareal (in-time parallelism). For (i), both the solver and the RandNet computation can be carried out on the GPU, see e.g.[31]. In-time parallelism would then be implemented across GPUs. Otherwise, under (ii), both model and time parallelism can be implemented concurrently on the same hardware [32]. The latter guarantees a better use of resources at the cost of increased implementation complexity.  

We hope these clarifications provide you with evidence of the expected impact/implications of our method in the context of generative AI models and beyond.

**References used above:**  
[24] J. Song, C. Meng, S. Ermon. Denoising diffusion implicit models. 2020.  
[25] C. Lu, Y. Zhou, F. Bao, J. Chen, C. Li, J. Zhu. DPM-solver: A fast ODE solver for diffusion probabilistic model sampling in around 10 steps. 2022.  
[26] T. Karras, M. Aittala, T. Aila, S. Laine. Elucidating the design space of diffusion-based generative models. 2022.  
[27] A. Shih, S. Belkhale, S. Ermon, D. Sadigh, N. Anari. Parallel sampling of diffusion models. 2023.  
[28] Z. Tang, J. Tang, H. Luo, F. Wang, T.-H. Chang. Accelerating parallel sampling of diffusion models. 2024.  
[29] J. Ho, A. Jain, P. Abbeel. Denoising diffusion probabilistic models. 2020.  
[30] Y. Song, J. Sohl-Dickstein, D.P. Kingma, A. Kumar, S. Ermon, B. Poole. Score-based generative modeling through stochastic differential equations. 2021.  
[31] A.N.Budko, M. Möller, C. W. J. Lemmens. Time integration parallel in time. 2020.   
[32] A.Q. Ibrahim, S. Götschel, D. Ruprecht. Parareal with a physics-informed neural network as coarse propagator.
**AQ1:**
I appreciate the authors clarification on the background and encourage the authors to add a corresponding section into the paper.

**AQ2:**
I think this should be defined more clearly in the experiment section.

**AQ3:** The y-axis label is incorrect, as it is denoted in hours but should be log10(hours). It would be better to scale the y-axis logarithmically while keeping the y-axis ticks in hours, as this is easier to interpret than log10(hours).

**AQ4:** The authors should not drop indices without mentioning it, and they should add a footnote for clarification.

**AW2:** While this may not be a NeurIPS requirement, the reviewers should still consider the potential impact of the method in AI when assigning ratings. The reviewers noted the method's potential applicability to ODE solvers used in machine learning, such as in diffusion models. The paper would benefit from a section discussing this potential application in more detail.

**Further Question:** Are PinT schemes used to solve the reverse ODE in Diffusion Models? If not, which ODE solvers are used? What would be necessary to apply the authors' method in this setting? Likely it would need to be implemented on a GPU?
 Demonstrating the benefit of the authors' method in the context of diffusion models would significantly increase the relevance of this paper to the NeurIPS community.