A Stable and Scalable Method for Solving Initial Value PDEs with Neural Networks
Authors: Marc Anton Finzi
Andres Potapczynski
Matthew Choptuik
Andrew Gordon Wilson
Venue: ICLR 2023 poster (2023)
Keywords: PDE
Neural PDE solvers
Initial value problems

Abstract: Unlike conventional grid and mesh based methods for solving partial differential equations (PDEs), neural networks have the potential to break the curse of dimensionality, providing approximate solutions to problems where using classical solvers is difficult or impossible. While global minimization of the PDE residual over the network parameters works well for boundary value problems, catastrophic forgetting impairs applicability to initial value problems (IVPs). In an alternative local-in-time approach, the optimization problem can be converted into an ordinary differential equation (ODE) on the network parameters and the solution propagated forward in time; however, we demonstrate that current methods based on this approach suffer from two key issues. First, following the ODE produces an uncontrolled growth in the conditioning of the problem, ultimately leading to unacceptably large numerical errors. Second, as the ODE methods scale cubically with the number of model parameters, they are restricted to small neural networks, significantly limiting their ability to represent intricate PDE initial conditions and solutions. Building on these insights, we develop Neural-IVP, an ODE based IVP solver which prevents the network from getting ill-conditioned and runs in time linear in the number of parameters, enabling us to evolve the dynamics of challenging PDEs with neural networks.

Published as a conference paper at ICLR 2023
A STABLE AND SCALABLE METHOD FOR SOLVING INITIAL VALUE PDES WITH NEURAL NETWORKS
Mare Finzi' Andres Potapczynski', Matthew Choptuik?, Andrew Gordon Wilson'
New York University' and University of British Columbia²
ABSTRACT
Unlike conventional grid and mesh based methods for solving partial differen- tial equations (PDEs), neural networks have the potential to break the curse dimensionality, providing approximate solutions to problems where using clas- sical solvers difficult or impossible. While global minimization of the PDE residual over the network parameters works well for boundary value problems, catastrophic forgetting impairs applicability to initial value problems (IVPs). In an alternative local-in -en--ii-teee approach, the optimization problem can be converted into an ordinary differential equation (ODE) on the network parameters and the solution propagated forward in time; however, we demonstrate that current meth- ods based on this approach suffer from two key issues. First, following the ODE produces an uncontrolled growth the conditioning of the problem, ultimately leading to unacceptably large numerical errors. Second, as the ODE methods scale cubically with the number model parameters, they are restricted to small neural networks, significantly limiting their ability to represent intricate PDE ini- tial conditions and solutions. Building on these insights, we develop Neural-IVP, ODE based IVP solver which prevents the network from getting ill-conditioned and runs in time linear in the number of parameters, enabling us to evolve the dy- namics of challenging PDEs with neural networks.
INTRODUCTION
Partial differential equations (PDEs) are needed to describe many phenomena in the natural sciences. PDEs that model complex phenomena cannot be solved analytically and many numerical techniques are used to computer their solutions. Classical techniques such as finite differences rely on grids and provide efficient and accurate solutions when the dimensionality is low (d 2). Yet, the computational and memory costs of using grids or meshes scales exponentially with the dimension, making it extremely challenging to solve PDEs accurately in more than 3 dimensions.
Neural networks have shown considerable success in modeling and reconstructing functions on high- dimensional structured data such images or text, but also for unstructured tabular data and spatial functions. Neural networks sidestep the "curse of dimensionality" by learning representations of the data that enables them perform efficiently. In this respect, neural networks have similar benefits and drawbacks as Monte Carlo methods. The approximation error e converges at a rate x 1/ v/vn from statistical fluctuations where is the number of data points or Monte Carlo samples. Expressed inversely. would need: x e2 1og 1/e samples to get error c, a compute grows exponentially in the number of significant digits instead of exponential in the dimension as it is for grids. For many problems this tradeoff is favorable and an approximate solution is much better than no solution.
Thus, it is natural to consider neural networks for solving PDEs whose dimensionality makes stan- dard approaches intractable. While first investigated in Dissanayake & Phan-Thien (1994) and La- garis et al. (1998), recent developments by Yu et al. (2018) and Sirignano & Spiliopoulas (2018) have shown that neural networks can successfully approximate the solution by forcing them sat- isfy the dynamics of the PDE on collocation points in the spatio-temporal domain. In particular, the global collocation approaches have proven effective for solving boundary value problems where the neural network can successfully approximate the solution. However, for initial value problems
"qqull contribution, order chosen by random coin flip. {ma 1820, ap6604 ap6604)enyu.es.

---

