Year: 2021

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 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 (53% (52%) better MAE than best neural operator (non-neural operator deep learning model)) compared to state-of-the-art forecasting models.

heading: 1 INTRODUCTION
text: Predicting the future states using the current conditions is a fundamental problem in machine learning. Such problems fall under the umbrella of a 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 (Lovelace, 2021)), neuroscience (Hodgkin-Huxley model (Zhang et al., 2020)), engineering (fluid dynamics (Wendt, 2008)), water waves (tsunami (Elbanna et al., 2021)), mean field games (Ruthotto et al., 2020), to list just a few. Within the current pandemic context, the applications areas like epidemiology (Kermack–McKendrick model (Kermack et al., 1991; Diekmann et al., 2021)) are of tremendous interest.
Neural Operators The use of 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 are suitable to solve problems with a large number of training samples. Moreover, because of 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 differential equations (PDEs) versions of the IVPs for a given instance. The works of (Kochkov et al., 2021) 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 a 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;c;a; Gupta et al., 2021; Bhattacharya et al., 2020; Patel et al., 2021). Most of these approaches tried to efficiently work with the integral kernel operators, for example, Graph Nyström sampling in (Li et al., 2020b), convolution approximation in (Li et al.,
2020a), multiwavelets compression in (Gupta et al., 2021). 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-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 as a solution of: dydt “ at, yp0q “ y0 as yptq “ e
aty0 (for more general examples, see Table 1). 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).
Padé 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-squaring (SSQ) methods are commonly suggested approaches to deal with the errors (Lawson, 1967). In addition to SSQ, the Padé 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 Padé approximation based approach to compute the matrix exponential expm (Al-Mohy & Higham, 2009). Matrix exponential via Padé representation requires dense matrix computations (for example, inverse and higher-order polynomials). Such operations are not numerically feasible, 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 Padé approximation.
Our strategy, in this work, is to 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 3.3). To the advantage of Padé approximation, the exponential of a 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 Padé approximation, we compute the exponential of the operator using a novel recurrent neural architecture that also eliminates the need for matrix inversion. (iii) We theoretically demonstrate that the proposed recurrent scheme, using the Padé coefficients, have bounded gradients with respect to (w.r.t.) the model parameters across the recurrent horizon. (iv) We demonstrate the data-efficiency on the synthetic 1D datasets of Korteweg-de Vries (KdV) and Kuramoto–Sivashinsky (KS) equations, where with less parameters we achieve state-of-the-art performance. (v) We formulate and investigate the epidemic forecasting as a 2D time-varying neural operator problem, and show that for real-world noisy and scarce data, the proposed model outperforms the best neural operator architectures by 53% and best non-neural operator schemes by 52%.

heading: 2 OPERATORS FOR INITIAL VALUE PROBLEM
text: 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 for space-discretization. Section 2.3 describes the proposed use of canonical exponential operators and presents a novel architecture using Padé approximation.

heading: 2.1 INITIAL VALUE PROBLEM
text: The initial value problem (IVP) for PDEs can be written in its general form as follows.
ut “ Fpt, uq, x P Ω upx, 0q “ u0pxq, x P Ω
(1)
where, ut is the first-order time derivative of u, F is a time-varying differential operator (non-linear in-general) such that F : R` Y t0u ˆ B Ñ B with B being a Banach space. Usually, the system in eq. (1) is required to satisfy a boundary condition such that Bupx, tq “ 0, x P BΩ @t in the solution horizon, and BΩ is the boundary of the computational region Ω 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 U as two Sobolev spaces Hs,p with s ą 0, p “ 2, an operator T is such that T : A Ñ U . For a given τ ą 0 and two functions u0pxq and upx, τq, in this work, we take the operator map as T u0pxq “ upx, τq with x P Ω. Table 1 summarizes 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., 2021) modeled the non-linear operators in one way or another by using multiple canonical integral operators along with non-linearity (for example, ReLU). In this work, we directly produce an exponential operator approximation. First, we discuss an efficient basis (multiwavelets) for space discretization of the input/output functions in Section 2.2.

heading: 2.2 MULTI-RESOLUTION ANALYSIS
text: The multi-resolution analysis (MRA) aims at projecting a 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 a finite interval (Alpert et al., 2002). The multiwavelets are 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, 1997). However, we do not rely on this assumption in this work. Here, we briefly introduce the MRA and refer the reader to Gupta et al. (2021) for a detailed picture.
Notation We begin by defining the space of finite interval polynomials as Vkn “ tf |f are polynomials of degree ă k defined over interval p2´nl, 2´npl ` 1qq for all l “ 0, 1, . . . , 2n ´ 1, and assumes 0 elsewhereu. The Vkn are contained in each other for subsequent n or,
Vk0 Ă Vk1 Ă . . . Ă Vkn´1 Ă Vkn Ă . . . . (2)
1Time-advection equation with linear operators L,N and non-linear function fp.q. A wide range of problems can be modeled, for example, Korteweg-de Vries, Kuramoto-Sivashinsky, Burgers’ Equation, Navier-Stokes (list not exhaustive).
2A non-linear integro-differential solution to the time-advection equation using semi-group approach (Beylkin & Keiser, 1997; Pazy, 1983; Yoshida, 1980). A slightly general version is discussed in (Beylkin et al., 1998).
The orthogonal component of these polynomial spaces is termed as multiwavelet space Wkn and are defined such that
Vkn à Wkn “ Vkn`1, Vkn K Wkn. (3) The orthonormal basis of Vk0 are OPs ϕ0, ϕ1, . . . , ϕk´1 and we have used appropriately normalized shifted Legendre Polynomials in this work. The basis for Vkn and W k n are ϕ n jlpxq “ 2n{2ϕjp2nx´ lq and ψnjlpxq “ 2n{2ψjp2nx´ lq, respectively, for l “ 0, 1, . . . , 2n ´ 1 and j “ 0, 1, . . . , k ´ 1.
Finally, an important trick for representing the operator T in the multiwavelet basis is called nonstandard (NS) form (Beylkin et al., 1991). The NS form decouples the interactions of the scales and is useful in obtaining an efficient numerical procedure. Using NS form, the projection of operator T is expanded using a telescopic sum as follows.
Tn “ ÿn
i“L`1 pQiTQi `QiTPi´1 ` Pi´1TQiq ` PLTPL, (4)
where, Pn : Hs,2 Ñ Vkn is the projection operator, Tn “ PnTPn, Qn : Hs,2 Ñ Wkn such that Qn “ Pn ´ Pn´1, and L is the coarsest scale under consideration pL ě 0q. Therefore, the NS form of the operator is a collection of the triplets tAi, Bi, Ciuni“L`1 and PLTPL with Ai “ QiTQi, Bi “ QiTPi´1 and Ci “ Pi´1TQi. In this work, we aim to model Ai, Bi, Ci as the exponential operators to better learn the IVP by explicitly embedding the non-linear operators into the multiwavelet transformation. This is not straightforward due to 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.

heading: 2.3 EXPONENTIAL OPERATOR APPROXIMATIONS
text: Due to the nature of first-order time-evolution equations, the exponential operators appear in the solution of IVP as discussed in Section 2.1. Being an analytic function, the exponential also assumes a Taylor series expansion. However, the approximation error by truncation is (Abramowitz & Stegun, 1965) and may require a large number of coefficients. We now discuss a better approximation for the non-linear functions.
Padé Approximation Given an analytic function fpzq, z P C at 0 and let p, q P N, the rp{qs Padé approximation of f at 0 is a rational polynomial rpqpzq “ Apqpzq{Bpqpzq, whereApqpzq, Bpqpzq are polynomials of degree p, q, respectively. For our work, fpzq “ ez , and the rp{qs Padé approximation to the exponential function is written as Apqpzq “ řp j“0 ajz j and Bpqpzq “ řq j“0 bjz j , where
aj “ pp` q ´ jq!p! pp` qq!j!pp´ jq! , 0 ď j ď p bj “ pp` q ´ jq!q!p´1qj pp` qq!j!pq ´ jq! , 0 ď j ď q. (5)
Note that, a0 “ b0 “ 1, and we now discuss the exponential of operators. Given an operator L (possibly non-linear), the rp{qs Padé approximation for eL can be written as
eL « rpqpLq “
¨
˝
q ÿ j“0 bj L ˝ L ˝ . . . ˝ Llooooooomooooooon
j´times
˛
‚
´1 ¨
˝
p ÿ j“0 aj L ˝ L ˝ . . . ˝ Llooooooomooooooon
j´times
˛
‚. (6)
For L being a linear transformation, the inverse in eq. (6) is computed as a matrix inverse, for example when evaluating a matrix exponential. Even convolution operator can be represented as the circulant matrix (shift-invariant kernel) operation. But storing the big circulant matrix is not a numerically efficient solution. Instead, we now discuss a neural architecture that emulates the Padé approximation for operators in eq. (6) while avoiding taking the inverse.
Padé Exponential Model We compute the operator exponential using the rp{qs Padé approximation via a recurrent neural architecture in Figure 1. First, the denominator polynomial BpqpLq is evaluated using the left recurrent network. Second, the output is passed through a non-linear layer which we implement as v “ σpWu` bq with σp.q 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-linear layer is passed through another recurrent network implementing the numerator polynomial ApqpLq. Note that, both polynomial recurrent networks use same operator L, or in other words, the parameters are shared across the network. With aj , bj fixed-beforehand, the total trainable parameters of the Padé exponential model are tθL,W, bu.
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 a 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 rp{qs Padé 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 “ LpθLq, a non-linearity layer v “ σpWu`bq, and p, q P N, at points of differentiability, the gradients of the operation x ÞÑ y “ F px; θL,W, bq :“ rp{qseLpxq using the rp{qs Padé network in Figure 1 are bounded in operator norm by∥∥∥∥ ByBθL ∥∥∥∥ ď expp}L}q p}b}2 ` }W }}x}2q ˜ nθ ÿ
j“1
∥∥∥∥ BLBθj ∥∥∥∥2 ¸1{2
, (7)∥∥∥∥ ByBW ∥∥∥∥ ď expp}L}q}x}2, (8)∥∥∥∥ByBb ∥∥∥∥ ď expˆ pp` q }L} ˙ . (9)
The detailed proof is provided in the Supplementary Materials (Appendix E). Note that, nθ “ |θL| is such that nθ ! 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 σpWu ` bq 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 ∥x∥2 which could be pre-normalized.
Finally, the Padé neural model is integrated with the multiwavelet transform by substituting each of A,B,C, from Section 2.2, with rp{qseLA , rp{qseLA , rp{qseLA . The complete flow-diagram by plugging-in the Padé model is shown in Appendix D.

heading: 3 EXPERIMENTS
text: Here, we empirically evaluate the proposed model and compare against the existing approaches. We consider several synthetic PDE datasets as 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 u0pxq and the output function upτ, xq for some
pre-specified τ ą 0, are evaluated at M discretized locations of the domain Ω. This yields a single training sample of pu0pxiq, upτ, xiq, xi P Ω, 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. Padé Model Implementation The operator L in Figure 1 is fixed as a single-layered convolution operator for 1D datasets, and 2-layered convolution for 2D datasets. For getting the input/output operator mapping, the multiwavelet transform is used only for discretizing the spatial domain. The Padé neural model easily fits into the sockets of the multiwavelet transformation based neural operator as sketched in Figure 7 (Appendix D). The multiwavelet filters are obtained using shifted Legendre OPs with degree k “ 4. In contrast to the work of Gupta et al. (2021), only a single cell of multiwavelet transform is used in the current work because the non-linearity in the operator is explicit via Padé neural operator. This saves a lot of trainable parameters and yields a 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 Padé approximation is fixed as pp, qq “ p5, 6q for 1D models and p4, 2q for 2D models. While the authors do not advocate that this is the best possible choice, an ablation study for 1D data is performed in Appendix B.1.
Benchmark Neural Operators We compare against the recently proposed neural operator works with the proposed Padé exponential model (Padé 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 et al., 2020c). LNO A 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 Padé Exp model delivers state-of-the-art (Sota) performance on a range of datasets, both synthetic and real-world. With less number of parameters, the proposed model shows promise for small datasets as shown in Section 3.3.
Training parameters All neural operator models are trained using Adam optimizer with a learning rate of 0.001 and decay of 0.95 after every 100 steps. The loss function is taken as the relative L2 error. For synthetic datasets we train for a total of 500 epochs and for real-world COVID-19 dataset we train for a total of 750 epochs. All experiments are done on an Nvidia A100 40GB GPUs.

heading: 3.1 KORTEWEG-DE VRIES EQUATION
text: The Korteweg-de Vries (KdV) equation is a one-dimensional non-linear PDE used to model the non-linear shallow water waves. For a given field upx, tq, the KdV PDE takes the following form:
ut “ ´0.5u Bu Bx ´ B3u Bx3 , x P p0, 1q, t P p0, 1s
u0pxq “ upx, t “ 0q. (10)
For Table 1 notations, we have L “ ´ B 3 Bx3 , fpuq “ ´0.25u 2, and N “ BBx . The neural operator learns the mapping of the initial condition u0pxq to the solutions upx, t “ 1q. The initial condition
is generated in Gaussian random fields according to u0 „ N p0, 74p´∆ ` 72Iq´2.5q 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 210, and datasets with lower resolutions are obtained by sub-sampling the highest resolution dataset.
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 sampling the models are evaluated. For consistency, each model is fed the same training sub-sample and the results are shown in Figure 2 (Left). We see that Padé Exp model has the sharpest decay compared to the other state-of-the-art neural operators on varying N . This shows that the proposed model is data-efficient and works well when less data is available. Next, we also evaluate the Padé Exp model on the complete data N “ 1000 but by varying the input resolution as shown in Table 2. The proposed model performs consistently better for all the input resolutions.

heading: 3.2 KURAMOTO-SIVASHINSKY (KS) EQUATION
text: The Kuramoto-Sivashinsky (KS) equation is a fourth-order non-linear PDE derived to model the diffusive instabilities in a laminar flame front. For a given field upx, tq, the KS PDE takes the following form:
ut “ ´u Bu Bx ´ B2u Bx2 ´ B4u Bx4 , x P p0, 1q, t P p0, 1s
u0pxq “ upx, t “ 0q. (11)
The KS equation is also time-advection and according to the Table 1 notations, we have: L “ ´ B 2
Bx2 ´ B4 Bx4 , fpuq “ ´0.5u 2, and N “ BBx . The neural operator learns the mapping of the initial condition u0pxq to the solutions upx, t “ 1q. Similarly to the KdV equation in Section 3.1, the initial condition is sampled from a Gaussian random field u0 „ N p0, 74p´∆ ` 72Iq´2.5q 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 in a similar experimental setup as for the KdV equation in Section 3.1. The reduced training experiment results are shown in Figure 2 (Right), where we have uniformly and randomly sub-sampled the training samples from the complete dataset. For consistency, we have evaluated all models on the same sub-sampled training set. We again observe that the proposed Padé Exp model with its compact
structure has the steepest decay of relative L2 error compared to the recent neural operator works. By using the exponential operators, the IVPs can be efficiently solved as we witness from Figure 2. Finally, we also evaluate the proposed model for the complete data but varying the input resolution s in Table 3. We again conclude that the proposed model attains better performance than state-of-the-art approaches for all the resolutions.

heading: 3.3 EPIDEMIC FORECASTING: COVID-19 STUDY
text: The epidemic forecasting problem refers to the prediction of future counts of infected individuals, recovered individuals and deaths using the current observational data. A variety of compartmental models exist to model the epidemic spread with certain assumptions (see related work, Appendix A). The dynamics of the epidemic spread by modeling their time-evolution behavior is not exactly known, and may not apply to any problem, for example, the recent COVID-19 pandemic. Neural operators, providing a complete data-driven approach, are capable to learn PDE agnostic maps. Consequently, we show that the epidemic forecasting can be formulated as an operator map learning problem.
Dataset The COVID-19 data set 3 from April 12th 2020 to August 28th 2021 is provided by Johns Hopkins University (Dong 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 ˆ 3. Deep learning Benchmarks In addition to the neural operators, we also compare against the stateof-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 et al., 2017) utilizing encoder-decoder structure. For Seq2Seq, we have used LSTM architecture for encoder and decoder. Neural ODE (Chen et al., 2019) 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 7 days data for each of the 50 US states. Let dt be the 50 ˆ 3 array on day-t, then the operator map can be written as follows.
T pd´14, d´13, . . . , d´1 looooooooooomooooooooooon
u0pxq
q “ pd0, d1, . . . , d6 loooooomoooooon
upτ,xq
q.
Forecasting The COVID19 forecasting benchmarks are presented in Table 4. Due to data scarcity (484 samples in total), we do a 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 Padé 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 Padé 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 3. We also show the corner cases (best and worst test sample) predictions in the Figure 8 (Appendix B.4) along with the best/worst prediction states.

heading: 4 ABLATION EXPERIMENTS
text: The ablation experiments are presented in the Supplementary Materials (Appendix B.1). The following experiments are performed; (i) Comparison between the proposed Padé approximation vs Taylor Series neural operator, (ii) the variations of p, q for the rp{qs Padé approximation model in Figure 1, and (iii) variation of the non-linearity module in the Padé neural model.

heading: 5 FUTURE DIRECTIONS AND CONCLUSION
text: Time-evolution analysis of PDEs and forecasting the future states from a set of observations is fundamental for studying a 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 a 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 Padé approximation based exponential neural operator architecture. The proposed model has order-magnitude fewer parameters (see Table 4) and attains data-efficiency.
The IVPs deal with first-order time evolution, and thus, the exponential operator appears. In a more general setup, the higher-order time derivatives, or even fractional time derivatives should be considered to model the non-Markovian 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.
3https://github.com/CSSEGISandData/COVID-19

heading: A RELATED WORK
text: Epidemic Forecasting In Epidemiology, the compartmental models are very common, where the population is assigned into discrete categories (for example, Susceptible, Infected, Recovered), also termed as SIR models which date back to 1927 (Kermack et al., 1991; Diekmann et al., 2021). Along the lines, several works exist that model the time-evolution of epidemic using various assumptions. The work of Chen et al. (2020) uses a ridge regression based prediction. An extension of SEIR (E for Exposed) compartments to sub-populations moving across different places in (Pei & Shaman, 2020). A convolution filter based approach in (Wang et 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 et al., 2020b) incorporating the unreported cases. The work of Davis et al. (2020) which simultaneously models the transmission rate among US states. A physics-based model in (Wang et al., 2021a) that uses auto-differentiation using Runga-Kutta methods to estimate the model parameters.

heading: B ADDITIONAL RESULTS
text: B.1 ABLATION STUDY
Taylor vs Padé Exponential Approximation We compare the Padé approximation based neural model in Figure 1 vs Taylor series approximation for evaluating the exponential of the operator L. The Taylor series version is written as
eL « Lt ÿ
j“0
1 j! pL ˝ L ˝ . . . ˝ Lq loooooooomoooooooon
j´times
,
where, Lt is the truncated length of the series. We evaluate the truncated Taylor series approximation vs rp{qs Padé approximation in Figure 1 for same number of maximum polynomial degree, i.e., Lt “ maxpp, qq. 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 3.1). The fluctuated
input is sampled from Random fields Filip et al. (2019) with controllable parameter λ with lower values attaining higher fluctuations and the input u0pxq reaches the Brownian motion limit with λ Ñ 0. We see in Table 5 that the Padé approximation model works better compared to Taylor series based expansion for incorporating the exponential operators for each λ.
L’Hospital Approximation The exponential exppxq can also be approximated as limnÑ8p1`x{nqn using L’Hospital rule (Krantz, 2004). For a given operator L, the L’Hospital based exponential operator approximation can be written as
eL « pI ` L{nq ˝ pI ` L{nq ˝ . . . ˝ pI ` L{nq looooooooooooooooooooooooomooooooooooooooooooooooooon
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 Padé model in Table 5 for the same degree of operator polynomials, i.e., n “ maxpp, qq. For the same experimental settings as taken in the Taylor approximation comparison, we see that the proposed Padé model performs better than the L’Hospital-based exponential approximation for all fluctuations strengths.
Varying Padé polynomials degree The numerator/denominator polynomials degree p, q in eq. (5) for rp{qs Padé approximation model in Figure 1 are hyper-parameters. We vary the p, q over a 2D grid of r1, 10s ˆ r1, 10s and evaluate for the KdV equation settings (s=1024) as mentioned in the Section 3.1 and the results are shown in Figure 4. We make the following observations, namely; (i) For p “ q “ 1, i.e., when there is no operator L and only the non-linear layer σpWu` bq, the model has the worst performance of 0.1063. This shows the importance of operator L in the Padé model. (ii) By increasing the p, q in either direction, the performance improves as we attain better approximation with higher degree polynomials. The best performance is achieved as pp, qq “ p10, 1q. (iii) The performance is roughly similar in the mid-range of p and q’s, as long as p ą 1 and q ą 3 for the 1D experiment on KdV. Note that, although larger values of p and q are preferable but they also incur the additional cost of run-time which scales linearly with p` q as evident from the recurrent structure in Figure 1.
Structure of the non-linear MLP We have used the non-linear layer σpWu` bq in the Padé neural model for all the experiments with σp.q=ReLU. Now, we vary the depth of this layer from 1 to 3 (each layer is followed by ReLU non-linearity) to see the individual contribution of this layer towards the final test performance. The experiments settings are same as for Table 2 and we fix s “ 1024. Upon increasing the depth we do not see any improvements, in fact, after having three layers, the performance starts to degrade.
A possible reason could be Denominator polynomial getting affected by the deeper non-linearity layer.
B.2 PREDICTION AT HIGHER RESOLUTIONS
B.3 TRAINING/EVALUATION WITH DIFFERENT pp, qq
Due to the recurrent structure of the Padé approximation model in Figure 1, with only trainable parameters tθL,W, bu, the model can be trained and tested on different values of pp, qq by just varying the recurrence loop lengths of the numerator and denominator polynomial. We setup this experiment on the KdV equation (same setting as in Section 3.1). The Padé model is trained for pp, qq “ p5, 6q and p4, 2q and then tested on r1, 10s ˆ r1, 10s grid of p, q tuples as shown in Figure 5.
We observe that the best performance (lowest relative L2 error) results when the train and test pp, qq matches. In addition, we observe that for a fixed q, the best performance results when p is such that p/q ratio is similar as the training ratio. For example, in Figure 5 (Left), we see that lowest value for each column is obtained when p{q is closest to 5{6 while in (Right) we see lowest value of each column when p{q “ 2. The boundary cases diverge because either p, q “ 1 does not have any operator L in the corresponding polynomial as evident from eq. (6). When trained for pp, qq “ p4, 2q, the best performance for q ą 5 saturates to p “ 10 because of the maximum limit on the experimental grid pp “ 10q. Overall, this experiment suggests that the relative length of the polynomials during evaluation should be similar to training setting for best performance.
B.4 COVID-19 FORECASTING
Additional forecasting results for COVID-19 by considering the corner cases are presented in Figure 8.
B.5 2D NAVIER-STOKES EQUATIONS
The Navier-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-art models (Li et al., 2020a; Gupta et al., 2021) under the same conditions, we use the same data sets that have been published in (Li et al., 2020a), where the NS equations take the following form:
wtpx, tq ` upx, tq ¨ ∇wpx, tq ´ ν∆wpx, tq “ fpxq, x P p0, 1q2, t P p0, T s ∇ ¨ upx, tq “ 0, x P p0, 1q2, t P r0, T s
w0pxq “ wpx, t “ 0q, x P p0, 1q2 (12)
where u is the velocity, w is the vorticity such that w “ ∇ ˆ u. The incompressible flow is modeled via divergence condition as ∇ ¨ upx, tq “ 0. We set the experiments to let the neural operator map the first 10 time units to last T´10 time units of vorticityw. The initial condition is generated in Gaussian random fields according to w0 „ N p0, 7 3 2 p´∆ ` 72Iq´2.5q with periodic boundary conditions and the forcing function is fpxq “ 0.1psinp2πpx1 ` x2qq ` cosp2πpx1 ` x2qqq. We experiment with different viscosities ν, final time T , and the number of training pairs N : piq ν “ 1e ´ 3, T “ 50, N “ 1000; piiq ν “ 1e ´ 4, T “ 30, N “ 1000; piiiq ν “ 1e ´ 4, T “ 30, N “ 10000; pivq ν “ 1e´ 5, T “ 20, N “ 1000 on a 64 ˆ 64 grid. The time-varying 2D data is modeled as 3D operator learning problem. We implemented the Padé exponential model by taking L as 2-layered 3D CNNs with ReLU non-linearity and p “ q “ 4 in
Figure 1. A total of 4 layers of multiwavelet skeleton (Figure 7) with k “ 3 are concatenated using ReLU nonlinearity. The results are reported in Table 8, and we observe that Padé model performs well compared to the recent state-of-the-art neural operator approaches. For the less data setup of N “ 1000 with ν “ 1e ´ 4, we see in Figure 6 that the Padé Exp model quickly converges to the lowest value compared to the other neural operator approaches.
B.6 NUMERICAL VALUES
The numerical values for the training samples variation experiment for the KdV and KS equation in the Figure 2 are shown in Table 9 and Table 10, respectively.
C NOTATIONS Operator Learning
T,L Operators between function spaces Hs,p Sobolev spaces such that constituent functions and their weak
derivatives upto order s have finite Lp norms Multiwavelets
à
Subspace addition
Vkn tf |f are polynomials of degree ă k defined over interval p2´nl, 2´npl ` 1qq for all l “ 0, 1, . . . , 2n ´ 1, and assumes 0 elsewhereu
Wkn Orthogonal space to V k n such that W k n
À
Vkn “ Vkn`1 Pn Projection operator such that Pn : Hs,2 Ñ Vkn Qn Projection operator such that Qn : Hs,2 Ñ Wkn L Coarsest scale of the multiwavelet transform
Padé Approximation
}.} Operator norm }.}2 Euclidean L-2 norm rp{qsf Padé approximation of function f with numerator/denominator
polynomial degree p/q, respectively
θL Parameters used to represent the operator L θn Number of parameters used to represent the operator L x ÞÑ y “ F px; Θq Mapping from x to y using function F with parameters set Θ

heading: D MULTIWAVELET EXPONENTIAL OPERATOR ARCHITECTURE
text: The Padé neural operator based multiwavelet transform model is shown in Figure 7. The input to the model is spn`1q, where n “ logpMq, and the output is U pnqs , where n is the finest scale (or the log of input resolution). For a detailed description of multiwavelet transform, we refer the reader to (Gupta et al., 2021).

heading: E PROOF OF THEOREM 1
text: First note that, for the Padé approximation rp{qsex, the polynomial coefficients in eq. (5), we have
aj “ 1
j! ppp´ 1q ¨ ¨ ¨ pp´ j ` 1q pp` qqpp` q ´ 1q ¨ ¨ ¨ pp` q ´ j ` 1q ď 1 j!
ˆ
p
p` q
˙j
. (13)
and similarly,
|bj | ď 1
j!
ˆ
q
p` q
˙j
.
For the simplicity of notation: let Apzq “ Apqpzq “ řp j“0 ajz j , A1pzq “ řp j“1 jajz j´1, Cpzq “ řq
j“0 |bj |zj , C 1pzq “ řq j“1 j|b|jzj´1. Note |Bpqpzq| ď Cp|z|q and |B1pqpzq| ď C 1p|z|q. Using eq. (13), we can further write that
|Apzq| ď exp ˆ p p` q |z| ˙ , |A1pzq| ď p p` q exp ˆ p p` q |z| ˙
(14)
|Cpzq| ď exp ˆ q p` q |z| ˙ , |C 1pzq| ď q p` q exp ˆ q p` q |z| ˙
(15)
Further let d and nθ denote the dimensions of x and θ respectively: x, y, u, v P Rd and θ :“ θL P Rnθ . For the operation x ÞÑ y, we denote y “ F px; θL,W, bq :“ rp{qseLpxq. Then y “ ApLqv, v “ σpWu` bq, u “ BpLqx, and∥∥∥∥ ByBθL ∥∥∥∥ “ maxgPRnθ ,}g}2“1 ∥∥∥∥limtÑ0 1t pF px; θ ` tg,W, bq ´ F px; θ,W, bqq ∥∥∥∥ 2
“ max g P Rnθ , }g}2 “ 1 h P Rd, }h}2 “ 1 |x∇θxy, hy, gy| (16)
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 Bg,θ “ řnθ j“1 gj B Bθj , the differential operator induced by g. The quantity of interest is |Bg,θxy, hy|. Next, by chain rule Bg,θxy, hy “ xpBg,θApLqqv, hy ` xApLqBg,θv, hy
“ xpBg,θApLqqv, hy ` xApLqDWBg,θBpLqx, hy (17) where Dg is a diagonal matrix of 1’s and 0’s depending upon the signatures of entries in Wu` b (see for instance the computation of Lipschitz constants in Zou et al. (2020a)). Hence }Dg} “ 1. Since σ (ReLU) is contractive with Lipschitz constant 1, }v}2 ď }W }}u}2 ` }b}2 ď }W }}BpLq}}x}2 ` }b}2. Substituting back in eq. (17) we obtain
|Bg,θxy, hy| ď }Bg,θApLq}}W }}BpLq}}x}2 ` }Bg,θApLq}}b}2 ` }ApLq}}W }}Bg,θBpLq}}x}2 (18)
Using eq. (14)-(15), the spectral norms are bounded further by
}ApLq} ď Ap}L}q ď exp ˆ p p` q }L} ˙ , (19)
}BpLq} ď Cp}L}q ď exp ˆ q p` q }L} ˙ . (20)
Next, by using the product rule, we obtain
}Bg,θApLq} ď p ÿ
j“1 |aj |}Bg,θpL ˝ L ˝ . . . ˝ Llooooooomooooooon
j´times
q}
paq ď
p ÿ j“1 jaj}L}j´1}Bg,θL}
“ A1p}L}q}Bg,θL}
ď p p` q exp
ˆ
p p` q }L} ˙ }Bg,θL}. (21)
Similarly,
}Bg,θBpLq} ď C 1p}L}q}Bg,θL} ď q p` q exp ˆ q p` q }L} ˙ }Bg,θL}, (22)
}Bg,θL} ď nθ ÿ
j“1 |gj |} BL Bθj } ď }g}2
˜
nθ ÿ
j“1 } BLBθj
} 2 ¸1{2 . (23)
Using eq. (18), (22), and (23), we obtain the bound for eq. (16) as∥∥∥∥ ByBθL ∥∥∥∥ ď ˆexpp}L}q}W }}x}2 ` pp` q exp ˆ p p` q }L} ˙ }b}2 ˙ ˜ nθ ÿ j“1 } BLBθj } 2 ¸1{2 , (24)
which is used to write the eq. (7) in Theorem 1. Similarly, the analysis is extended as∥∥∥∥ ByBW ∥∥∥∥ “ max
ΦPRdˆd,}Φ}2“1 ∥∥∥∥limtÑ0 1t pF px; θ,W ` tΦ, bq ´ F px; θ,W, bqq ∥∥∥∥ 2
“ max ΦPRdˆd,}Φ}2“1 }ApLqDΦΦBpLqx}2 ď }ApLq}}BpLq}}x}2, (25)
and ∥∥∥∥ByBb ∥∥∥∥ “ max
hPRd,}h}2“1 ∥∥∥∥limtÑ0 1t pF px; θ,W, b` thq ´ F px; θ,W, bqq ∥∥∥∥ 2
“ max hPRd,}h}2“1 }ApLqDhh}2 ď }ApLq}, (26)
where, DΦ and Dh are diagonal matrices with 1’s and 0’s. Using eq. (25) and (26), we obtain the eq. (8) and (9) in the Theorem 1, respectively.

Review summaries:
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

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.

This paper study the problem prediction in time-evolving partial differential equations.  Inspired 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.  The paper proposes a nice idea and approach to solve the mentioned problem.