Non-Linear Operator Approximations for Initial Value Problems
Authors: Gaurav Gupta
Xiongye Xiao
Radu Balan
Paul Bogdan
Venue: ICLR 2022 Poster (2022)
Keywords: exponential operators
initial value problem
pade approximation
multiwavelets
partial differential equations

Abstract: Time-evolution of partial differential equations is the key to model several dynamical processes, events forecasting but the operators associated with such problems are non-linear. We propose a Padé approximation based exponential neural operator scheme for efficiently learning the map between a given initial condition and activities at a later time. The multiwavelets bases are used for space discretization. By explicitly embedding the exponential operators in the model, we reduce the training parameters and make it more data-efficient which is essential in dealing with scarce real-world datasets. The Padé exponential operator uses a $\textit{recurrent structure with shared parameters}$ to model the non-linearity compared to recent neural operators that rely on using multiple linear operator layers in succession. We show theoretically that the gradients associated with the recurrent Padé network are bounded across the recurrent horizon. We perform experiments on non-linear systems such as Korteweg-de Vries (KdV) and Kuramoto–Sivashinsky (KS) equations to show that the proposed approach achieves the best performance and at the same time is data-efficient. We also show that urgent real-world problems like Epidemic forecasting (for example, COVID-19) can be formulated as a 2D time-varying operator problem. The proposed Padé exponential operators yield better prediction results ($\textbf{53\%} (\textbf{52\%})$ better MAE than best neural operator (non-neural operator deep learning model)) compared to state-of-the-art forecasting models.

Published as a conference paper at ICLR 2022
NON-LINEAR OPERATOR APPROXIMATIONS FOR INI- TIAL VALUE PROBLEMS
Gaurav Gupta Xiongye Xiao Radu Balan 2, Paul Bogdan
1 Ming Hsich Department of Electrical and Computer Engincering University of Southern California
Los Angeles, CA 90089, USA
2 Department of Mathematics and the Norbert Wiener Center for Harmonic Analysis and Applications University of Maryland College Park, MD 20742, USA
{ggaurav, xiongyex, pbogdar)eusc. edu, rvbalan@math. wmw.aanmmaaau.u.mdde..... umd
ABSTRACT
Time-evolution partial differential equations is fundamental for modeling several complex dynamical processes and events forecasting, but the operators associated with such problems are non-linear. We propose a Padé approximation based exponential neural operator scheme for efficiently learning the map between a given initial condition and the activities at a later time. The multiwavelets bases are used for space discretization. By explicitly embedding the exponential operators in the model, we reduce the training parameters and make it more data-efficient which is essential in dealing with scarce and noisy roal-world datasets. The Padé exponential operator uses a recurrent structure with shared parameters to model the non-linearity compared recent neural operators that rely using multiple linear operator layers in succession. We show theoretically that the gradients associated with the recurrent Padé network are bounded across the recurrent horizon. We perform experiments on non-linear systems such as Korteweg-de Vries (KdV) and Kuramoto Jurammot--Svashinskky (KS) equations to show that the proposed approach achieves the best performance and at the same time is data-efficient. We also show that urgent real-world problems like epidemic forecasting (for example, COVID- 19) can be formulated as a 2D time-varying operator problem. The proposed Padé exponential operators yield better prediction results (53 % (52% better MAE than best neural operator (non-neural operator deep learning model)) compared state-of-the- forecasting models.
INTRODUCTION
Predicting the future states using the current conditions is fundamental problem in machine learning. robotics, autonomous aerial ground underwater systems and cyber-physical systems (Xue & Bogdan (2017). Such problems fall under the umbrella of common term, the "Initial Value
Problems (IVPs). The basic structure of IVP involves a first-order time-evolution along with non- linear operators. The class of IVPs spans the domain of physics (modeling gravitational waves Boodgkin-Huxleds model (Zhang et al. (2200), engineering (fluid dynamics Wendt 2008)). water waves (tsunami (Elbanna al. [221]). mean field games (Ruthotto et al. 2020 Bogdan & Marculescu] (01T), to list just a few. Within the current pandemic context, the
(Lovelace 2021) neuroscience
applications areas like epidemiology Beermack-MeKendaee model (Kermack et al. 1991 Diekmann
al. 2021) are of tremendous interest.
Neural Operators The use deep learning to solve the IVP like problems for predictions has been exploiting within the framework of convolutional neural networks (CNNs) (Bhatnagar et al.) 2019 Guo et al. 2016). and time-evolution by employing multiple layers (Khoo et al.) 2020). The multi-layered deep networks with CNNs suitable to solve problems with large number of training samples. Moreover, because the image-regression like structures, such models are restricted to the specifications of the input size. Another research direction aims at solving and modeling the partial

---

Published as a conference paper at ICLR 2022
differential equations (PDEs) versions of the IVPs for given instance. The works of (Kochkov et al. 202I) model the IVP solution as NNs for modeling the turbulent flows. Along the same lines, we have physics-informed neural networks (PINNs) (Raissi et al. 2019. Wang et al.) 2021b) that utilize PDE structure for defining the loss functions. Such models are not applicable within the context of complete data-driven scenario, or for the setups where the exact PDE structure is not known, for example, modeling the climate, epidemic, or unknown physical and chemical phenomena. Finally, we have the works of Neural Operators that are completely data-driven and input-resolution independent
schemes (Li et al. 2020b a Gupta et al. [2021] Bhattacharya et al. 2020 Patel et al. [202T] Most of these approaches tried to efficiently work with the integral kernel operators, for example, Graph
Nyström sampling in (Li et al.) (020b), convolution approximation in (Li et al. (020a), or the multiwavelets compression
(Gupta et al., 202T). Apart from solving non-homogeneous linear differential equations, the PDE operators are mostly non-linear. To tackle the non-linear behavior, these prior works use a multi-cell architecture with non-ling non-linearity (for example, ReLU). To work with the IVP like problems, and be data-efficient, we aim to adopt explicitly the non-linear operator (exponential) that appears in the IVP solutions.
Exponential Operators The exponential of linear transformation has been a subject of research for the last 150 years (Laguerre, 1898). In the simplest form, the exponential operator appears ae at, y(O) = yo as y(t) = yo (for more general examples, see Table With applications ranging from control systems theory (converting continuous to discrete systems) to solving partial differential equations (Cox & Matthews, 2002) Kassam & Trefethen] (2005), the exponential function of operators is a subject of active research. In deep learning, the exponential function to model non-linearity is used in (Andoni et al. 2014). Recently, the exponential operators have also been explored in the field of computer vision for generative flows in (Hoogeboom et al., 2020
as a solution of: dy
Padé Approximation Although one approach to implementing an exponential operator could be attained through a Taylor series representation, this operator function is prone to errors (Abramowitz) & Stegun| 1965). Scale-and aqu-rin (SSQ) methods are commonly suggested approaches to deal with the errors (Lawson 1967). In addition to SSQ, the Padé approximation which represents an analytic function as the ratio of polynomials achieves state-of-the-art accuracy in computing exponential operators (Fasi & Higham/ 2019). Industry standard numerical toolboxes (for example, MATLAB, SciPy) use the Padé approximation based approach to compute the matrix exponential expm (AI-Mohy & Higham| 2009). Matrix exponential via Padé representation requires dense matrix computations (for example, inverse and higher-order polynomials). Such operations are not numerically feasible, in- in-general, for the inputs with large size. However, the commonly used operators like convolution (possibly, multi-layered) have parameters that are fixed beforehand and are much less than the input dimension. A suitable approach, therefore, is a neural architecture based Padé approximation.
Our strategy, in this work, is explicitly embed the exponential operators in the neural operator architecture for dealing with the IVP like datasets. The exponential operators are non-linear, and therefore, this removes the requirement of having multi-cell linear integral operator layers. While with sufficient data in-hand, the proposed approach may work similarly to the existing neural operators with a large number of training parameters. However, this is seldom a feasible scenario for the expensive real-world experiments, or on-going recent issues like COVID19 prediction. Here, the current work is helpful in providing data-efficiency analytics, and is useful in dealing with scarce and noisy datasets (see Section To the advantage of Padé approximation, the exponential of given operator can be computed with the pre-defined coefficients (see Section 2.3) and a recurrent polynomial mechanism.
Our Contributions The main novel contributions of this work are summarized as follows: (i) For the IVPs, we propose to embed the exponential operators in the neural operator learning mechanism. (ii) By using the Padé approximation, we compute the exponential of the operator using novel recurrent neural architecture that also eliminates the need for matrix inversion. (ii) We theoretically demonstrate that the proposed recurrent scheme, using the Padé coefficients, have bounded gradients with respect to (w.r.t)) the model parameters across the recurrent horizon. data-efficiency on the synthetic ID datasets of Korteweg-de Vries (KdV) and Kuramoto Sivassiisii (KS) equations, where with less parameters we achieve sfate-of-the-arr performance. (v) We formulate and investigate the epidemic forecasting as 2D time-varying neural operator problem, and show
(iv) We demonstrate the
2

---

Published as a conference paper at ICLR 2022
Equation
Solution
du Au, u(t 0) Uo (Linear ODE) dt
u(t T) c/Au..
8-u da u(x 0 0) 240 (Heat equation)
u(z, r) c" up(x))
Lu + WI[ uxx 0) xx)
u(x,T) = e0 (x) + fe((-- M((la,, dt
Table Initial Value Problem examples along with their time-evolution solutions. The exponential function of operators appears in the IVP solutions. For Linear ODE, A is the linear transformation.
that for world noisy and scarce data, the proposed model outperforms the best neural operator architectures by 53% and best non-neural operator schemes by 52%.
OPERATORS FOR INITIAL VALUE PROBLEM
We formalize the partial differential equations (PDEs) version of the Initial Value Problem studied in this work in Section 2.1 Section 2.2
summarizes the multi-resolution analysis using multiwavelets deseribes the proposed use of canonical exponential operators and presents a novel architecture using Padé approximation.
for space-discretization. Section 2.3
2.1
INITIAL VALUE PROBLEM
The initial value problem (IVP) for PDEs can be written in its general form as follows.
Ut F(t, E
(I)
u(x, 0) = uo(x), E
where, u, is the first-order time derivative of u. F is time-varying differential operator (non-linear in-general) such that F _- U [0) x B - B with B being a Banach space. Usually, the system in eq. (1) required to satisfy boundary condition such that Bu(x,
0, An Vt in the solution horizon, and 2n is the boundary of the computational region I with B some linear function. Pertaining to our work, the operator map problem for IVP can be formally defined as follows.
Operator Problem Given A and 22 as two Sobolev spaces HSP with s >> 0,p 2, an operator T is and two functions ug and u(x, T). in this work, we take
such that T : AA U. For given the operator map as T (I) u(x, T) with x: Ω.
Ta summarizzs a few examples of the IVP and their solutions. The exponential operators are ubiquitous in the IVP solutions and, therefore, are important to study. One issue, however, is that the exponential operators are non-linear and unlike convolution like operators, there does not exist a general way to diagonalize them (Fourier transform diagonalizes convolution operator) for an efficient
representation. Previous work on neural operators (Li et al., 2020c a. Gupta et al. [222]) modeled the non-linear operators in one way another by using multplee canonical integral operators along with non-linearit (for example, ReLU). In this work, we directly produce an exponential operator approximation. First, discuss efficient basis (multiwavelets) for space discretization of the input/ output functions in Section 2.2
2.2 RULTIIRESOOI ANALYSIS
The multi-resolution analysis (MRA) aims projecting function to a basis over multiple scales. The wavelet basis (e.g., Haar, Daubechies) are some popular examples. Multiwavelets further this operation by using the family of orthogonal polynomials (OPs), for example, Legendre polynomials for an efficient representation over finite interval (Alpert et al. 2002). The multiwavelets are
Time-advection equation with linear operators L.N and non lonelin function f(.). A wide range of problems modeled, for example, Korteweg-de Vries, Srrrmmtt-Sivvssinskky. Burgers' Equation, Navier-Stokes (list not exhaustive).
A mon-linear infegroddiffereniall solution the time-advection equation using semi-group approach Beylkin Keiser 1997 Pazy 1983 Yoshida 1980). A slightly general version is discussed in (Beylkin al. 1998
3

---

Published as a conference paper at ICLR 2022
useful in the sparse representation of the integral operators with smooth kernels. In addition, the multiwavelets also sparsify the exponential functions of the strictly elliptic operators (Beylkin & Keiser (997). However, do not rely this assumption in this work. Here, briefly introduce the MRA and refer the reader to Gupta et al. (202T) for a more detailed discussion.
Notation We begin by defining the space of finite interval polynomials as Va (flf are polynomi- als of degree < k defined over interval (2-1,2201 (1 + 1)) for all / = 0, 1..... 2" 1, and assumes 0 elsewhere). The V are contained in each other for subsequent n or,
V C Vf c... V v c V c....
(2)
The orthogonal component of these polynomial spaces is termed as multiwavelet space w and are defined such that
W
V n++
V L W
The orthonormal basis of V are OPs @u. 01..... Φk-1 and we have used appropriately normalized shifted Legendre Polynomials in this work. The basis for V and W
a are ())()) = 2 6, (2"x -
and (z) 2, (2"x I), respectively. for l 0, 1, ........2. 2 and k
k 1
Finally, an important trick for representing the operator 7 in the multiwavelet basis is called non- standard (NS) form (Beylkin et al.) 1991). The NS form decouples the interactions of the scales and is useful in obtaining efficient numerical procedure. Using NS form, the projection of operator T expanded using a telescopic sum as follows.
Tp 2 (Q:QQ: Q.T.,
P.-) PITPI
where, Pn H.s.2 - V is the projection operator, Tn = P.TPn, On : H=2 W n such that Qn Pr Pa-1. and L is the coarsest scale under consideration (L * 0). Therefore, the NS form of the operator is a collection of the triplets (A,, B₁ C.).-++1 and P,TP, with A; Q,TQ; B, Q,TP:-1 and C; P,-ITQ TQ;- In this work, aim to model A, B,, C, the exponential operators to better learn the IVP by explicitly embedding the non-linear operators into the multiwavelet transformation. This is not straightforward due the non-linearity of exponential functions. We are now in shape to present the main contribution of the current work in the Section 2.3 where we discuss an implementable neural approximation of the exponential operators.
2.3 EXPONENTIAL OPERATOR APPROXIMATIONS
Due to the nature of first-order time-evolution equations, the exponential operators appear in the solution of IVP as discussed Section Being an analytic function, the exponential also assumes Taylor series expansion. However, the approximation error by truncation is (Abramowitz & Stegun 1965) and may require large number of coefficients. We now discuss a better approximation for the non-linear functions.
Padé Approximation Given an analytic function f(2), 2 e C at and let p. q E N, the [p/4] Padé approximation of f at is rational polynomial (m) (2) By (2), where Ap, and (z) are polynomials of degree p and respectively. For our work, f(z) e², and the [p/4] Padé approximation to the exponential function is written as App(2) = Σj az and (2g
N3=0b, where
0j = (p+g-))-p) 0 D<j<p by (p+q j)lا!{ 1)3
(p + g)15 j)! 0 EEję.
(5)
(p (1):(p--))
Note that, Op = by = 1. and we now discuss the exponential of operators. Given an operator L (possibly non-linear), the [p/g] Padé approximation for ex can be written as
of (p) = Σ by LoLo...
Σ ag Co Lo
(6)
j -imess
j times
For L being a linear transformation, the inverse in eq. 6 is computed as a matrix inverse, for example when evaluating matrix exponential. Even convolution operator can be represented the circulant
4

---

Published as a conference paper at ICLR 2022
matrix (shift-invariant kemnel) operation. But storing the big circulant matrix is not numerically efficient solution. Instead, we now discuss a neural architecture that emulates the Padé approximation for operators in eq. (6) while avoiding taking the inverse.
Padé Exponential Model We compute the operator exponential using the [p/q] Padé approximation via recurrent neural architecture in F First, the denominator polynomial Bpg (C) is evaluated using the left recurrent network. Second, the output is passed through a non-linear layer which we implement as aW + b) with o(.) being the ReLU function. Note that the non-linearity layer is applied to the input channels, and therefore, its size is independent of the input spatial dimension. Finally, the output of the non-linewt layer is passed through another recurrent network implementing the numerator polynomial Apg (C). Note that, both polynomial recurrent networks use same operator f. or in other words, the parameters are shared across the network. With aj, bj fxed-beffrehhnd, the total trainable parameters of the Padé exponential model are (0z. W, b)-
One issue with the recurrent architectures is the possibility of gradients explosion when evaluated over a large horizon. In such cases, techniques like gradient clipping is used as workaround. For a given p. q. the proposed network runs the recurrent loops for a total of p q 2. We show that the proposed [p/g] Padé network does not suffer from the issue of gradient explosion, and the boundedness of the gradients is established through the following result:
Theorem 1. Given a linear operator L = [(0c), non-linearity layer o (Wu+ b), and p. q N. al points of differentiabiliyy the gradients of the operation x - F(x; 0z. W; b) = [p/ (p/g]=((x) using the (p/g) Padé network in Figure are bounded in operator norm by
1/2
ne
àc 202
< exp((I<1) (1b1l2 + ||
Σ
z02
θy
< exp(II |=|l]>
(8)
eW
θy
(p 1 11-11
< exp
(9)
2b
The detailed proof is provided in the Supplementary Materials (Appendix E Note that, no 0 is such that ny < M,, where M is the dimension of the input. For example, the convolution operator has a fixed kernel independent of the size of the input, the Fourier-based convolution in (Li et al. 2020a) has km Fourier modes independent of the input dimension. Next, as noted earlier, the non linearity layer a aWWu + b) is applied to the input channels (instead of input spatial dimensions), therefore, the gradient bounds do not scale with high input-resolutions, apart from the dependence from ||>|||2 which could pre-normalized.
Finally, the Padé neural model is integrated with the multiwavelet transform by substituting each of A, B, C, from Section 2.2 with [p/g]ef* [p/ lp/qlefa, [p/g]eC<. The complete flow-diagram by plugging-in the Padé model is shown in Appendix
3 EXPERIMENTS
Here, we empirically evaluate the proposed model and compare against the existing approaches. We consider several synthetic PDE datasets well as a real-world example of pandemic prediction (COVID-19). For the data, a similar input/output structure is used as in the recent works of neural operator architectures. Specifically, the input function u1(x)) and the output function u(r, x) for some pre-specified T > 0, are evaluated at M discretized locations of the domain Ω. This yields a single training sample of (uo(x;), u(T, 2;), 2y E 22, 1 & i & M. In total, we take N training samples, and unless stated otherwise, N = 1000 and we test on 200 samples for the synthetic datasets.
Padé Model Implementation The operator C in Figure I is fixed as a single-layered convolution operator for ID datasets, and 2-layered convolutionn for 2D datasets. For getting the input/output operator mapping, the multiwavelet transform is used only for discretizing the spatial domain. The Padé neural model easily fits into the sockets of the multiwavelet transformation based neural operator as sketched in Figure Tiguree[(Appentt (Appendix The multiwavelet filters are obtained using shifted Legendre OPs with degree k 4. In contrast to the work of Gupta et al. (202T). only single cell of
5

---

Published as a conference paper at ICLR 2022
C
C
11
b2
shared parameters A..
L
C
a1
02
ap
Figure 1: Padé Exponential Model. A recurrent neural architecture for computing the exponential of an operator using [p/g] Padé approximation. The multiplicative scalar coeflicients a, by are beed-beforehand using eq. 5 The non-linear fully-connected layer is used to mimic the inverse polynomial operation.
multiwavelet transform used in the current work because the non-lineanee non-linearity in the operator is explicit
via Padé neural operator. This saves lot of trainable parameters and yields compact model which is suitable when dealing with scarce noisy real-world data as we see in Section 3.3 The numerator / denominator polynomials degrees p/q for the Padé approximation is fixed (p, ② (5, 6) for 1D models and (4, 2) for 2D models. While the authors do not advocate that this is the best possible choice, an ablation study for ID data is performed in Appendix
Benchmark Neural Operators We compare against the recently proposed neural operator works with the proposed Padé exponential model (Padé Exp). The graph neural operator (GNO) was proposed in Li et al. 2020b). A multi-level version of the graph neural network (MGNO) in (Li jet a 2020c LNO low-rank representation of the integral operator kernel and then using multiple layers with non-linearity, which also emulates unstacked DeepONet (Lu et al. 2020). A convolution approximation to the canonical integral kernel and then Fourier transform to diagonalize in (Li et al.) 2020a) as FNO. MWT Leg utilize multiwavelets for spatial projections using Legendre OPs and uses multi-cell structure along with ReLU non-linearity. The Padé Exp model delivers state-of-the-art (Sota) performance parameters, the proposed model shows promise for small datasets as shown Section
range of datasets, both synthetic and real-world. With less number of
Training parameters All neural operator models are trained using Adam optimizer with learning rate of 0.001 and decay of 0.95 after every 100 steps. The loss function is taken as the relative L.2 error. For synthetic datasets train for total of 500 epochs and for real-worke COVID-19 dataset we train for a total of 750 epochs. All experiments are done on an Nvidia A100 40GB GPUs.
3.1 KORTEWEG-DE VRIES EQUATION
The Korteweg-de Vries (KdV) equation one-dimensional non-linear PDE used to model the non-linear shallow water waves. For given field u(x, t), the KdV PDE takes the following form:
Su
tu x (0, 1), t (0,
4 00.5u
δx
(10)
uo(x) u(xx 0).
For Table notations, we have L learns the mapping of the initial condition (x) to the solutions u xz 1]. The initial condition is generated in Gaussian random fields according to un (0, 74 (- 721) -2.5) with periodic boundary conditions. The data set is obtained by solving the equation using the fourth-order stiff time-stepping scheme known as ETDRK4 (Cox & Matthews 2002) with a resolution of 200 and datasets with lower resolutions obtained by sub-sampling the highest resolution dataset.
f(u)
0.2502 and N Z The neural operator
For evaluation, firstly, we vary the total training samples (N) for all the operator models. We sample randomly and uniformly 5 times the training subset from the complete data N 1000, and for each

---

Published as a conference paper at ICLR 2022
10"
Padé Exp MWT Leg FNO
Padk Exp
10
MWT Leg
FNO
Relative aror
Relative
6x
4x 10
4x 10
3
200
400
600
800
1000
200
400
600
800
1000
N
Figure 2: (Left) Number of training samples N performance (relative 1.2 error) for neural operators evaluated on the KdV equation with s=024. For N 1000, each smaller dataset sampled uniformly randomly times from the complete dataset (N = 1000) and mean + std.dev (shaded region) results are shown across the sampling experiments. (Right) Same analysis for KS equation with s=1024.
Networks
64
128
$ 256
512
1024
Padé Exp
0.00301
0.00308
0.00311
0.00298
0.00295
MWT Leg
0.00372
0.00369
0.00391
0.00408
0.00392
FNO
0.00663
0.00676
0.00657
0.00649
0.00672
MGNO
0.12507
0.13610
0.13664
0,15042
0.13666
LNO
0,04234
0.04764
0,04303
0.04465
0.04549
GNO
0.13826
0.12768
0.13570
0.13616
0. 0.12521
Table Korteweg-de Vries (KdV) equation benchmarks different input resolution The relative errors are shown shown for each model.
sampling the models are evaluated. For consistency, each model is fed the same training sub-sample and the results are shown in (Lefft)(Left) We see that Padé Exp model has the sharpest decay compared to the other state-of-the-att neural operators model dltt-efficieet and works well when less data is available. Next, we also evaluate the Padé Exp model on the complete data N 1000 but by varying the input resolution shown in Table The proposed model performs consistently better for the input resolutions.
varying N. This shows that the proposed
3.2 KURAMOTO-SIVASHINGEE (KS) EQUATION
The Kuramoto Surramot-SSvvshinssyy (KS) equation is fourth-order non-linear PDE derived to model the diffusive instabilities in laminar flame front. For given field u(x, (), the KS PDE takes the following form:
2u
2²u
Au dat x (0,1), (0,
We -
dr
dx
(11)
m0(x:)
ul xr,
0).
The KS equation is also time-advection and according to the Table notations, have: L 3
A. T. f(u) -0.5u² and N à The neural operator learns the mapping of the initial condition (1( to the solutions u(z, t 1). Similarly to the KdV equation in Section 3.1 the initial condition is sampled from Gaussian random field ug (O. (4 521)
22
2.5) with periodic boundary conditions. The equation is numerically solved using chebfun package (Driscoll et al. 2014) with a resolution of 210, and datasets with lower resolutions are obtained by sub-sampling the highest resolution data set.
We evaluate the proposed model on the KS equations and compare with the existing works a similar experimental setup as for the KdV equation in Section 3.1 The reduced training experiment results are shown in 2 (Right), where we have uniformly and randomly sub-sampled the training samples from the complete dataset. For consistency, sub-sampled training set. We again observe that the proposed Padé Exp model with its compact structure has the steepest decay of relative 1.2 error compared the recent neural operator works.
have evaluated all models the same
7

---

Published as a conference paper at ICLR 2022
MAE
Networks
Relative 1.2 error
Net. vs FC
R
Padé Exp
1219 1
1752 666
0.0155
0.0034
82.14% (+652K)
MWT Leg
3554
1157
2928
1338
284
209
0.0245
0.0043
62.0% (+18M)
FNO 3D
4213
3391 + 1233
592
157
0.0301
0.0045
54.0% (+1.0MM)
LNO 3D
28502
12698
6586 + 3442
1465
965
0.1056
0.0394
-105.0% (+238K)
Neural ODE
4339
1174
3443 1 1408
443
192
0.0310
0.0069
53.8% (+172K)
Seq2Seq
2798 +: 456
3317 1690
346 83
0.0273
0.0058
63.7%(+188))) 8M)
Transformer
7087
+ 972
6613
+ 2853
1722 + 320
0.0501 + 0.0094
13.4% (+15.2K)
FC
10305 2818
5885 + 1609
1634 + 686
0.0609 0.0111
(37.2K)
Table COVID-19 prediction benchmarks for different networks using 10-fold resampling with mean 1:: std. dev, across folds. The Mean Average Error (MAE) is presented for Confirmed (C), Recovered (R), and Deaths (D) counts averaged across 7 days of prediction for 50 US states. The relative L2 error is the test error for cach model. The last column compares each network vs FC in terms of the total MAE improvement and total model parameters difference.
Confirmod
Recovered
Deaths
350000
actual Pade Exp
gctuol Pade Exp sog22og
actal Pade Exp seg2sog
325000
acq2seg
800
7500
600
7000
Counts
275000
6500
250000
225000
2001
6000
200000
14
Eie (days)
time (days)
time (days)
Confirmed
Recoverod
Deaths
100000
8500
astual Exp
112000
80000
8400
111000-
8500
Couzzs
110000-
8200
40000
8100
109000
20000
8000
105090
7900
14
Eme (days)
time (dayss)
(days)
Figure 3: COVID19 Forecasting. Confirmed, Recovered, and Deaths count forecasting results for the 07/07/20 07/13/20 (chosen arbitrarily) using previous 2 weeks as the input. The Padé Exp prediction and the best non-neural operator scheme from (seqq[[s is shown. Top row: Most populous US state California with population 39.77 M (2018 8 census). Botiom row: Same results for Massachusetts with moderate population of 6.89 M.
By using the exponential operators, the IVPs can be efficiently solved as we witness from Figure Finally, we also evaluate the proposed model for the complete data but varying the input resolution (see results in Table Appendix(B.4] We again conclude that the proposed model attains better performance than state-of-the-at approaches for all the resolutions.
3.3 EPIDEMIC FORECASTING: COVID-19 STUDY
The epidemic forecasting problem refers to the prediction of future counts of infected individuals, recovered individuals and deaths using the current observational data. variety of compartmental models exist model the epidemic spread with certain assumptions (see related work, Appendix AA) The dynamics of the epidemic spread by modeling their time-evolution behavior is not exactly known, and may not apply any problem, for example, the recent COVID-19 19 pandemic. Neural operators, providing a complete data-driven approach, are capable to learn PDE agnostic maps. Consequently, show that the epidemic forecasting can be formulated
operator map learning problem.
8

---

Published as a conference paper at ICLR 2022
Dataset The COVID-19 data set from April 12th 2020 August 28th 2021 is provided by Johns Hopkins University (Dogg et al. 2020). We take the data of 50 US States, and for each state, we have the total counts of daily reported confirmed (C), recovered (R), and deaths (D). We normalize the data of each state by their total population. Therefore, we have a daily collection of 2-dimensional data of size 50 x 3.
Deep learning Benchmarks In addition to the neural operators, we also compare against the state- of-the-art deep learning techniques: an auto-regressive fully connected (FC) network, Sequence to Sequence (Seq2Seq) (Salinas et al.) 2019 Rangapuram et al. 2018) and Transformers (Vaswani al. 2017) utilizing encoder-decoder structure. For Seq2Seq. we have used LSTM architecture for encoder and decoder. Neural ODE (Chen et al. (219) utilizing latent ODE for time-series forecasting.
Operator Map The operator task is to learn the map between the 14 consecutive counts (C, R. D) to next days data for each of the 50 US states. Let de be the x 3 array on day-t, then the operator map can be written as follows.
T(d-14, dd..... ....d....,dd))) (do, d1, .ih....,d!!
(o(
u(T,2)
Forecasting The COVID19 forecasting benchmarks are presented Table Due to data scarcity (484 samples in total), we do 10-fold resampling of the dataset to obtain train/test samples and the averaged results are presented for all models. We see that the proposed Padé exponential model achieves better performance than existing approaches especially in the presence of scarce and noisy data setup. The Seq2Seq performs best among all non-neural operator models. In terms of the total mean averaged error (MAE) (for C, R, and D counts), the Padé exponential achieves a 53% improvement over the best neural operator (MWT Leg), and 52% over the best non-neural operator model (Seq2Seq). The percentage improvement over the FC model and the difference between the total model parameters w.r.t. FC are shown in the last column. A sample forecasting is shown for 2 US states in Figure We also show the corner cases (best and worst test sample) predictions in the Figure (Appendix B.5 along with the best/worst prediction states.
4 ABLATION EXPERIMENTS
The ablation experiments are presented in the Supplementary Materials ( (Appendix B.1) The follow- ing experiments are performed: (i) Comparison between the proposed Padé approximation vs Taylor Series neural operator. (ii) the variations of p, q for the [p/g] Padé approximation model in Figure]] and (iii) variation of the non-linearity module in the Padé neural model.
FUTURE DIRECTIONS AND CONCLUSION
Time-evolution analysis of PDEs and forecasting the future states from a set of observations is fundamental for studying wide range of complex systems in physics, chemistry, geoscience, neuroscience, system biology, social, political, and climate sciences. In many such contexts, we need to solve an initial value problem consisting of first-order time-evolution of several nonlinear spatial operators. To efficiently solve these initial value problems while also overcoming data science challenges (e.g., scarcity in the number of samples, samples corrupted by unknown noise types and sources), we proposed a combined multiwavelet and Padé approximation based exponential neural operator architecture. The proposed model has magnitue fewer parameters (see Table 3) and attains data-efficiency.
The IVPs deal with first-order time evolution, and thus, the exponential operator appears. In more general setup, the higher-order time derivatives, or even fractional time derivatives should be considered to model the Man-Manoovian dynamics. In such cases, Reisz transform like approach helps. This is an interesting future direction for neural operator design, where data helps identifying the time-evolution order.
https://githhii

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGEMENT
We are thankful to the anonymous reviewers for providing their valuable feedback which improved the manuscript. We gratefully acknowledge the support by the National Science Foundation Career award under Grant No. CPS/CNS-145386) the NSF award under Grant CCF- 183-18 MCB-1936775, CNS-1932620, the U.S. Army Research Office (ARO) under Grant No. W91 INF- 17-1-0076, the Okawa Foundation award, the Defense Advanced Research Projects Agency (DARPA) Young Faculty Award and DARPA Director Award under Grant No. N66001-17-1-4044, an Intel faculty award, Northrop Grumman grant, and Google cloud. A part of this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548561 The third author has been supported in part by NSF award under grant DMS- 100900 and by the Simons Foundation. The views, opinions, and/or findings contained in this article those of the authors and should not be interpreted representing the official views policies, either expressed or implied by the Defense Advanced Research Projects Agency, the Army Research Office, the Department of Defense the National Science Foundation.
10

---

Published as a conference paper at ICLR 2022
REFERENCES
M. Abramowitz and LA. Stegun. Handbook of Mathematical Fumctions: With Formulas, Graphs, and Mathematical Tables. Applied mathematics series. Dover Publications, 1965. ISBN 9780486612720.
Awad Al-Mohy and Nicholas Higham. new scaling and squaring algorithm for the matrix exponen- tial. SIAM Journal on Matrix Analysis and Applications, 31, 01 2009, doi: 10. 177009999077
B. Alpert, G. Beylkin, D. Gines, and L. Vozovoi. Adaptive solution of partial differential equations in 2002. ISSN 0021-9991.
multiwavelet bases. Journal of Computational Physics,
Bradley K. Alpert and Vladimir Rokhlin. A fast algorithm for the evaluation of legendre expansions. SIAM Journal on Scientific and Statistical Computing, 12(1):158-179, 1991. doi: 10.1137/ 191200
Kevin Amaratunga and John Williams. Wavelet based green's function approach 2d pdes, Engi- neering Computations, 10, 07 2001. doi: 10.1 108/eb0233113.
Alexandr Andoni, Rina Panigrahy, Gregory Valiant, and Li Zhang. Learning polynomials with neural networks. In Eric P. Xing and Tony Jebara (eds.), Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pp. 1908-1916, Bejing, China, 22-24 Jun 2014. PMLR.
M. Arioli, B. Codenotti, and C. Fassino. The pade method for computing the matrix exponential. Linear Algebra and its Applications, 1996.
G. Beylkin. On the representation of operators in bases of compactly supported wavelets. SIAM Journal on Numerical Analysis, 17661117, 1992. doi: 10. 10.1137/7729097.
G. Beylkin, R. Coifman, and V. Rokhlin. Fast wavelet transforms and numerical algorithms i. Communications on Pure and Applied Mathematics, 44(2):141-18,, 1991. doi: https://doi.org/10. (002/cpa.1100440000
Gregory Beylkin and James M. Keiser. On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases. Journal of Computational Physics, 132(2):233-255 1997. ISSN 0021-999) doi: https://doi.org// wttpg//.000////p///.oooo
1006/jeph. 1996.556..
Gregory Beylkin, James M. Keiser, and Lev Vozovoi. A new class of time discretization schemes for the solution of nonlinear pdes. Journal of Computational Physies, 147(2):362-3877 1998. ISSN 0021-9991. doi: https://doi.org/ wtw.r///...000660 1006/jchh. 19w.ph...99...0. 1998.6093.
Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik. Predic- tion of aerodynamic flow fields using convolutional neural networks. Computational Mechanics, 64(2):525-544 Jun 2019. ISSN 1432-0924. doi: 10. (007///046-019-0072700
Kaushik Bhattacharya, Bamdad Hosseini, Nikola B. Kovachki, and Andrew M. Stuart. Model reduction and neural networks for parametric pdes, 2020.
Paul Bogdan and Radu Marculescu. fractional calculus approach to modeling fractal dynamic games. In 20/ 50th IEEE Conference Decision and Control and European Control Conference, pp. 255-260. IEEE, 2011.
Nicolas Boulle, Yuji Nakatsukasa, and Alex Townsend, Rational neural networks, In H. Larochelle, M. Ranzato, Hadsell, M. Balcan, and H. Lin (eds.), Advances in Neural Information Process- ing Systems, volume 33, pp. 14243-14253. Curran Associates, Inc., 2020.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differentiall equations, 2019.
Yi-Cheng Chen, Ping-Ei En Lu, Cheng-Shang Chang, and Tzu-Hsuan Liu. A time- sir model for covid-19 19 with undetectable infected persons. IEEE Transactions on Network Science and Fnggneering, 7(4):3279-327 3294, Oct 2020. ISSN 2334-329X.

---

Published as a conference paper at ICLR 2022
Steven M Cox and Paul C Matthews. Exponential time differencing for stiff systems. Journal of Computational Physics, 176(2):430-333 2002.
Jessica Davis, Matteo Chinazzi, Nicola Perra, Kunpeng Mu, Ana Pastore Piontti, Marco Ajelli, Natalie E. Dean, Corrado Gioannini, Maria Litvinova, Stefano Merler, Luca Rossi, Kaiyuan Sun, Xinyue Xiong, M. Elizabeth Halloran, Ira M. Longini, Cécile Viboud, and Alessandro Vespignani. Estimating the establishment local transmission and the eryptic phase of the covid- 19 pandemic in the usa. medRxiv, 2020.
Odo Diekmann, Hans G. Othmer, Robert Planqué, and Martin C. J. Bootsma. The discrete-tine kermack-mckendriik model: versatile and computationally attractive framework for modeling epidemics. Proceedings of the National Academy of Sciences, 118(39), 2021. ISSN 0027-8424. doi: 10.10733pnas:210633211
Ensheng Dong, Hongru Du, and Lauren Gardner. An interactive web-based dashboard to track covid-19 in real time. The Lancet infectious diseases, 20(5):533-534, 2020.
A Driscoll, N. Hale, and L. N. Trefethen. Chebfin Guide. Pafnuty Publications, 2014.
Ahmed Elbanna, Mohamed Abdelmeguid, Xiao Ma, Faisal Amlani, Harsha S. Bhat, Costas Synolakis. and Ares J. Rosakis. Anatomy of strike-slip fault tsunami genesis. Proceedings of the National Academy of Sciences, 118(19), 2021. ISSN 0027-8424. doi: 10.1073//naa..2255531111 2025632118
Wyman Fair and Yudell L. Luke. Padé approximations to the operator exponential. Numer: Math., 14 (4):379-382, March 1970. ISSN 0029-599X. doi: 10.1007/BF02116552.
Yuwei Fan, Jordi Feliu-Faba, Lin Lin, Lexing Ying, and Leonardo Zepeda-Nunez. A multiscale neural network based on hierarchical nested bases, 2019.
Massimiliano Fasi and Nicholas Higham. An arbitrary precision scaling and squaring algorithm for the matrix exponential. SIAM Journal on Matrix Analysis and Applications, 40: 1256- 01 2019. doi: 10.1137/18M122887
Jordi Feliu-Fabà, Yuwei Fan, and Lexing Ying. Meta-learning pseudo seeudo-differential operators with deep neural networks. Journal of Computational Physics, 108:109309, May 2020. ISSN 0021-9991. doi:
Silviu-loan Filip. Aurya Javeed, and Lloyd Trefethen. Smooth random functions, random odes, and gaussian processes. SIAM Review, 31:185-205, 01 2019. doi: 10. 0.1133//7M11 161853.
Xiaoxiao Guo, Wei Li, and Francesco lorio. Convolutional neural networks for steady flow ap- proximation. Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD 16. pp. 481-496 Association for Computing Machinery, 2016.
Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavelet- operator learning for differential equations, 2021.
Emiel Hoogeboom, Victor Garcia Satorras, Jakub Tomczak, and Max Welling. The convolution exponential and generalized sylvester flows. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 18249-18260 Curran Associates, Inc., 2020.
Aly-Khan Kassam and Lloyd Trefethen. Fourth-order time-stepping for stiff pdes. SIAM J. Scientific Computing, 26: 1233, 01 2005. doi: 10. 10.1137751064875524116333
W. O. Kermack, A. G. McKendrick, W. O. Kermack, and A. G. McKendrick. Contributions to the mathematical theory of epidemics-I. 1927. Bull Math Biol, 53(1-2):33-55 1991.
Yuehaw Khoo, Jianfeng Lu, and Lexing Ying. Solving parametric pde problems with artificial neural networks. European Journal of Applied Mathematics, 32(3):421-435, Jul 2020. ISSN 1469-4425.
Dmitrii Kochkov, Jamie A. Smith, Ayya Alieva, Qing Wang. Michael P. Brenner, and Stephan Hoyer. Machine learning-acceleratee accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences, 118(21), 2021. ISSN 0027-8424, doi: 10.1073/pnas.210 21017841
12

---

Published as a conference paper at ICLR 2022
S.G. Krantz. A Handbook of Real Variables: With Applications Differential Equations and Fourier Analysis. A Handbook of Real Variables: With Applications to Differential Equations and Fourier Analysis. Birkhääsee,, 2004.
E. Laguerre. Le calcul des systémes linéaires, extrait d' d'une lettre adressée m. hermite. Extrait du Journal de I' l'École Polytechnique, LXII. Cahier, pp. 221-267, 1898. URL http:://gglll// bnf. Er Er/ark ark / 6k000//0000000000000/ Ek90210p/E242
Douglas Lawson Generalized runge-kutta processes for stable systems with large lipschitz constants. SIAM Journal on Numerical Analysis, 4(3):372-38)) 1967. ISSN 00361429.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametriic partial differential equations, 2020a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations, 2020b.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. Advances in Neural Information Processing Systems, volume 33, pp. 6755-6766, 2020c.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Markov neural operators for learning chaotic systems, 2021.
Geoffrey Lovelace. Computational challenges in numerical relativity in the gravitational wave era.
Nature Computational Science, 1(7):450-4522 Jul 2021. ISSN 2662-8457.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem operators, 2020.
Alejandro Molina, Patrick Schramowski, and Kristian Kersting. Padé activation units: End-to-end learning of flexible activation functions in deep networks. In International Conference Learning Representations, 2020. URL /ttpt::t/////tto /ttttp::// net wttp://wwwwee/oeuuneueme wtw.orummmm....1........
Ravi G. Patel, Nathaniel A. Trask, Mitchell A. Wood, and Eric C. Cyr. A physics-informeee operator regression framework for extracting data-driven continuum models. Computer Methods in Applied Mechanics and Engineering, 373:1 13500, 2021. ISSN 0045-7825. doi: http:://do..og/1011100 cma.2020.11 1.300..13500..... 13500.
A. Pazy. Semigroups of Linear Operators and Applications to Partial Differential Equations. Applied mathematical sciences. Springer, 1983. ISBN 9783540908456.
Sen Pei and Jeffrey Shaman. Initial simulation of sars-cov2 spread and intervention effects in the continental us. medRxiv, 2020. doi: 10. 00.1101/222000322222400000000000
M. Raissi, Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019. ISSN 0021-999)
Syama Sundar Rangapuram, Matthias W Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, and Tim Januschowski. Deep state space models for time series forecasting. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
Lars Ruthotto, Stanley J. Osher, Wuchen Li, Levon Nurbekyan, and Samy Wu Fung. A machine learning framework for solving high-dimensional mean field game and mean field control problems. Proceedings of the National Academy of Sciences, (1117)))918399999 2020. ISSN 0027-8424. doi: 10.1073/pnas. 1922204117.
David Salinas, Valentin Flunkert, and Jan Gasthaus. Deepar: Probabilistic forecasting with autore- gressive recurrent networks, 2019.
13

---

Published as a conference paper at ICLR 2022
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention all you need, 2017.
Qinxia Wang. Shanghong Xie, Yuanjia Wang, and Donglin Zeng. Survival-convolution models for predicting covid-19 cases and assessing effects of mitigation strategies. medRxiv, 2020. doi: 10.1 1e.11000222004.1..222
Rui Wang, Danielle Maddix, Christos Faloutsos, Yuyang Wang, and Rose Yu. Bridging physics-based and data-driven modeling for learning dynamical systems, 2021a.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with physics-informet deeponets, 2021b.
J. Wendt. Computational Fluid Dynamics: An Introduction. von Karman Institute book. Springer Berlin Heidelberg, 2008. ISBN 9783540850557.
Yuankun Xue and Paul Bogdan. Constructing compact causal mathematical models for complex dynamies. In 2017 ACM/IEEE 8th International Conference on Cyber-Physical Systems (ICCPS), pp. 97-108, April 2017.
K. Yoshida. Functional Analysis. Classics in mathematics Springer. World Publishing Company, 1980. ISBN 7780387102106
Mingliang Zhang, Menghua Man, Guilei Ma, Meiyu Ye, and Shanghe Liu. Research on action behavior of neuron system in case of single pulse stimulus. Scientific Reports, 10(1):1240, Jan 2020. ISSN 2045-2322.
D. Zou, R. Balan, and M. Singh. On Lipschitz Bounds of General Convolutional Neural Networks. IEEE Trans.on on Info.Theor Theory, 6638))1738-1755) 2020a. doi: 10.11099/1..2211222611111
Difan Zou, Lingxiao Wang, Pan Xu, Jinghui Chen, Weitong Zhang, and Quanquan Gu. Epidemic model guided machine learning for covid-19 forecasts in the united states. medRxivv 2020b. doi: 10.1 101/2020.05 05. 2012201199
14

---

Published as a conference paper at ICLR 2022
Networks
X 0.05
X 0.04
x= 0.03
0.02
Padé Exp (p=5, q=
0.00199
0.00237
0.00219
0,00215
Padé Exp (p=4, q=2)
0.00222
0.00231
0.00291
0.0159
Taylor (L,=6)
0.00730
0.00690
0.00846
0.00704
Taylor (L,=4)
0.00938
0.00842
0.00866
0.0193
L'Hospital (Л==)
0.00775
0.00699
0.00752
0.0194
L'Hospital (n=4)
0.00731
0.00798
0.00844
0.0199
Table 4: Taylor Padé Approximation for exponential operator evaluated the KdV equation for s=1024 with varying input (uo (x)) fluctuation strength (low À implies high fluctuation). The highest polynomials degree are similar for each model by setting Lt max(p.g) q) for Taylor and similarly n max q) for L'Hospital.
0.0036
0.0034
0.0032
0,0030
7
0.0028
0.0026
7
Figure 4: Polynomials degree p; grid variation for [p/9] Padé neural model of Figure evaluated KdV equation with s=1024.
The performance
A RELATED WORK
Epidemic Forecasting In Epidemiology, the compartmental models very common, where the population is assigned into discrete categories (for example, Susceptible, Infected, Recovered), also termed SIR models which date back 1927 (Kermack al., 1991 Diekmann et al. 202T). Along the lines, several works exist that model the time-evolution of epidemic using various assumptions. The work of Chen al. (020)) uses ridge regression based prediction. An extension of SEIR (E for Exposed) compartments to sub-populations moving across different places in (Pei & Shaman 2020). convolution filter based approach in (Wang al. 2020), where the total infected cases are modeled as convolution of novel cases with proportion of infected counts over time. A susceptible EIR model in (Zou al. 2020b) incorporating the unreported cases. The work of Davis al. (2020) which simultaneously models the transmission rate among US states. A physics-based model in (Wang et al. 2021a) that uses dito-diifernntiiaioon using Runga-Kun Kutta methods to estimate the model parameters.
B ADDITIONAL RESULTS
B. ABLATION STUDY
Taylor vs Padé Exponential Approximation We compare the Padé approximation based neural model in Figure vs Taylor series approximation for evaluating the exponential of the operator The Taylor series version is written as
L
e (CoLo...of)
j-mmms
where, Lg is the truncated length of the series. We evaluate the truncated Taylor series approximation vs [p/q] Padé approximation in Figure for same number of maximum polynomial degree, i.e., L, = max(p. q). For both models, we give input with varied degree of fluctuations (to make the data more challenging) and the output is governed by the KdV equation (see Section The fluctuated
15

---

Published as a conference paper at ICLR 2022
input is sampled from Random fields Filip al. (2019) with controllable parameter À with lower values attaining higher fluctuations and the input un ()) reaches the Brownian motion limit with A - O. We see in Table that that the Padé approximation model works better compared to Taylor series based expansion for incorporating the exponential operators for each X.
L'Hospital Approximation The exponential exp(x:) also be approximated as limm +a/n)" using L' L'Hospital rule (Krantz) 2004). For a given operator L, the L' L'Hospital based exponential operator approximation can be written as
ef [ C/n) (I C/n) C/n),
n-times
where, I is the identity operator and n is the length of truncation of operator compositions. We compare L'Hospital approximation with the Padé model Table/ the same degree operator polynomials, i.e., max(p, q). For the same experimental settings taken in the Taylor approxi- mation comparison, we see that the proposed Padé model performs better than the L. L'Hospital-basee exponential approximation for all fluctuations strengths.
Varying Padé polynomials degree The numeratorrdeenoiinator polynomials degree in eq. (5) for [p/g] Padé approximation model in Figure are hyper-parammters.. We vary the p, q over 2D
grid of 1, 10] x 10] and evaluate for the KdV equation settings (s=1024) as mentioned in the Section 3.1 and the results are shown Figure We make the following observations, namely; (i)
For p = q= 1. i.c., when there is no operator L and only the non-linear layer rWWu + b), the model has the worst performance of 0.1063. This shows the importance of operator By increasing the 9 in either direction, the performance improves as we attain better approximation with higher degree polynomials. The best performance is achieved as (p, q) = (10, 1). (iii) The performance is roughly similar in the mid-range of and g's, as long as p >> 1 and experiment on KdV. Note that, although larger values of p and are preferable but they also incur the additional cost of run-time which scales linearly with p g as evident from the recurrent structure in F
the Padé model. (ii)
for the ID
Structure of the non-linear MLP We have used the non-linear layer ( + b) in the Padé neural model for all the experiments with or. (.)==eLU.. Now, vary the depth of this layer from 1 to 3 (each layer is followed by ReLU non- non-linearity) to see the individual contribution of this layer towards the final test performance. The experiments settings are same as for Table and we fix 1024. Upon increasing the depth we do not see any improvements, in fact, after having three layers, the performance starts to degrade.
Depth 1.2 error
0.00295
0.00339
0.00628
Table 5: Non-linear layer depth test error for the Padé neural model.
possible reason could be Denominator polynomial getting affected by the deeper non- non- con-linearity layer.
B.2 PREDICTION AT HIGHER RESOLUTIOON
Test
2048
S 4096
8192
Train
s= 128
0.0423 (0.0473)
0.0440 (0.0511)
0.0450 (0.0544)
s=256
0.0229 (0.0315)
0.0250 (0.0374)
0.0263 (0.0427)
s=512
0.0124 (0.0230)
0.0148 (0.0305)
0.0162 (0.0372)
Table 6: Padé exponential (MWT Leg) model trained lower resolutions can predict the output at higher resolutions.
B.3 TRAINING/EVALUATION WITH DIFFERENT (p. q)
Due to the recurrent structure of the Padé approximation model in Figure with only trainable parameters {0z, W. b). the model can be trained and tested on different values of (p. q) by just varying the recurrence loop lengths of the numerator and denominator polynomial. We setup this experiment on the KdV equation (same setting as Section The Padé model trained for (p,4) (5, 6) and (4, 2) and then tested on [1, 10] [1, 10] grid of p, q tuples shown in Figure
16

---

Published as a conference paper at ICLR 2022
0.8
0.8
3
3
-0.6
0.6
5
5
0..
0.4
0.2
0.2
9
9
3
7 9
1
5
7
9
Figure 5: Relative L. error when Padé model trained and tested on different (p. g). (Left) The KdV data using (p,g) (5, 6) for training and other (p.g) for testing. (Right) The KdV data using (p, g) (1, 2) for training and other (p. 7) for testing. For each plot, the Orange/Green color denotes best performance (lowest value) in the corresponding row/column, respectively.
We observe that the best performance (lowest relative L2 error) results when the train and test (p, g) matches. In addition, we observe that for a fixed g. the best performance results when is such that p/q ratio is similar the training ratio. For example, in Figure (Left), value for each column is obtained when p/q is closest to 5/6 while (Right) we see lowest value of each column when p/q 2. The boundary cases diverge because either p. does not have any operator C in the corresponding polynomial as evident from eq. (6). When trained for (p, q) (4, 2), the best performance for
that lowest
saturates to p 10 because of the maximum limit on the experimental grid (p 10). Overall, this experiment suggests that the relative length of the polynomials during evaluation should be similar to training setting for best performance.
B.4 KURAMOTO-SIVASIINEE
EQUATION
Comparison of relative error for different neural operators evaluated over KS equation is presented in To
B.5
COVID-19 FORECASTING
Additional forecasting results for COVID-19 by considering the corner cases are presented in Figure
B.6 2D NAVIER-STOKEE EQUATIONS
The Ravier-Stokes (NS) equations are 2D time-varying PDEs describing the motion of viscous fluid substances. The NS equations can describe many physical processes and have wide range of practical uses. In this paper, to compare with the state-of-the-an models (Li et al. 2020a Gupta et al. 2021 under the same conditions, we the same data sets that have been published in (Li et al. 2020a). where the NS equations take the following form:
wr + u(x, t) Vw(a, t) - vAw(a, f) = f(x),
x (0,1) (0,7)
V V.ul u(x, f) = 0, wo(x) w(I, 1 0),
x (0, 1) [0,
(12)
JE (0, 1)²
where u is the velocity, w is the vorticity such that V u. The incompressible flow is modeled We the experiments to the neural operator map the first 10 time units to last T 10 time units of vorticity The initial condition is generated in Gaussian random fields according wy (0, 7 ( + 721)- 227-2..5 with periodic boundary conditions and + x2)) + cos(2 + Ø2))). We experiment with
via divergence condition as V u(x,
the forcing function is f(z) = 0. (sin
different viscosities v, final time T, and the number of training pairs N: (i) le 1e-3. 3, T 50,
17

---

Published as a conference paper at ICLR 2022
V 1e 4
1e-4
Networks
T = 50
T 30
T = 30
T = 20
N 1000
1000
N 10000
N 1000
Padé Exp
0.00621
0.1427
0.0619
0.1533
MWT Leg
0.00625
0.1518
0.0667
0,1541
FNO-3D
0.0086
0.1918
0.0820
0.1893
FNO-2D
0.0128
0.1559
0.0973
0.1556
U-Net
0.0245
0,2051
0.1190
0.1982
TF-Net
0.0225
0,2253
0.1168
0.2268
Res-Net
0.0701
0,2871
0.2311
0.2753
Table 7: Navier-Stokes Equation validation various viscosities and prediction horizon T.
Networks
N 200
N 400
N 600
N =800
N 1000
Padé Exp
0.0084415.1e-4
0.0043911 2.8e-4
0.00365 2.2e-4
0.003221 1.7e-4
0.00295
MWT Leg
16.19811611e-4
0.00463 0.77-4
0.00420+2.7e-4
0.00392
FNO
0.00
0..37
0.007
0.0067911.2e-
0,00672
Table 8: Korteweg-de Vries (KdV) equation benchmarks for different numhers of training samples N. Top: Our method. Bottom: the state-of-the-art methods.
Networks
N 200
N 400
N 600
N =800
1000
Padé Exp
0.007
0.00489 2.2e-4
0.004161 1.8e-4
0.00376 1.0e-4
0.00338
MWT Leg FNO
0.00849::5.7-44
0.00612155.e-4
0,00496:4.22-4
0,00478::33e-4
0.00445
1.11024:1.1e-3
0.00771
0.00625+1.7e-4
0.00508:1144e
0.00457
Table 9: Kuramoto-Sivashinsky (KS) equation benchmarks for different numbers of training samples N. Top:
Our method. Bottom: the state-of-the-art methods.
N = 1000; (ii) // le 4,7 = 30, 1000; (tii) v = le - 4,T 30, N = 10000; (iv) 1e - 5, 20,N 1000 on 64 64 grid.
The time-varying 2D data is modeled as 3D operator learning problem. We implemented the Padé exponential model by taking C as 2-layered 3D CNNs with ReLU non-linearity and p q 4 in
Figure A total of 4 layers of multi- wavelet skeleton (Figure with & 3 are concatenated using ReLU non- linearity. The results are reported in Table and we observe that Padé model performs well compared the recent -fate-oflie sattee-of-hh-aate art neural operator approaches. For the less data setup of N 1000 with in F that the Padé Exp model quickly converges to the lowest value compared to the other neural operator
10-1
Relative error
3 10
Padé Exp
MWT Leg
FNO
x
1e 4, we see
100
200
300
400
500
approaches.
epochs
B.7 NUMERICAL VALUES
Figure 6: Relative 1.2 error vs epochs for MWT Leg with different number OP basis
The numerical values for the training samples variation experiment for the KdV and KS equation in the Figure are shown Tnd Table respectively.
18

---

Published as a conference paper at ICLR 2022
Networks
64
s= 128
$ 256
s=
$ 1024
Padé Exp
0.00359
0.00326
0.00347
0.00348
0.00338
MWT Leg
0.00445
0.00414
0.00436
0.00485
0.00445
FNO
0.00461
0.00451
0.00469
0.00491
0.00457
MGNO
0.10362
0.12038
0.13361
0.13343
0.13799
LNO
0.04133
0.04020
0.04498
0.04341
0.04360
GNO
0.14037
0.14277
0.13862
0.14525
0.14363
Table Kuramoto-Siin Kuramoto-Sivashinsky (KS) equation benchmarks for different input resolution emros are shown shown for each model.
relative 1.2
C NOTATIONS
Operator Learning
T,C C
Operators between function spaces
H
Sobolev spaces such that constituent functions and their weak derivatives upto order have finite L norms
Multiwavelets
+
Subspace addition
V
{fli are polynomials of degree k defined over interval (2-1,22 "(I + 1)) for all I 0.1...., and assumes 0 elsewhere)
w
Orthogonal space V such that w + Vh V V
Pp
Projection operator such that Pn H V
Qn
Projection operator such that Qn 74.2 W
L
Coarsest scale of the multiwavelet transform Padé Approximation
1.11
Operator norm
1.12
Euclidean L-2 norm
ip/q1f
Padé approximation of function f with numeratoor/denominator polynomial degree p/q, respectively
Oc
Parameters used to represent the operator
On
Number of parameters used to represent the operator C
y F(:
Mapping from using function F with parameters set θ
D MULTIWAVELET EXPONENTIAL OPERATOR ARCHITECTURE
The Padé neural operator based multiwavelet transform model is shown in Figure The input the model is (n++)) where n = (og(M),, and the output (!!)) where n finest scale (or the log of input resolution). For detailed description of multiwavelet transform, we refer the reader to (Gupta et al. 202I)
E PROOF OF THEOREM 1
First note that, for Padé approximation [p/q]e", the polynomial coefficients in eq.
5
we have
1 p(p (p
(
(13)
pp
19

---

Published as a conference paper at ICLR 2022
U
U(n)
8 ((++1)
s5(2)
H,
U U.
H G
p
L))
T
p//gle_
U Uln-1) 1) U(n-1)
1) U(m
(()
(p/gleCa
Figure 7: Multiwavelet Padé Exponential Model. The Padé exponential neural operator from Figure is used in the skeleton of the multiwavelet transform based neural operator model. The inputs and outputs recursively updated using the decomposition cell (Left) and reconstruction cell (Right).
I
and similarly,
1 q
1by1 & j! + q
For the simplicity of notation: let A(2) = App(2) = E)= az2d, A(22) = Ej= C(2)
1by|=,, (²)) و Note |By(2)] C([=l) and [B2g(z)] & C'([=l). Using eq.
can further write that
|A(2) & exp
p [2l
p P_|z|
[A'(2) &
exp
P
|2|
(14)
p p+q
p
[C(=] <
exp
4+ |2l)
|C'()) < p+
exp
(15)
Further let d and na denote the dimensions of and @ respectively: x, y, VE Rd and 0 = 0c E RR
For the operation x + y, we denote y Fc b,b) = (//]e=(x). Then y AL)),,
oW u + 0), u- (()),, and
Fl lim [F(x; 0 + tg. W,b) F(z: 0, 0, W,
max 9ER .1oll
20c
maxx
KVoK,, 1)
(16)
9 E I 1g
h Rd. |
All matrix norms are spectral (i.e., the largest singular value), unless otherwise indicated as in equation
25 whereas all vectors norms are Euclidean (i.e., the L2-norm). Denote by ēg Ei= 9j an, the
differential operator induced by g. The quantity of interest is (0 hi Next, by chain rule yy h) ((,g, h) (A(C), h) = ((0), h) + O(CCD,
(17)
where Dg is diagonal matrix of I's and 0's depending upon the signatures of entries in W b (see
for instance the computation of Lipschitz constants in Zou et al. (2020a)). Hence ||D || = 1. Since o
(ReLU) is contractive with Lipschitz constant 1, 1|0//2
+16b2 & |[W]] B(2) + 1b/2-
Substituting back in eq. (I7) we obtain
10 h>l & 16 (C))A(((]]]I +
+ ||A(CC)]]
(18)
Using eq.
15 the spectral norms are bounded further by
P [|]
|A(C]] & A(ICI) & exp
(19)
p q
[[(] & C(I exp
q pp+ q
+11
(20)
20

---

Published as a conference paper at ICLR 2022
Canfirmod
Recoward
sctuol
900
actual wag/oway
25
actual Exp eogilog
Count (Wymmia))
$50
1050
1000
750
700
650
Sene
Lanc Recoveicd
Deaaltt
8000
actual
actral seq2.00)
actual seguinn
Conta (Nrginia
7500
1600
6000
Cun fimmoo
Recovente
Dicaiicc
65000
2500
Counts (leee))
100000
Pade Exp
60000
Pode Exp woql-wa
Fode Exp sesplooy
90000
$5.00
2000
50000
50000
70000
45000
60000
40000
3000
-14
-5
0
3
-14
-9
-6
-3
0
-14
-9
-6
-3
lime (day)4
Lime (Bas5)
Hise (Jays)
Contlime
Recovesed
Deaths
accual Pad: Fap
64000
Park Hxp
Churt (Abii(A (Abeia)
s002900
1250
Seq2so0
$2000
1000
50000
-2000
7.0
500
58000
250
14
5
14
6
14
Time (day
time (days)
time (dyy
Comnte Nebmman)))
210000
sctual
Pade imp
4000
acoual Pade Exp
207500
100000
205000
50000
2000
acttall Pode Exp s002
titte (daps)
time (days) Recovered
tizse (days)
2.5
actual
50000
2.70
se02500
2.0
1.5
1.0
2.60
0.5
actual
Pode Exp Mog
-50000
actual Pade Uxp
0.0
14
(day
(days)
(day
Figure 8: COVID-19 Forecasting. The Padé Exp prediction and the best non-neural operator scheme from Table (seq2seq) is shown. (Top 3 rows) Prediction for the week 06/12/20 06//18/0 using previous 2 weeks data. This week is the test sample with least MAE averaged across Confirmed, Recovered, and Death counts Padé Exp: or test sample with best averaged prediction. Next, row are US states the order of best-middle-wor's averaged prediction for this test sample. (Bottom rows) Same analysis for the week 03/07/21 03/13/21 which is having largest MAE, or the test sample with worst averaged prediction. Similarly, row 4-6 US states in the order of best-middle-wond averaged prediction for this worst test sample. The missing values are substituted with zeroes which one the reason higher errrr..
21

---

Published as a conference paper at ICLR 2022
Next, by using the product rule, we obtain
(0,,0A(E)Il E
times
(a) jay |21 10
A' ((I
P exp
p q |1| 110
p
(21)
p+g
Similarly,
((C)) & C'(1141)110.00444 q cxp
&
9 1pg
121)
1l,0fl,
(22)
p+
p
1/2
dc ||
11e
dc
2
gl
(23)
;03
00,
Using eq.
(18..
and
23
we obtain the bound for eq.
(66)
as
1/2
exp (
dL
2
(cxp(I
p+g p q 121 lb|2
(24)
a0c
20,
which is used to write the eq. in Theorem T Similarly, the analysis is extended
dy
lim 1 - (F(x; 0, W (0, F(x; F(;;0, W. W)) 1-0
max bERdxd 1 1
S
max deRdxd 1 1-1-1
(A(C))a |A(())
(25)
and
oy
1 lim (F(x; 0, W th) F(x;0, W. (- / |A([))ahl] < |A(4),
max hER |h|
0b
max hER" 1ER%,Ih,_-
(26)
where, Dz. and Dr are diagonal matrices with Iss and 0's. Using eq. 25 and 26 we obtain the eq. and in the Theorem I respectively.
22


Reviews:
[{"review_id": "jtCLeY7aid", "reviewer": "Reviewer_tvEi", "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": "This paper uses a rational polynomial to approximate the exponential function and proposed a Pade network. The idea is reasonable and the experimental result is convincing. Here are several minor points.\n- The Pade approximation need the inverse, the paper use  a non-linear layer which we implement as a simple fully connected layer to approximate the inverse. Can the author justify this approximation, as an example justify the network is powerful enough?\n-exp(x) can also be approximated by (1+(1/n)x)^{n}, can the author justify the benefit of the Pade approximation beyond the simple (1+(1/n)x)^{n} approximation. (the approximation doesn't need an inverse) I guess [1] is doing such approximation.\n- After learning the operators in the Pade approximation, can the learned operator be generalized to other selections of (p,q)? As an example using one set of (p,q) for training and another (p,q) for testing?\n- To justify an operator is learned, an experiment to show the generalization ability cross different grids is needed.\n\n[1] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew\nStuart, and Anima Anandkumar. Markov neural operators for learning chaotic systems, 2021.", "summary_of_the_paper": "This paper proposed a recurrent Pade network fro learning non-linear operator approximations for IVP.  The Pade exponential operator uses a recurrent structure with shared parameters to model the non-linearity compared to recent neural operators that rely on using multiple linear operator layers in succession. The paper showed that Pade network does not suffer from the issue of gradient explosion and the boundedness of the gradients can be established", "summary_of_the_review": "My first main concern is that exp(x) can also be approximated by (1+(1/n)x)^{n}, the paper should take this as an baseline. Second is the nonlinear layer is not powerful enough to approximate the inverse.", "strength_and_weaknesses": "", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "XDCWTVjC3L", "reviewer": "Reviewer_uucn", "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": "The idea in this paper is novel and promising. The approach follows the basics in PDEs, resulting in an efficient method. \n\nThere are multiple issues with this paper that I can not recommend acceptance at its current stage.\nThe first is related to the notation and explanation, and the second is the experimental study. \n\n1) It seems there is room to improve notation and explanation.\n1a) In eq 3, what does the O plus sign do?\n1b)  Is L an operator between function spaces or finite spaces?\nIt is not clarified.\n1c ) the problem construction is missing in the paper. \n\n1d) how from Pade approximation on z in C do we have derivations of A applied on opeators? \n1e) This paragraph \"Pade ́Approximation Givenananalyticfunctionfpzq\" needs clean up. The sentences are not easy to follow.\n\n1f) what are the norms? what are the input-output spaces of operators? what are the spaces?\n1g) In the proof of Theorem, A is applied on z and also L. Is A a spectral function? \n1h) If everything here is finite-dimensional, then a more clarifying notation and problem structure are needed. It seems the authors jump between operator between function spaces and finite-dimensional vectors. The authors are encouraged to be consistent. \n1i) what is F in the proof of the theorem? \n1j) fix eq 19\n1k) This list is longer, and the authors are encouraged to improve their problem construction.\n1l) For what case Pade's approximation is valid? what are the a and b coefficients when dealing with operators?\n2) The empirical study is not extensive compared to the prior works.\n2a) are the numbers in the table and plots for one-step prediction? if yes, what happens if the operators are composed to predict future time steps?\n2b) how the baselines are trained?\n2c) Prior works study 2d Navier Stocks, the authors are encouraged to provide an extended study on that front to provide more compelling results on the benefits of their methods.\n2d) study of further statistics is missing\n3) Why not directly learn the exp(L) instead of using such thing as Pade approximation?\n", "summary_of_the_paper": "This paper study the problem prediction in time-evolving partial differential equations. \nInspired by the nature of solutions in PDEs where the solution often time can be written in terms of exponent of the operator and Pade approximation exponent of operators, the authors propose a recurrent articture to learn solution operator in PDEs.\n\nThe paper proposes a nice idea and approach to solve the mentioned problem. ", "summary_of_the_review": "A great work, and a great idea. \n\nHowever the notation needs a lot of work, the explanation needs a lot of work, the empirical study requires more work and explanation. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made 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": "SlHmOMeCdo", "reviewer": "Reviewer_QGHZ", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "The recently introduced neural operator framework exhibits several advantageous properties, and works that extend and generalize these methods are of interest. In this context, the current paper proposes to model the inherent nonlinearity in the data via exponential operators. From a differential equations viewpoint, this is a reasonable assumption as many IVPs could be solved using exponential maps. The theoretical stability guarantee is a nice addition, especially since RNN are known to be challenging to train. Finally, while the evaluation is not extensive with respect to state-of-the-art baselines and datasets, it does convey the message that the proposed method attains SOTA or potentially even beyond SOTA results.\n\nMy main concern regarding the paper and the reason for my lukewarm score is the unclear novelty in the method. Specifically, the expRNN (https://arxiv.org/abs/1901.08428) and dtriv (https://arxiv.org/abs/1909.09501) architectures considered the exponential map in the context of Lie groups. Essentially, these approaches use the scaling-squaring and a Pade approximant of Al-Mohy & Higham for which an analytic Jacobian can be derived and computed. In particular, no matrix inversion is explicitly calculated (see https://github.com/Lezcano/expRNN/blob/master/expm32.py). Also, this matrix exponential computation is *not* recurrent, and its stability is mostly governed by the modulus of the eigenvalues of the generating operator ($\\mathcal{L}$ in your notation). Given the expRNN/dtriv line of work, it is not clear why the authors design a new recurrent architecture for computing the matrix exponential, and not simply using expm in torch. This approach will eliminate the recurrence, use less weights, will probably be faster, and will actually compute the matrix exponential. Regarding the last bit---since your approach is based on a learning component whose role is to learn the inverse polynomial, it is not entirely clear what is the actual space of operators your method learns. Please clarify.\n\nMinor comments:\n\t- Is $\\tau$ fixed in the operator problem? Why, and why is this problem so significant? Typically, one is interested in the prediction of one or more steps to the future (as in your COVID-19 example)\n\t- Is your $T$ finite- or infinite-dimensional? It seems like it is infinite-dimensional, and thus some discussion is required regarding its approximation as the $A_i, B_i, C_i$ are finite-dimensional.\n\t- I am somewhat puzzled about Eq. (4). Specifically, $T$ appears on the left and right sides of the equation. In addition, some of the operators are not specified ($P_i$?). How is this form used in practice?\n\t- Why $A_i, B_i, C_i$ are modeled as exponential operators?\n\n", "summary_of_the_paper": "This paper introduces a new approach for solving the operator map problem which is based on the non-standard form and its approximation via exponential operators. The authors integrate their method in the recently introduced neural operator framework. The technique is evaluated on two synthetic PDE problems and one real life example. In all cases, the proposed approach improves over the other baselines on the relative L2 metric.", "summary_of_the_review": "The proposed method is interesting as it extends a recently introduced framework (neural operators). However, the main technical novelty seems not necessary in the context of previous work.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "1: The contributions are neither significant nor novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "7A1zz6tERnw", "reviewer": "Reviewer_EkGn", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "The paper proposes a novel and natural model for operator learning. \n\nPros:\n1. The paper writes out the time update as in Table 1.\n2. The paper uses the Pade approximation for the exponential operator.\n3. The paper gets the state-of-art results on both PDEs and real-life benchmarks.\n\nCons:\nIt seems to me that the exponential operator only addresses the linear term exp(tA), while the most challenging and interesting part of the non-linear equation is the time integral term addressing the non-linear term N. This design may have a limitation on harder non-linear equations such as the Navier-Stokes equation. Would it be possible to add an additional term to approximate the time integral?", "summary_of_the_paper": "This paper proposes to learn the exponential operator in time evolution using the Pade approximation. The proposed method shows the states of art results on 1d Burgers, KdV equation, as well as the covid 19 benchmark.", "summary_of_the_review": "In my opinion, this work is novel and meaningful. However, it is less satisfactory not to address the more interesting non-linear term N. I recommend borderline acceptance.", "strength_and_weaknesses": "", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper provides a unique contribution that uses Padde approximations to approximate non-linear operators for solving initial value problems in PDEs. The paper contains also a non-trivial experiment with a real-world dataset that showcase the impact of the proposed model. The authors have provided a strong rebuttal and therefore I recommend Accept.
Decision: Accept (Poster)