Published as a conference paper at ICLR 2023
(IVPs), treating time merely another spatial dimension results in complications for the neural network like catastrophic forgetting. Some heuristics have been developed to ameliorate this latter problem, such increasing the collocation points time progresses, but then the computational cost of training the neural network becomes impractical.
Recently, Du & Zaki 2021) and Bruna et al. (2022) have provided two methods that follow novel local-in-tine time approach for training neural networks to solve IVPs by updating the network param- eters sequentially through time rather than by having some fixed of parameters model the whole spatio-temporal domain. These methods have proven successful for a variety of PDEs, but they currently suffer from two shortcomings. First, the conditioning of the linear systems required to follow the ODE on the network parameters degrades over time, leading to longer solving times and ultimately to complete breakdown of the solution. Second, the current methodologies lack the capacity to represent difficult initial conditions and solutions their runtime scales cubically in the number of network parameters, limiting their ability use large neural networks. In this work we provide a local-in-time IVP solver (Neural-IVP) that circumvents the shortcomings of Du & Zaki 202T
and Bruna et al. (2022) and thus enable us to solve challenging PDEs. In particular:
Leveraging fast matrix vector multiplies and preconditioned conjugate gradients, we develop an approach that scales only linearly in the number of parameters, allowing us to use considerably larger neural networks and more data.
We further improve the representational power and quality the fit to initial conditions through the use of last layer linear solves and sinusoidal embeddings.
We show how following the parameter ODE leads the network parameters to an increasingly poorly conditioned region of the parameter space, and we show how this relates exact and approximate parameter symmetries in the network.
Using regularization, restarts, and last layer finetuning, we are able to prevent the parameters from reaching these poorly conditioned regions, thereby stabilizing the method.
We provide a code implementation at attpps:/ggtaa...m///ttta///eeeee
BACKGROUND
Given spatial domain X C RDD, we will consider the evolution of time-dependent function : X' x [0, [7] + R which at all times belongs to some functional space U and with dynamics governed by
d.u (x, t) L [u] (x,
for (x,t) X' [0,
u (I. 0) = 14 (x) (I, = h (x,t)
for x E
for 0X x [0,
where uo 2/ is the initial condition, h is spatial boundary condition, and L is the (possibly non- linear) operator containing spatial derivatives. We can represent PDEs with higher order derivatives in time, such the wave equation 0 = Ad. by reducing them system of first order in time equations == [ó, 0 where in this example L[uo, u1] [uI. Auo].
Global Collocation Methods The first approaches for solving PDEs via neural networks are based on the idea of sampling uniformly on the whole spatio-temporal domain and ensuring that the neural network obeys the PDE by minimizing the PDE residual (or a proxy for it). This approach (194)) and Lagaris et al. (1998). which used was initially proposed by Dissanayake & However, recent advances automatic differentiation, neural networks as approximate solutions. compute, and neural network architecture have enabled successful applications such as the Deep Galerkin Method (Sirignano & Spiliopoulos) 2018), Deep Ritz Method (Yu al. 2018). and PINN Raissi et al. 2019), which have revitalized interest using neural networks to solve PDEs.
Learning From Simulations Not all approaches use neural networks as basis function to repre- sent the PDE solution. Some approaches focus on directly learning the PDE operator as in Lu et al. (2019) or Kovachki et al. (2021). where the operator can be learned from simulation. However, as these methods typically use grids, their purpose is to accelerate existing solvers rather than tackling new problems. Other approaches that do not rely on collocation points exploit specific information of elliptic and semi-linear parabolic PDEs, such as E. et al. (017) and Han et al. (2018).
2

---

Published as a conference paper at ICLR 2023
2.1 GLOBAL PDE RESIDUAL MINIMIZATION
The most straightforward method for producing neural network solution to initial values PDEs similar to the approach used for boundary value problems of treating the temporal dimensions as if they were spatial dimensions and parameterizing the solution simultancously for all times u(x,t) = No(x, t). The initial and boundary conditions can he enforced through appropriate pa- rameterization of the network architecture (Berg & Nyström 2018), whereas the PDE is enforced through minimization the training objective:
/ xx(TT.T)
ro(a, (((x)))d J Jxx1o.T1 [(ruul( () L [uo] (z, ())))d(()))
S(0)
where the expectation estimated via Monte Carlo samples over a chosen distribution of μu and times t.
Initial value PDEs have a local temporal structure where only values on the previous spatial slice are necessary to compute the next; however, global minimization ignores this property. Moreover, as the weights of the neural network must be used to represent the solution simultancously at all times, then must ensure that the neural network approximation does not forget the PDE solution learnt at earlier times (catastrophic forgetting). While SSriinnano Spiliopoulos/ (2018) and [Sizzmann) et al. (2020) take this approach, the downsides of avoiding catastrophic forgetting is to increase the computation spent by ensuring the presence of data from previous times.
2.2 LOCAL-IN-TIME METHODS
To circumvent the inherent inefficiency of the global methods, Du & Zaki (021) and Bruna et al. (0222) propose a local-in-time method whereby the minimization problem gets converted into an ODE that the parameters satisfy at each point in time. In this approach, the PDE solution is given by u(x,t) (ax a(t)) for a neural network N, where the time dependence comes from the parameter vector (()) rather than as an input to the network. Thus, the network only represents the solution at single time, rather than simultaneously at all times, and subsequently O(t) can be recorded and no representational power or computational cost is incurred from preserving previous solutions. Assuming the PDE is one-dimensional, the PDE residual at single time can be written as,
L(,,t) / r(a, ddxt)))(((=== / (0 (ØT VoN(a, ((t) ([)
(1)
since the time derivative is 0,u(z, = 0T (eN(x,
Choosing the dynamics ê of the parameters minimize the instantaneous PDE residual error L(ê,t) yields the (implicitly defined) differential equation
argmin min / (N(x, - u(x, 0))
M = F(0) and 0o = arg
(2)
where M()) J Vz, 0) du(x) and F(0) Jx VoN(a, ()] l)dp(x).
Once we find 6g to fit the initial conditions, we have a fully specified system differential equations, where we can advance the parameters (and therefore the solution u(x; o()))) forward in time.
Since both (0) and F(0) involve integrals over space, we can estimate them with n Monte Carlo samples, yielding (0) and F(B). We then proceed to solve the linear system M (0) = F(0) at each timestep for the dynamics 0 and feed that into ODE integrator such as RK45 (Dormand & Prince 1980). For systems of PDEs such as the Navier-Stoker equations, the method can be extended in straightforward manner by replacing the outer product of gradients with the Jacobians of the multi-output network N: (0) = fx D2N(x, 0) DaN(x, 0)d)(=) and likewise for F, which results from minimizing the norm of the PDE residual fx [|r( 1) (ddj(x).
Introducing some additional notation, illuminating way. Defining the Jacobian matrix the network for different input points Ju = 39,
write the Monte Carlo estimates M and F in more
0 N(v;, 0), and defining f
(2,, 0), the PDE residual estimated via the n-sample
3

---

Published as a conference paper at ICLR 2023
Monte Carlo estimator is just the least squares objective t(),t) [| - fl?.. The matrices M (0) = _J J and F(0) = J reveal that the ODE dynamics is just the familiar least squares solution = = = (J (J
DIAGNOSING LOCAL-IN-TING NEURAL PDE SOLVERS
The success local-in-tines methods hinges on making the PDE residual L(ô, t) close to 0 as we follow the dynamics of 0 = M (0) (10)-1fpp (()) The lower the local error, the lower the global PDE residual S(Ø) = f L(B, t) and the more faithfully the PDE is satisfied.
Even though á directly minimizes [(0, t). the PDE residual is not necessarily small and instead the value of r(x, t) depends nontrivially on the network architecture and the values of the parameters themselves. While local-in-time methods have been applied successfully in several cases, there are harder problems where they fail unexpectedly. For example, they fail with unacceptably large errors in second-order PDEs or problems with complex initial conditions. In the following section, we identify the reasons for these failures.
3.1 REPRESENTATIONAL POWER
The simplest reason why enough representational power. Having enough degrees of freedom and inductive biases in the network matters for being able to find a the span of J which can match the spatial derivatives. The spatial derivatives (1,, 0) of the PDE must be able to be expressed (or nearly expressed) a linear combination of the derivative with respect to each parameter: 0 Nx, 0), which is 89, different task than neural network is typically designed for. The easiest intervention is to simply increase the number of parameters p. yielding additional degrees of freedom.
methods fail is because the neural networks do not have
Increasing the number of parameters also improves the ability of the network reconstruct the initial conditions, which can have knock-on effects with the evolution later in time. However, increasing the number of parameters and evolving through time following Du & Zaki (202T) and Bruna et al. (2022) quickly leads to intractable computations. The linear solves used to define the ODE dynam- ies require time Opp" + p²n) and use Opp² + pn) memory, where p represents the neural network parameters and n. is the number of Monte Carlo samples used to estimate the linear system. There- fore the networks with more than around p = 000 parameters cannot be used. Neural networks of these sizes are extremely small compared to modern networks which often have millions or even billions of parameters. section 4.1 we show how our Neural-IVP method resolves this limitation, allowing us to use large neural networks with many parameters.
3.2 STABILITY AND CONDITIONING
In addition lacking sufficient representational power, there are more subtle reasons why the local- in-time methods fail.
Even when the solution is exactly representable, a continuous path 0* (t)) between the solutions may not exist. Even if a network is able to faithfully express the solution at a given time u(I, t) N [2, 0*) for some value of 0* in the parameter space, there may not exist continuous path between these 6* for different times. This fact is related to the implicit function theorem. With the multi- output function (() N Nz 0) u(x,, t). even if we wish satisfy H; 0 only at a finite collection of points 2i, the existence of a continuous path 0" (4) = g(() in general requires that the Jacobian matrix D, H-
invertible. Unfortunately the Jacobian is not invertible because there exist singular directions and nearly singular directions in the parameter space,
now argue.
There exist singular directions of J and M as a result of symmetries in the network. Each continuous symmetry of the network will produce a right singular vector of J, regardless of how many points n are used in the Monte Carlo estimate. Here we define a continuous symmetry as parameterized transformation of the parameters Ta R D" defined for neighborhood of the identity To = Id, and T, has a nonzero derivative with respect to a at the identity. For convenience, consider reparametrizing be unit speed so that | (!)] 1.
c), in
4

---

Published as a conference paper at ICLR 2023
10
Adveciion
102
Have
Advection Wave (SGD)
Nave Adveciion Viav (5G0)
10
Espervalies
10
Max Elgenveluv
Iterators
110
30->
N000
10a
10-1
200
400
600
800
2
1.0
Figure 1: The conditioning of the linear systems needed solve the ODE on the network param- eters increases for challenging PDEs like the wave equation, but not for others like the advection equation. (Left): Eigenvalue spectrum of M(0) matrix initialization. (Middle): Growth of largest eigenvalue of (0)) over time. (Right): Number of preconditioned CG iterations required to solve the linear system a specified tolerance of 10-?
Theorem 1. Suppose the network N(x, 0) has continuous parameter symmettrrrymmmeerry Ta which preserves the oulputs of the function: v0, x : Næx Ta (0)) = [Tx, 0), then
v(0) (a) a=
is singular vector both J and M.
Proof: Taking the derivative with respect to a at O, from the chain rule we have: 0. (ax Ta (6)) VoN(x,, (aae,,)TTTTaaa (aTa())) As this expression holds for all x, J(0)u(0) and (()vv(()) 0.
As Dinh et al. (2017) demonstrated, multilayer perceptrons using ReLU nonlinearities have a high- dimensional group of exact parameter symmetries corresponding to rescaling of weights in alter- nate layers. Furthermore even replacing ReLUs with alternate activation functions such as Swish (Ramachandran al. (017) does solve the problem, these will have approximate symmetries which will produce highly ill colcinnit M and J matrices.
Theorem 2. An approximate symmetry Vx [| Ta (0) 0x10))11= Ea2 will produce nearly singular vectors v()) = T. Ta (0) (0) for which
MM <,
and therefore the smallest eigenvalue M is less than e.
Proof: See Appendix A
Additionally, the rank of the Monte Carlo estimate MT = there is p dimensional manifold of parameters which match the function values the sample points = =,, 1,..., n N(r, 0) = N, rather than over the whole domain. In Figure (left), we show empirically that the eigenspectrum of M is indeed deficient and highly ill ill-condiiooned, with a long tail of small eigenvalues.
2. LJJJ using n samples is at most n, and
Hence some form of regularization such as [M (0) + pT] ê necessary have a bounded condition number, & (M (0) pu) = (A1 + ()) (0 (t).
Furthermore, as seen in Figure (middle) the conditioning of the linear system (equation 2 dete- riorates over time. This deterioration worsens in more challenging PDEs like second-order wave equation contrast to the easier advection equation. Even when using dense solver, the quality of the solves will degrade over time, Icading to increased error in the solution. When using iterative method like CG shown in Figure I (right), the runtime of the method will increase during the evo- lution and eventually not be able to meet the desired error tolerance. In contrast, if we instead take snapshots of the solution at different times and fit it directly with SGD, we find that the conditioning much better as shown by the green curve in Figure
Making sense of this observation, we can learn from the ways that neural networks are typically used: in conjunction with stochastic gradient descent. In general when training neural networks with
5

---

Published as a conference paper at ICLR 2023
SGD, must chose carefully the initialization to make the problem well-conditioned (Mishkin & Matas 2015). Many initializations, such as setting the scale of the parameters with drastically dif- ferent values between layers, will lead either to diverging solutions or no progress the objective. With the right initialization and a good choice of hyperparameters, the optimization trajectory will stay in well-conditioned region of the parameter space. However, many bad regions of the pa- rameter space exist, and number of architectural improvements in deep learning such as batch normalization (loffe & Szegedy) 2015) and skip connections (He al. (016) were designed with the express purpose of improving the conditioning of optimization while leaving the expressive power unchanged. Unfortunately, while SGD optimizers stay in a well-conditioned regions of the parameter space, equation (2) does not.
To see this, consider a singular vector v of which has a very small singular value oa (due to an approximate symmetry or otherwise). Analyzing the solution, we see that the projection of the parameters along the singular vector evolves like: v v = a J o, f. where Jr 0 A small singular value leads to large change in that subspace according to the evolution of the ODE, Considering the approximate rescaling symmetry for swish networks, this points the dynamics directly into amplifying the difference between size of weights in neighboring layers and worsening the conditioning, precisely what was identified sharp minima from Dinh et al. (201)). We describe how this problem can be circumvented by our method as described section 4.2
4 NEURAL IVP
Drawing on these observations, we introduce Neural IVP, a method for solving initial value PDEs that resolves the scalability and numerical stability issues limiting current local-in-time methods. We also enhance the Neural IVP, through projections to different regions of the parameter space, fine- tuning procedures, and mechanisms increase the representation power of the neural networks.
4.1 IMPROVING REPRESENTATIONAL POWER AND SCALABILITY
To evaluate the representational power of the network, we examine fitting complex initial condi- tion typically present in the later evolution of an entangled system, and which will contain com- ponents at many different frequencies. We fit sum of two Gaussians of different sizes mod- ulated by frequencies pointing in different directions in 3-dimensions within the cube 1|s, with the slice through 2 = shown in Figure 2 (left). The function is defined as the sum of two Gaussian like wave packets, modulated by spatial frequencies pointing in different directions: up (1)) 30 (2 (lv1111/((2e)) cos (2t fx T f) 24 (2T e -1w (1 cos (2T T m)) where v = x 0.5(2 + its,, w = x x (01 + 82 + l///66 also î = (01 + 22) / //2 and
m 2-1 (i- 21 $3)).
By changing the frequency parameter f. we can investigate how well the network is able to fit fine-scale details. We introduce three substantial improvements
the models over those used in Evolutional Deep Neural Networks (HDNN) of Du & /aki (202T) in order to improve their repre- sentational power, and evaluate the impact of these changes in Figure 2 (middle), starting with the exact 4-layer 30 hidden unit tanh nonlincarity MLP architecture used in EDNN.
Increasing Number of Parameters and Scalability As shown in Figure 2 (right), the number of network parameters has a large impact on its representational power, not just for solving the initial conditions but also for finding a that achieves low PDE residual. Increasing the number of parameters substantially reduces the approximation error, especially at high frequencies which are challenging for the model. While dense solves prohibit scaling past 5,000 parameters, we show that our solves can be performed much faster making use of the structure of M. Matrix vector multiplies with the matrix M (0) IJ can be implemented much more efficiently than using the dense matrix. Making use of Jacobian- vector-produee implemented using automatic differentiation (such as in JAX (Bradbury et al. 2018)), can implement matrix vector multiply using 2 Jacobian- vector- productss
(5) which takes Onn + p) time and memory for a single matrix-vector product, sidestepping memory bottlenecks that are usually the limiting factor. Then, with these efficient matrix-vector products,
(0) e 2m Jel1?

---

Published as a conference paper at ICLR 2023
Relaeive Enor
10
Base EDMN
Sinuvaddal Embedding
10-2
Improved
Scaling Head furning
10
-2
202
102
Un(x) with f=
Frequency
Frequency
Figure 2: (Left): Example wave packet initial condition with varying levels of fine details param- eterized by the frequency f. (Middle): Impact of Neural-IVP improvements in the model on the initial condition fit relative error across different levels of difficulty of the solution as parameterized by the frequency f. yielding improvement of 1 orders of magnitude. (Right): Initial condition fit with all Neural-IVP interventions but with varying number of parameters in the model (shown by the colors). Note that the largest networks which can be used by the dense method are only 5000 parameters.
we can use a scalable Krylov subspace routine of conjugate gradients (CG). To further accelerate the method, we construct Nyström preconditioner (Frangella al., 2021 drastically reducing the number of CG iterations. Using this approach, each solve takes time (( + p) where x(P-IM) is the condition number of the preconditioned matrix M, rather than O(p² + p²n) time of dense solves. These improvements to the runtime using the structure of M mirror the sparse and Kronecker structures used by finite difference methods.
Sinusoidal Embedding We make several architectural enhancements to the networks used in Du & Zaki (2021) and Bruna et al.) (2022), which improve the quality of the initial condition fit and the error when evolving forward time. Notably, find that using the sinusoidal embed- ding (Mildenhall et al. 2021) substantially improves the ability of the network to represent higher frequency details in the solution. In contrast to the original form, we use the featurization
(6) which scales the magnitude of the high frequency (large w) components down by 1/wa. While = 0 (the original sinusoidal embedding) works the best for fitting an initial signal (the only requirement needed for Neural Radiance Fields (Mildenhall al. 202T)) the derivatives of will not be well behaved as the magnitude of the largest components of ? (z) will scale like 2 and )(xx) will scale like 22r. We find setting
(x)) [sin(2 x112=xx))) ) - IK=0 + [cos(2* x ak k
be the most effective for both first order and 2nd
order PDEs. Figure 2 (middle) shows the sinusoidal embedding helps the model represent complex functions.
Last Layer Linear Solves To further improve the quality of the initial condition fit, after training the network on the initial condition, we recast the fitting of the last layer of the network as solving linear least squares problem. The network up until the last layer can be considered as features N(z) = w de(x)) + b over fixed set of collocation points X. We can then solve the minimization problem with respect to the final layer weights w.
min I|w (X) + (x))? m
which can be solved to higher level of precision and achieve lower error when compared to the values of the last layer obtained without tuning from the full nonlinear and stochastic problem.
Combining these three improvements of scalability, sinusoidal embeddings, and last layer linear solves (head tuning), we are able to reduce the representation error of the networks by 1-2 orders of magnitude across different difficulties of this challenging dimensional problem.
4.2 STABILITY AND CONDITIONING
Preconditioning In section 3.2 we discussed how even for easier PDEs, the symmetries in the
neural networks generate badly conditioned linear systems for the ODE the parameters. To counteract this negative effect on our CG solver, we use the highly effective and scalable ran- domized Nyström preconditioner (Frangella et al. 2021) As discussed in Frangella al. (2021)

---

Published as a conference paper at ICLR 2023
10
Neurakiv
Neuralop
10-1
10-1
PDE residual
iters
PDE
10
10-4
000 025 0.2
1.00 3.75
2.0 0.2
0.4
0.5 03.8 1.0
0.0
0.2
0.4
3..6
0.5
1.0
Time
Figure 3: (Left): Restarts improve the condition of the linear systems. Here we cap the number of CG iterations at 1000, but without restarts the number required to reach a desired error tolerance will only continue to increase. Neural-IVP achieves order of magnitude better PDE residual than EDNN on the Fokker- Plank equation (middle) and on the Vlasov equation (right).
this preconditioner is close to the optimal truncated SVD preconditioner and it is empirically impressive. To use this preconditioner, we first construct a Nyström approximation of (0): Mays (0) (M (0) 52) (2 (0)
(M ()) using a Gaussian random subspace projec- tion Ω E pxxx where denotes the subspace rank. Then, using the SVD of the approximation Mays (0) = U we can construct a preconditioner (where is a small regularization) as follows:
P= UAAUAA+Dee + UI-
This preconditioner closely approximates the optimal truncated SVD preconditioner when the eigen- spectrum Å resembles that of the original problem. The cost of using this preconditoner are matrix-vector-multinn (MVMs) and Cholesky decomposition of O (?*), where in problems é {100, 200, 300}.
Projection to SGD-like regions Following the ODE on the parameters leads to linear systems whose condition worsens over time seen the middle panel of Figure | In that plot it is also visible how the conditioning of the systems is much lower and does not increase as rapidly when the neural network is trained fit the PDE solution using SGD. In principle would opt for the SGD behavior but we do not have access to fitting the ground truth PDE solution directly. refit against
neural predictions using SGD and thus start from another location in the parameter space. That is, every so often, we solve 8SGD (t) arg min, Sx. (N (0, .) - N (0(t), x²)² du (x) where 0 (t) is our current parameters time t. Performing restarts this way considerably improves the conditioning. As seen in Figure 3 (left) the number of CG iterations increases substantially slower when using the SGD restarts.
Instead, when the condition number grows
large,
4.3 VALIDATING OUR METHOD
We validate our method with three different PDEs: (6+1) and the Pokker-Planck equation (8+1). For more details on these equations, see Appeenpxx B For the wave equation we make comparisons against its analytic solution, while for remaining equations we compare using the PDE residual (evaluated tional details about the experiment setup can be found in appendix B For the wave equation, Neural- IVP achieves the lowest eтror beating EDNN and finite differences evaluated grid as seen Figure 6 For the remaining equations Neural-IVP achives an order magnitude lower PDE residual compared to EDNN
the wave equation (3+1), the Vlasov equation
different samples from training), addi-
100 x 100 100
in the middle and right panels of Figure
While we would like to henefit from the improved sampling distribution proposed by Bruna al. 2002 unfortunately since the solution is prohahility distribution there is no clear measure sample from.
8

---

Published as a conference paper at ICLR 2023
10-1
0.4
0.4
10
0.2
0.2
0.0
N
0.0
0.2
-.2
0.4
0.4
Time
y
Figure 4: (Left): Neural-IVP's PDE residual over time for wave maps compared to EDNN. (Mid- dle): Neural-IVP solution for the wave maps at 0.36 at a slice. (Right): Finite difference solution for the wave maps at the same slice.
SOLVING CHALLENGING HYPERBOLIC PDES
We now turn to challenging PDE: the wave maps equation. This equation is second-order hy- perbolic PDE that often arises general relativity and that describes the evolution of scalar field in a curved (3+1) dimensional spacetime. Following Einstein's 's tensor notation from differential geometry, the wave maps equation can expressed as
g" V V g" 0,0000- 1" 110000
where g is the metric tensor expressing the curvature of the space, 1 the Christoffel symbols, combinations of the derivatives of the metric, and 0, are derivatives with respect to components of both space and time. For the metric, choose the metrie from Schwarzschild black hole located at coordinate value = - 2. in a Cartesian- like coordinate system:
gou dx = -1 - -rr/re)dt² + [8eg + rg(ro/r--1) (I; - c)) - cj)]dr'dr) dx",
where No = || - cll E, C2 and rs 2M is the radius of the event horizon of the black hole, and choose the mass M 1/2 so that r, = 1. We choose wave packet initial condition and evolve the solution for time 7' .5 1M inside the box []1, with the artificial Dirichlet boundary conditions the boundary 1J3 Here the event horizon of the black hole lies just outside the computational domain and boundary conditions meaning that we need not worry about complications on and inside the horizon, and instead the scalar field only feels the effect of the gravity and is integrated for time short enough that it is not yet pulled inside.
While we do not have an analytic solution compare to, we plot the relative error of the PDE residual averaged over the spatial domain in Figure 4 which is consistently small. We also compare the solution the time spatial grid size of 150 150 x 150 which is the largest we were able to run with optimized sparse finite difference solver before rumning out of memory. Despite the challenging nature of the problem, Neural-IVP is able to produce consistent solution for this task. Finally, present an ablation study in Figure 6 showing the gains of using the sinusoidal embedding and of scaling the neural network size and grid for this experiment.
0.36 of our solver solution against the finite difference solution run at
6 DISCUSSION
There are many PDEs of interest that are massively complex to simulate using classical methods due to the scalability limitations of using grids and meshes. At the same time, neural networks have shown promise for solving boundary value problems but the current methods for solving initial value problems can be unstable, deficient scale and in representation power. To ameliorate these defi- ciencies, we presented Neural-IVP, local-in-time method for approximating the solution to initial value PDEs. Neural-IVP a compelling option for problems that are computationally challenging for classical methods like the (3+1) dimensional wave maps equation in section 5) Continuous effort on this front will empower researchers and engineers to simulate physical PDEs which lie at the boundary of what is currently possible solve, allowing prototyping and experimentation without the massive complexity of modern large scale grid-based solvers involving mesh generation, mesh refinement, boundaries, excision, parallelization, and communication.

---

Published as a conference paper at ICLR 2023
REFERENCES
Jens Berg and Kaj Nyström. A Unified Deep Artificial Neural Network Approach To Partial Differ- ential Equations Complex Geometries. Neurocomputing, 1(117)22 41, 2018.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: composable transformations of Python+NumPPy programs. Version 0.3. 17, 2018.
URL httptttttgtttubbbocttto http:///ttthuubbccoooo /aapa//gooole///aaa
Joan Bruna, Benjamin Peherstorfer, and Eric Vaden-Eijnden. Neural Galerkin Scheme with Active Learning for High-Dimensional Evolution Equations. Preprint arXiv 2203.01360vl, 2022.
Laurent Dinh, Razvan Pascanu, Samy Bengio, and Yoshua Bengio. Sharp minima can generalize for deep nets. In International Conference on Machine Learning, pp. 1019-1028. PMLR, 2017.
M. W. M. G. Dissanayake and N. Phan-Thien. Neural-network- based approximations for solving partial differential equations. Communications in Numerical Methods in Engineering, 10(3):195- 201. 1994.
J. Dormand and P. Prince. A Familiy of Embedded Runge-Kutta Formulae. Journal of Compula- tional and Applied Mathematics, 6(1):19-26, 1980.
Yifan Du and Tamer A Zaki. Exolutional Deep Neural Network. Physical Review E, 104(4):045303, 2021.
Weinan E., Jiequn Han, and Amulf Jentzen. Deep Learning Baarni Based Numerical Methods for High- Dimensional Parabolic Partial Differential Equations and Backward Stochastic Differential Equa- tions. Communications in Mathematics and Statistics, 5(4):349-380 2017.
Zachary Frangella, Joel A. Tropp, and Madeleine Udell. Randomized Nyström Preconditioning. Preprint arXiv 2110.02820v2, 2021.
Jiequn Han, Arnulf Jentzen, and Weinan E. Solving High DigeensDmmiii Partial Differential Equations using Deep Learning. Communications in Mathematics and Statistics, Procedigs the National Academy of Sciences, 2018.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Proceedings of the IEEE conference computer vision and pattern recognition, pp. 770-778, 2016.
Sergey loffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In International conference on machine learning, pp. 448 456. PMLR, 2015.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Neural Operator: Learning Maps Between Function Spaces. Preprint arXiv 2108.08481 2021.
I. E. Lagaris, Likas, and D. Fotiadis. Artificial Neural Networks for Solving Ordinary and Partial Differential Equations. IEEE Transactions on Neural Networks, 95):987-1000, 1998.
L. Lu, P. Jin, and G. E. Karniadakis. DeepONet: Learning Nonlinear Operators for Identifying Differential Equations Based on the Universal Approximation Theorem of Operators. Preprint arXiv 1910.03/99 2019.
Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. Nerf: Representing scenes as neural radiance fields for view synthesis. Communications of the ACM, 65(1):99-106, 2021.
Dmytro Mishkin and Jiri Matas. All you need is a good init. arXiv preprint arXiv:1511.066 2015.
10

---

Published as a conference paper at ICLR 2023
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-Informen Neural Networks: A Deep Learn- ing Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differen- tial Equations. Journal of Computational Physics. 2019.
Prajit Ramachandran, Barret Zoph, and Quoc V. Le. Swish: Self-Gated Activation Function. Preprint arXiv 1710.0594lv1, 2017.
Justin Sirignano and Konstantinos Spiliopoulos. DGM: A Deep Learning Algorithm for Solving Partial Differential Equations. Journal of computational physics, 375:1339-1$33 1364, 2018.
Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit Neural Representations with Periodic Activation Functions. Advances in Neural Information Processing Systems, 33:7462-7473, 2020.
Bing Yu et al. The Deep Ritz Method: Deep Learning Based Numerical Algorithm for Solving Variational Problems. Communications in Mathematics and Statistics, 6(1):1-12, 2018.
Learning-Based
A APPROXIMATE SYMMETRIES YIELD SMALL EIGENVALUES
Suppose that the network N has an approximate symmetry in the parameters, meaning that there exists a value for which
Vx [N, T. ((x, 0)1 ca?
holds for a in a neighborhood of 0. If this is the case, we can rearrange the inequality and take the limit as a + 0:
lim [N N(z,))/) a+0
N(x,
(10)
a (2, exists, we can interchange the limit and As the limit lima-00 (1).
norm to get
IVAN u11 -V.VoNN
(11)
since 8a
N(r, (a)))
VoNTo v for v0) = 0, (a) Recalling that (0)
=
I N dp(x), we can take the expectation of both sides of the inequality with respect to μ. producing
T Mu <.
(12)
B EXTENDED EXPERIMENTAL RESULTS
In this section expose additional experimental details that not fully covered in section 4.3
B. WAVE EQUATION
For this experiment use the 3 dimensional wave equation 0. Au, that has few well known analytic solutions. Even this equation computationally taxing for finite difference and finite ele- ment methods. We the radially symmetric outgoing wave solution u(x; () ʃ(L ||=l])/ ||=]||
with f() e 2000.2 2 and integrate the initial condition forward in time by T .5 seconds. In Figure w/ compare to the analytic solution each of the solutions produced by Neural-IVP, EDNN (Du & Zaki 2021). and finite differences evaluated on a 100 x 100 x 100 grid with RK45 set to 10 tolerance. Despite the fact that this initial condition does not contain fine scale detail that Neural-IVP excels at, Neural-IVP performs the best among the three solvers, and faithfully repro-
duces the solution shown in Figure 6 (left).

---

Published as a conference paper at ICLR 2023
B.2 VLASOV EQUATION
The Vlasov equation is PDE that describes the evolution of the density of collisionless but charged particles in electric field, expressed both as function of position and velocity. This equation has spatial dimension and takes the following form
u,u (ar, v, + V (x,0,()+ + E(x, Vuu (V, (I, =
where the vector Rd represents the position of the particles and E IR represents the velocity, and u represents a normalized probability density over I, v. Here is the charge of the particles and is their mass (both of which set 1). In self contained treatment of the Vlasov equation, the
electric field (x, t) is itself induced by the density of charged particles: E(x, t) = V 2 6 t) and the potential the solution to the Poisson equation Aø -p where p(x, L) f qu(a, v, l)..
However, to simplify the setting to a pure IVP, we assume that E(x, 1) is some known and fixed electric field.
For this particular example, we choose E (x) = V exp( |x)12 and the initial condition is product two Gaussians
No (z, = W 3 N) 0, .2).
corresponding to the Maxwell- -ootelll-Blttaannn distribution over velocities and a standard Gaussian dis- tribution over position, and we solve the problem on the cube |-1, 11%
B.3 FOKKER-PLANE
For the Fokker-Planck equation, we choose the harmonic trap for a collection of d = 8 interacting particles from Bruna et al. (2022). giving rise to the equation
0.u (x, t) DAu(z, t) V (hu),
where h2x) = (a x) + a(11 /d T)x. We choose a = (0.2)1 along with constants D = _01 and a 1/4.
We solve this equation dimensions with the initial condition
uo (x) (2) IT - x
which normalized probability distribution on [-1, 1jd. conditions for this problem.
We use the same Dirichlet boundary
Additionally, since we can naturally increase the dimensionality the Bokker-Plandee equation, we explore up to what dimension Neural-IVP can give reasonable PDE residuals. As seen in Figure Neural-IVP can still give solutions for = 20. For dimensions higher that 20, the lincar system solvers cannot converge the desire tolerance needed to evolve the parameters. We warn that these higher dimensional solutions are not guaranteed to be of high quality since the PDE residual estimation also worsens increase the spatial dimension.
C HYPERPARAMETERS
For experiments used the following hyperparameters for Neural-IVP in experiments unless otherwise specified:
RK45 Integrator with rtol: le-4 (for all equations except wave maps, which uses RK23)
2. Number of Monte Carlo samples: 50K for wave maps and 10-20K for all other PDEs
3. Maximum CG iterations: 1,000
4. CG tolerance: 1e-8
Nyström preconditioner rank: 200-350
Linear system regularization: 1e-6
7. Initial fit iterations, optimizer and learning rate: 50K, ADAM, le-3
12

---

Published as a conference paper at ICLR 2023
10-1
10-1
d-5
d-10
d=20
10-2
-0-2
PDE residual
PDE residual
10-3
10-3
EDNN
Neural-IVP
10a5
10-5
0.0
0.2
0.4
0.6 Time
0.8
1.0
0.0
0.1
0.2 Time
0.3
0.4
Figure 5: (Left): Neural-IVP able to provide solutions up to dimension 20 + 1 for Fokker- Planck equation. Higher dimensions break-down the linear systems to evolve the PDE parameters. (Right): Interventions to transform EDNN to Neural-IVP. (II) First, add a sinusoidal embedding before the MLP. (I2) Second, head finetuning (in this case there is notable improvement as the initial condition does not posses finer details). Finally, scale the neural network and the grid size. Note that this is only possible due scalable and efficient construction.
Neural-IVP
EDNN
10-1
0.0
0.2
0.4
Time
10-3
0.0
0.1
0.2
0.3
0.4
0.5
Time
Figure 6: (Left): Neural-IVP fit of the wave equation through time. (Right): Neural-IVP performs slightly better than a finite difference method on the 3D wave equation.
8. Floating point precision: double
9. Number of restarts: 10
The neural network architecture we use is simple MLP with hidden layers, each with 100 hidden units, L for the highest frequency power in the sinusoidal embedding, and the network uses swish nonlinearities. The initial sinusoidal embedding values (p) are scaled by 1.5 before feeding into the network.
13

---

Published as a conference paper at ICLR 2023
Algorithm NEURAL-IVP
1: Input:
1. IVP: Initial condition uo (z), PDE rhs operator C(u) (ar, t), integration time T
2. Design Choices: neural network architecture N (z, 0), sampling distribution μ
3. Hyperparameters: r number of Monte Carlo samples, ODE_TOL, CG_TOL, regularization
J., preconditioner rank r
2: Output: Solution u(x, t) specified times 11,
3: function NEURAL-IVP
4:
4 FitFunction(uy
5:
A/ + 20TODE.TOL
6: while
do
7:
if Sufficient time since last restart then
8:
< FitFunctionnN( 0))
9:
0, At + 1 1 A/ wtt.eeeeeet..ttttptaaaeaee wtt.sssttpDDDaa.css..mmmom A At, ODE.TOL)
RK23
return (N(, O,,) N(:, ONN)
11:
function FITFUNCTION()
12:
arg ming E (x, 0) u(2)ll minimized with Adam
13:
Separate last layer weights W from N(x, 0) = W (x(2) (including bias)
14:
Solve for W from regularized least squares over samples X:
15:
W <
(x) Φ (X) + (7)-11((x)) ((X)
16:
Assemble 0 [-1] W]
return
17:
function DYNAMICS(F)
18:
Let the Jacobian vector product of N (I;, 0) with v taken with respect to 0 be DN(x;, 0)v
19:
Construct efficient MVM M (0) == V 2m Ei= DN 0)vl1
Construct RHS F (8) V 1 [N] B) DN A)e
21:
Construct rank-r Nystrom predonditioner /' using M (0) MVMs
22:
Solve (M (0) + p/)0 f (0) for 0 using conjugate gradients with preconditioner P
23:
return /
D
NEURAL-IVO -VP PSEUDO-CODE
The pseudo-code for Neural-IVP is present in Algorithm 1 For reduced clutter, have omitted the logic for setting step sizes so that will be output by the integrator the specified times 4 22, IN- Sampling RNG state is updated outside the RK23 step but inside the time evolution while loop.
14


Reviews:
[{"review_id": "Nr-MDViuJfZ", "reviewer": "Reviewer_X8sx", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper describes a method which helps improve scalability of the local-in-time ODE solvers.\n", "summary_of_the_review": "The reviewer thinks that while the introduction and the summary of the challenges of scalability of the local-in-time solvers are really clear and well written, the details of the proposed method and  be experimental assessment could be improved (see comments above).\n\n==\nUpdating the scores to recommend acceptance as the authors, in my view, sufficiently addressed the comments", "strength_and_weaknesses": "Strengths:\n- The reviewer really enjoyed the problem statement in sections 1-3; it gives a good grasp on what the problem is and why the current (ML and non-ML based) solvers struggle to achieve good quality solutions for high dimensional problems, as well as complex initial conditions and for higher-order PDEs\n- The problem is really interesting and in many aspects understudied\n\nWeaknesses:\n- Some of the details of the descriptions are failing to give the consistent picture of what the method looks like; see comments 3, 4\n- Experimental assessment: the current experimental assessment lacks detail on why it is better than the existing solvers; see comment 5\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "Nkv-tSWmyT", "reviewer": "Reviewer_qtnF", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a method, called Neural-IVP, for solving partial differential equations (PDEs) using neural networks. They use various helpful techniques to stabilize the PDE solutions, increase the scalability, and improve the representation power of the neural networks. Neural-IVP is a local in time method, meaning that the time dependency is induced by having time-variant parameters in the neural networks. They improve the representation power by using sinusoidal embedding and solving the last layer of the neural networks linearly. In order to increase the scalability and the number of parameters in the neural networks, they propose to use Jacobian-vector-products implemented using automatic differentiation, which makes the complexity of time and memory grow linearly by the number of parameters. ", "summary_of_the_review": "I believe that using neural networks to solve PDEs is an interesting and important problem to pursue. This paper has done a good job in showing the drawbacks of the previous methods for solving PDEs and has tried to address those by using various techniques. However, one main shortcoming of the paper is the lack of high-dimensional problems/examples. One of the selling points of this paper is the improved scalability. So, one expects to see its scalability to high-dimensional settings as well. In the section 5 of the paper, the title say \"high-dimensional hyperbolic PDEs\" but the paper only shows a (3+1) dimensional spacetime example. I am very curious to see if this method can work in hundreds-dimensional problems where neural ODE methods can work without problem? Or, if there is a limit in achievable dimensions, the authors should be clear about it. Also, reporting the training time versus dimension could also be helpful to see. \n\nMinor comments:\n\nThere are some typos throughout the paper. These are some of them:\n\n- be emploted to computer --> be employed to compute\n- due the challenges --> due to the challenges", "strength_and_weaknesses": "Strength:\n\n- Showing the deficiency of the PDE solvers.\n- Suggesting new techniques to overcome the shortcoming of the PDE solvers.\n- Increasing the scalability and conditioning of the PDE solvers.\n- Proving (via theorems) why singular solutions happen (leading to ill conditioning) due to continuous parameter symmetry.\n\nWeakness:\n\n- Limited applications and examples.\n- Lack of evaluation in high-dimensional PDE problems.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "71ew1Vil_0", "reviewer": "Reviewer_ZCbx", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper developed Neural-IVP, an ODE based IVP solver, to solve initial value PDEs. Neural-IVP is based on the local in time methods, but the authors show that local in times methods have issue of representational power and numerical stability. Neural-IVP is developed by addressing these two issues.", "summary_of_the_review": "The paper proposed extensions of local in time method, but the experiments are not convincing.", "strength_and_weaknesses": "Strengths\n- The scalability and numerical stability issues of local in time methods are analyzed.\n- New techniques are proposed to address those issues.\n\nWeaknesses\n- A few techniques are used together to improve the local in time methods. Each of these techniques is not totally new. It is also unclear how important of each technique.\n- There is only one example in this paper. Many more examples are needed.\n- Even for this example, there are not sufficient convincing results. There are no quantitative errors.\n- In the paper, it claims about “high dimensional” PDEs. However, the example in Section 5 is a 3D problem, not high dimensional at all. Examples of at least 10 dim are required.\n- If time dependent 3D problems are “high dimensional”, then many existing papers have solved time dependent 3D problems by deep learning methods, such as physics-informed neural networks (PINNs). A lot of relevant literature review is missing. \n- There are no comparisons with other methods (either deep learning or conventional methods) in terms of accuracy and speed.\n- In Introduction, “Rather than compute costs which scale exponentially in the dimension, they instead typically have error …”. There are no definitions of N and epsilon. This sentence is not a rigorous mathematical statement. Also, references are required.\n- In Introduction, “more that 3 dimensions”. A typo here.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster