# NON-LINEAR OPERATOR APPROXIMATIONS FOR INI### TIAL VALUE PROBLEMS

**Gaurav Gupta** [1], Xiongye Xiao [1], Radu Balan [2], Paul Bogdan [1]

1
Ming Hsieh Department of Electrical and Computer Engineering
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, pbogdan}@usc.edu, rvbalan@math.umd.edu,_

ABSTRACT

Time-evolution of 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 Pade 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 real-world datasets. 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. We show theoretically that the gradients associated
with the recurrent Pade 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, COVID19) can be formulated as a 2D time-varying operator problem. The proposed Pade´
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.

1 INTRODUCTION

Predicting the future states using the current conditions is a 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 a common term, the “Initial Value
Problems” (IVPs). The basic structure of IVP involves a first-order time-evolution along with nonlinear 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; Bogdan & Marculescu, 2011), 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
Nystrom sampling in (Li et al., 2020b), convolution approximation in (Li et al., 2020a), or the¨
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_

With applications ranging from control systems theory (converting continuous to discrete systems)[“][ at, y][p][0][q “][ y][0][ as][ y][p][t][q “][ e][at][y][0][ (for more general examples, see Table 1).]
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
Pad´e 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


-----

|Equation|Solution|
|---|---|


|d du t “ Au, upt “ 0q “ u 0 (Linear ODE) u t “ B B2 xu 2, upx, 0q “ u 0pxq (Heat equation) u t “ Lu ` Nfpuq, upx, 0q “ u 0pxq 1|upt “ τq “ etAu 0 x2 upx, τq “ eτ BB 2 u 0pxq τ upx, τq “ eτLu 0pxq `şepτ´tqLNfpupx, tqqdt 2 0|
|---|---|


Table 1: 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 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%.

2 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
for space-discretization. Section 2.3 describes the proposed use of canonical exponential operators
and presents a novel architecture using Pad´e approximation.

2.1 INITIAL VALUE PROBLEM

The initial value problem (IVP) for PDEs can be written in its general form as follows.

_ut_ _t, u_ _,_ _x_ Ω
“ Fp q P (1)
_u_ _x, 0_ _u0_ _x_ _,_ _x_ Ω
p q “ p q P

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 H[s,p]** with s ą 0, p “ 2, an operator T is
such that T : . For a given τ 0 and two functions u0 _x_ and u _x, τ_, in this work, we take
_A Ñ U_ ą p q p q
the operator map as T u0 _x_ _u_ _x, τ_ with x Ω.
p q “ p q 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.

2.2 MULTI-RESOLUTION ANALYSIS

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

1Time-advection equation with linear operators L, N and non-linear function f p.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).


-----

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 more detailed discussion.

**Notation We begin by defining the space of finite interval polynomials as Vn[k]**
als of degree ă k defined over interval p2[´][n]l, 2[´][n]pl ` 1qq for all l “ 0, 1, . . .,[“ t] 2[n] _[f]´[|] 1[f][ are polynomi-], and assumes_
0 elsewhereu. The Vn[k] [are contained in each other for subsequent][ n][ or,]

**V0[k]** [Ă][ V]1[k] [Ă][ . . .][ Ă][ V]n[k]´1 [Ă][ V]n[k] [Ă][ . . . .] (2)

The orthogonal component of these polynomial spaces is termed as multiwavelet space Wn[k] [and are]
defined such that
**Vn[k]** **Wn[k]** _n_ 1[,] **Vn[k]** _n[.]_ (3)

[“][ V][k]` [K][ W][k]

The orthonormal basis of V0[k] [are OPs]à[ ϕ][0][, ϕ][1][, . . ., ϕ][k][´][1] [and we have used appropriately normalized]
shifted Legendre Polynomials in this work. The basis for Vn[k] [and][ W]n[k] [are][ ϕ]jl[n] [p][x][q “][ 2][n][{][2][ϕ][j][p][2][n][x][ ´][ l][q]
and ψjl[n] [p][x][q “][ 2][n][{][2][ψ][j][p][2][n][x][ ´][ l][q][, respectively, for][ l][ “][ 0][,][ 1][, . . .,][ 2][n][ ´][ 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.

_n_
_Tn_ (4)
“ _i“L`1[p][Q][i][TQ][i][ `][ Q][i][TP][i][´][1][ `][ P][i][´][1][TQ][i][q `][ P][L][TP][L][,]_
ÿ

where,thatthe NS form of the operator is a collection of the triplets Qn P “n : P Hn ´[s,][2] PÑn´ V1, andn[k] [is the projection operator,] L is the coarsest scale under consideration[ T][n] [“][ P]Ai[n], B[TP]i, C[n][,][ Q]i _i_ _[n] pL[:]L[ H]1 ě[and][s,] 0[2][ Ñ][ P]q. Therefore,[L][ W][TP]n[L][k]_ [ with][such]
t u[n]“ `
_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.

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 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 f pzq, z P C at 0 and let p, q P N, the rp{qs
Pade approximation of´ _f at 0 is a rational polynomial rpq_ _z_ _Apq_ _z_ _Bpq_ _z_, where Apq _z_ and
p q “ p q{ p q p q
_Bpq_ _z_ are polynomials of degree p and q, respectively. For our work, f _z_ _e[z], and the_ _p_ _q_
p q p q “ r { s
Pade approximation to the exponential function is written as´ _Apq_ _z_ _j_ 0 _[a][j][z][j][ and][ B][pq][p][z][q “]_
_q_ p q “ “
_j_ 0 _[b][j][z][j][, where]_
“

[ř][p]
ř _p_ _q_ _j_ !p!

_aj_ p ` ´ q _bj_ (5)
“ _p_ _q_ !j! _p_ _j_ ! _[,][ 0][ ď][ j][ ď][ p]_ “ [p][p]p[ `][ q] q[ ´]![ j]j![q][!]q[q][!][p´] j[1][q]![j] _[,][ 0][ ď][ j][ ď][ q.]_

p ` q p ´ q p ` q p ´ q

Note that, a0 _b0_ 1, and we now discuss the exponential of operators. Given an operator
(possibly non-linear), the “ “ rp{qs Pad´e approximation for e[L] can be written as _L_

_q_ ´1 _p_

_e[L]_ _rpq_ _bj_ _. . ._ _aj_ _. . ._ _._ (6)
« pLq “ ¨jÿ“0 _L ˝ Lj´ ˝times_ ˝ L˛ [¨]jÿ“0 _L ˝ Lj´ ˝times_ ˝ L˛

˝ ‚ ˝ ‚

For L being a linear transformation, the inverse in eq. (6) is computed as a matrix inverse, for examplelooooooomooooooon looooooomooooooon
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 Bpq is evaluated
pLq
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 Apq . Note that, both polynomial recurrent networks use same operator
pLq
_L, or in other words, the parameters are shared across the network. With aj, bj fixed-beforehand, the_
total trainable parameters of the Pad´e exponential model are _θ_ _, W, b_ .
t _L_ u

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** _θ_ _, a non-linearity layer v_ _σ_ _Wu_ _b_ _, and p, q_ N,
_L “ Lp_ _Lq_ “ p ` q P
_at points of differentiability, the gradients of the operation x_ _y_ _F_ _x; θ_ _, W, b_ : _p_ _q_ _e[L]_ _x_
ÞÑ “ p _L_ q “ r { s p q
_using the rp{qs Pad´e network in Figure 1 are bounded in operator norm by_
_y_ _nθ_ 2[¸][1][{][2]
B 2 2[q] BL _,_ (7)

BθL [ď][ exp][p}][L][}q p}][b][}] [` }][W] [}}][x][}] ˜jÿ“1 Bθj

_y_
B 2[,] (8)

_W_

[ď][ exp][p}][L][}q}][x][}]

B

_y_ _p_
B _._ (9)
_b_ _p_ _q_
B [ď][ exp] ˆ ` [}][L][}]˙

The detailed proof is provided in the Supplementary Materials (Appendix E). Note that,such that nθ _M_, where M is the dimension of the input. For example, the convolution operator has nθ “ |θL| is
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{qse[L][A] _, rp{qse[L][A]_ _, rp{qse[L][A]. The complete flow-diagram by_
plugging-in the Pad´e model is shown in Appendix D.

3 EXPERIMENTS

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 u0 _x_ and the output function u _τ, x_ for some
p q p q
pre-specified τ ą 0, are evaluated at M discretized locations of the domain Ω. This yields a single
training sample of _u0_ _xi_ _, u_ _τ, xi_ _, xi_ Ω, 1 _i_ _M_ . In total, we take N training samples, and
unless stated otherwise, p p Nq “ 1000p and we test onq P ď 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


-----

shared parameters

+ + +

+ + +


Figure 1: Pade Exponential Model´ . A recurrent neural architecture for computing the exponential of an operator
_L using rp{qs Pade approximation. The multiplicative scalar coefficients´_ _ai, bi are fixed-beforehand using eq._
(5). The non-linear fully-connected layer is used to mimic the inverse polynomial operation.

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.

3.1 KORTEWEG-DE VRIES EQUATION

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_ [B][u]x _x[3][, x][ P p][0][,][ 1][q][, t][ P p][0][,][ 1][s]_ (10)

B [´ B]B[3][u]

_u0_ _x_ _u_ _x, t_ 0 _._
p q “ p “ q

For Table 1 notations, we have L “ ´ B[B]x[3][3][,][ f] [p][u][q “ ´][0][.][25][u][2][, and][ N][ “] BBx [. The neural operator]

learns the mapping of the initial condition u0 _x_ to the solutions u _x, t_ 1 . The initial condition
p q p “ q
is generated in Gaussian random fields according to u0 0, 7[4] ∆ 7[2]I with periodic
boundary conditions. The data set is obtained by solving the equation using the fourth-order stiff „ N p p´ ` q[´][2][.][5]q
time-stepping scheme known as ETDRK4 (Cox & Matthews, 2002) with a resolution of 2[10], 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


-----

10

6 × 10


10

6 × 10


4 × 10

3 × 10


4 × 10


200 400 600 800 1000


200 400 600 800 1000


Padé Exp
MWT Leg
FNO


Padé Exp
MWT Leg
FNO


Figure 2: (Left) Number of training samples N vs performance (relative L2 error) for neural operators evaluated
on the KdV equation with s=1024. For N ă 1000, each smaller dataset is sampled uniformly randomly 5
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 s = 64 s = 128 s = 256 s = 512 s = 1024

Pad´e 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.12521


Table 2: Korteweg-de Vries (KdV) equation benchmarks for different input resolution s. The relative L2 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 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.


3.2 KURAMOTO-SIVASHINSKY (KS) EQUATION

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_ [B][u]x _x[2][ ´ B][4]x[u][4][, x][ P p][0][,][ 1][q][, t][ P p][0][,][ 1][s]_ (11)

B [´ B]B[2][u] B

_u0_ _x_ _u_ _x, t_ 0 _._
p q “ p “ q

The KS equation is also time-advection and according to the Table 1 notations, we have: L “
´ B[B]x[2][2][ ´] BBx[4][4][,][ f] [p][u][q “ ´][0][.][5][u][2][, and][ N][ “] BBx [. The neural operator learns the mapping of the initial]

condition u0 _x_ to the solutions u _x, t_ 1 . Similarly to the KdV equation in Section 3.1, the initial
p q p “ q
condition is sampled from a Gaussian random field u0 0, 5[4] ∆ 5[2]I with periodic
boundary conditions. The equation is numerically solved using „ N p _chebfunp´_ ` package (Driscoll et al.,q[´][2][.][5]q
2014) with a resolution of 2[10], 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.


-----

|MAE Networks Relative L2 error Net. vs FC C R D|Col2|
|---|---|
|Pade´ Exp 1219 ˘ 130 1752 ˘ 666 211 ˘ 31 0.0155 ˘ 0.0034 MWT Leg 3554 ˘ 1157 2928 ˘ 1338 284 ˘ 209 0.0245 ˘ 0.0043 FNO 3D 4213 ˘ 391 3391 ˘ 1233 592 ˘ 157 0.0301 ˘ 0.0045 LNO 3D 28502 ˘ 12698 6586 ˘ 3442 1465 ˘ 965 0.1056 ˘ 0.0394 Neural ODE 4339 ˘ 1174 3443 ˘ 1408 443 ˘ 192 0.0310 ˘ 0.0069 Seq2Seq 2798 ˘ 456 3317 ˘ 1690 346 ˘ 83 0.0273 ˘ 0.0058 Transformer 7087 ˘ 972 6613 ˘ 2853 1722 ˘ 320 0.0501 ˘ 0.0094 FC 10305 ˘ 2818 5885 ˘ 1609 1634 ˘ 686 0.0609 ˘ 0.0111|82.14% (+652K) 62.0% (+18M) 54.0% (+1.02M) -105.0% (+238K) 53.8% (+172K) 63.7%(+1.8M) 13.4% (+15.2K) (37.2K)|


Table 3: COVID-19 prediction benchmarks for different networks using 10-fold resampling with mean ˘ 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 each
model. The last column compares each network vs FC in terms of the total MAE improvement and total model
parameters difference.


|Confirmed|Col2|
|---|---|
|||
|actual Pade Exp seq2seq||
|||


6 3

time (days)


|Recovered|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||


6 3

time (days)

Recovered


Deaths


350000

325000

300000

275000

250000

225000

200000

112000

111000

110000

109000

108000


7500

7000

6500

6000

8500

8400

8300

8200

8100

8000

7900


800

600

400

200

0

100000

80000

60000

40000

20000


14

|Confirmed|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||



14


14

14


|actual Pade Exp seq2seq|Col2|
|---|---|


14

14


time (days)

Deaths


time (days)


time (days)


time (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 Pade Exp prediction and the best´
non-neural operator scheme from Table 3 (seq2seq) is shown. Top row: Most populous US state California with
population 39.77 M (2018 census). Bottom row: Same results for Massachusetts with a moderate population of
6.89 M.

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
(see results in Table 10, Appendix B.4). We again conclude that the proposed model attains better
performance than state-of-the-art 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. 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 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
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´1q “ pd0, d1, . . ., d6q.

_u0_ _x_ _u_ _τ,x_
p q p q

looooooooooomooooooooooon loooooomoooooon

**Forecasting The COVID19 forecasting benchmarks are presented in Table 3. 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.5) along with the best/worst prediction states.

4 ABLATION EXPERIMENTS

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 Pad´e neural model.

5 FUTURE DIRECTIONS AND CONCLUSION

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 3) 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](https://github.com/CSSEGISandData/COVID-19)


-----

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-1453860, the NSF award under Grant CCF-1837131, MCB-1936775,
CNS-1932620, the U.S. Army Research Office (ARO) under Grant No. W911NF-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, a
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-1548562. The third author has been supported in part by a NSF award under grant
DMS-2108900 and by the Simons Foundation. The views, opinions, and/or findings contained in
this article are those of the authors and should not be interpreted as representing the official views or
policies, either expressed or implied by the Defense Advanced Research Projects Agency, the Army
Research Office, the Department of Defense or the National Science Foundation.


-----

REFERENCES

M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions: With Formulas, Graphs,
_and Mathematical Tables._ Applied mathematics series. Dover Publications, 1965. ISBN
9780486612720.

Awad Al-Mohy and Nicholas Higham. A new scaling and squaring algorithm for the matrix exponential. SIAM Journal on Matrix Analysis and Applications, 31, 01 2009. doi: 10.1137/09074721X.

B. Alpert, G. Beylkin, D. Gines, and L. Vozovoi. Adaptive solution of partial differential equations in
multiwavelet bases. Journal of Computational Physics, 182(1):149–190, 2002. ISSN 0021-9991.

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/_
0912009.

Kevin Amaratunga and John Williams. Wavelet based green’s function approach to 2d pdes. Engi_neering Computations, 10, 07 2001. doi: 10.1108/eb023913._

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, 29(6):1716–1740, 1992. doi: 10.1137/0729097._

G. Beylkin, R. Coifman, and V. Rokhlin. Fast wavelet transforms and numerical algorithms i.
_Communications on Pure and Applied Mathematics, 44(2):141–183, 1991. doi: https://doi.org/10._
1002/cpa.3160440202.

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–259, 1997.
ISSN 0021-9991. doi: https://doi.org/10.1006/jcph.1996.5562.

Gregory Beylkin, James M. Keiser, and Lev Vozovoi. A new class of time discretization schemes for
the solution of nonlinear pdes. Journal of Computational Physics, 147(2):362–387, 1998. ISSN
0021-9991. doi: https://doi.org/10.1006/jcph.1998.6093.

Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik. Prediction of aerodynamic flow fields using convolutional neural networks. Computational Mechanics,
64(2):525–545, Jun 2019. ISSN 1432-0924. doi: 10.1007/s00466-019-01740-0.

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. A fractional calculus approach to modeling fractal dynamic
games. In 2011 50th IEEE Conference on 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, R. Hadsell, M. F. 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
differential equations, 2019.

Yi-Cheng Chen, Ping-En Lu, Cheng-Shang Chang, and Tzu-Hsuan Liu. A time-dependent sir model
for covid-19 with undetectable infected persons. IEEE Transactions on Network Science and
_Engineering, 7(4):3279–3294, Oct 2020. ISSN 2334-329X._


-----

Steven M Cox and Paul C Matthews. Exponential time differencing for stiff systems. Journal of
_Computational Physics, 176(2):430–455, 2002._

Jessica T. Davis, Matteo Chinazzi, Nicola Perra, Kunpeng Mu, Ana Pastore y Piontti, Marco Ajelli,
Natalie E. Dean, Corrado Gioannini, Maria Litvinova, Stefano Merler, Luca Rossi, Kaiyuan Sun,
Xinyue Xiong, M. Elizabeth Halloran, Ira M. Longini, Cecile Viboud, and Alessandro Vespignani.´
Estimating the establishment of local transmission and the cryptic phase of the covid-19 pandemic
in the usa. medRxiv, 2020.

Odo Diekmann, Hans G. Othmer, Robert Planque, and Martin C. J. Bootsma. The discrete-time´
kermack–mckendrick model: A versatile and computationally attractive framework for modeling
epidemics. Proceedings of the National Academy of Sciences, 118(39), 2021. ISSN 0027-8424.
doi: 10.1073/pnas.2106332118.

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.

T. A Driscoll, N. Hale, and L. N. Trefethen. Chebfun 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/pnas.2025632118._

Wyman Fair and Yudell L. Luke. Pade approximations to the operator exponential.´ _Numer. Math., 14_
(4):379–382, March 1970. ISSN 0029-599X. doi: 10.1007/BF02165592.

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:1233–1256, 01
2019. doi: 10.1137/18M1228876.

Jordi Feliu-Faba, Yuwei Fan, and Lexing Ying. Meta-learning pseudo-differential operators with deep`
neural networks. Journal of Computational Physics, 408:109309, May 2020. ISSN 0021-9991.
doi: 10.1016/j.jcp.2020.109309.

Silviu-Ioan Filip, Aurya Javeed, and Lloyd Trefethen. Smooth random functions, random odes, and
gaussian processes. SIAM Review, 61:185–205, 01 2019. doi: 10.1137/17M1161853.

Xiaoxiao Guo, Wei Li, and Francesco Iorio. Convolutional neural networks for steady flow approximation. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge
_Discovery and Data Mining, KDD ’16, pp. 481–490. Association for Computing Machinery, 2016._

Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavelet-based 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 H. 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:1214–1233, 01 2005. doi: 10.1137/S1064827502410633._

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–accelerated computational fluid dynamics. Proceedings of the National Academy
_of Sciences, 118(21), 2021. ISSN 0027-8424. doi: 10.1073/pnas.2101784118._


-----

S.G. Krantz. A Handbook of Real Variables: With Applications to Differential Equations and Fourier
_Analysis. A Handbook of Real Variables: With Applications to Differential Equations and Fourier_
Analysis. Birkh¨auser, 2004.

E. N. Laguerre. Le calcul des syst´emes lin´eaires, extrait d’une lettre adress´ee a. m. hermite. Extrait
_du Journal de l’Ecole Polytechnique, LXII. Cahier[´]_ [, pp. 221–267, 1898. URL http://gallica.](http://gallica.bnf.fr/ark:/12148/bpt6k90210p/ f242)
[bnf.fr/ark:/12148/bpt6k90210p/f242.](http://gallica.bnf.fr/ark:/12148/bpt6k90210p/ f242)

J. Douglas Lawson. Generalized runge-kutta processes for stable systems with large lipschitz
constants. SIAM Journal on Numerical Analysis, 4(3):372–380, 1967. ISSN 00361429.

Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew
Stuart, and Anima Anandkumar. Fourier neural operator for parametric 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. In 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–452, 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 of operators, 2020.

Alejandro Molina, Patrick Schramowski, and Kristian Kersting. Pade activation units: End-to-end´
learning of flexible activation functions in deep networks. In International Conference on Learning
_[Representations, 2020. URL https://openreview.net/forum?id=BJlBSkHtDS.](https://openreview.net/forum?id=BJlBSkHtDS)_

Ravi G. Patel, Nathaniel A. Trask, Mitchell A. Wood, and Eric C. Cyr. A physics-informed operator
regression framework for extracting data-driven continuum models. Computer Methods in Applied
_Mechanics and Engineering, 373:113500, 2021. ISSN 0045-7825. doi: https://doi.org/10.1016/j._
cma.2020.113500.

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.1101/2020.03.21.20040303.

M. Raissi, P. 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-9991.

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, 117(17):9183–9193, 2020. ISSN 0027-8424._
doi: 10.1073/pnas.1922204117.

David Salinas, Valentin Flunkert, and Jan Gasthaus. Deepar: Probabilistic forecasting with autoregressive recurrent networks, 2019.


-----

Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz
Kaiser, and Illia Polosukhin. Attention is 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.1101/2020.04.16.20067306.

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-informed deeponets, 2021b.

J. Wendt. Computational Fluid Dynamics: An Introduction. A von Karman Institute book. Springer
Berlin Heidelberg, 2008. ISBN 9783540850557.

Yuankun Xue and Paul Bogdan. Constructing compact causal mathematical models for complex
dynamics. 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 9780387102108.

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 Info.Theory, 66(3):1738–1759, 2020a. doi: 10.1109/TIT.2019.2961812._

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. medRxiv, 2020b. doi:
10.1101/2020.05.24.20111989.


-----

Networks _λ “ 0.05_ _λ “ 0.04_ _λ “ 0.03_ _λ “ 0.02_

Pad´e Exp (p=5, q=6) 0.00199 0.00237 0.00219 0.00215
Pad´e Exp (p=4, q=2) 0.00222 0.00231 0.00291 0.0159

Taylor (Lt=6) 0.00730 0.00690 0.00846 0.00704
Taylor (Lt=4) 0.00938 0.00842 0.00866 0.0193

L’Hospital (n=6) 0.00775 0.00699 0.00752 0.0194
L’Hospital (n=4) 0.00731 0.00798 0.00844 0.0199


Table 4: Taylor vs Pade Approximation for exponential operator evaluated on the KdV equation for s=1024 with´
varying input (u0pxq) fluctuation strength λ (low λ implies high fluctuation). The highest polynomials degree
are similar for each model by setting Lt “ maxpp, qq for Taylor and similarly n “ maxpp, qq for L’Hospital.

1

0.0036

0.0032

p 5

0.0030

7

0.0028

9

0.0026


Figure 4: Polynomials degree p, q grid variation for rp{qs Pade neural model of Figure 1. The performance is´
evaluated on the KdV equation with s=1024.

A RELATED WORK


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

B ADDITIONAL RESULTS


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


_Lt_


_e[L]_ «


1

_e[L]_ _,_
« _j!_

_j_ 0

ÿ“ [p][L][ ˝][ L]j´[ ˝]times[ . . .][ ˝][ L][q]

where, Lt is the truncated length of the series. We evaluate the truncated Taylor series approximationloooooooomoooooooon
vs rp{qs Pade approximation in Figure 1 for same number of maximum polynomial degree, i.e.,´
_Lt_ 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 3.1). The fluctuated “ p q


-----

input is sampled from Random fields Filip et al. (2019) with controllable parameter λ with lower
values attaining higher fluctuations and the input u0 _x_ reaches the Brownian motion limit with
p q
_λ Ñ 0. We see in Table 4 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 exp** _x_ can also be approximated as limn 1 _x_ _n_
p q Ñ8p ` { q[n]
using L’Hospital rule (Krantz, 2004). For a given operator L, the L’Hospital based exponential
operator approximation can be written as

_e[L]_ « pI ` L{nq ˝ pI ` L{nq ˝ . . . ˝ pI ` L{nq,

_n´times_

where, I is the identity operator andlooooooooooooooooooooooooomooooooooooooooooooooooooon n is the length of truncation of operator compositions. We
compare L’Hospital approximation with the Pade model in Table 4 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** Depth L2 error
layer σpWu ` bq in the Pade neural model for all the experiments´ 1 0.00295
with σp.q=ReLU. Now, we vary the depth of this layer from 1 to 3 2 0.00339
(each layer is followed by ReLU non-linearity) to see the individual 3 0.00628
contribution of this layer towards the final test performance. The
experiments settings are same as for Table 2 and we fixUpon increasing the depth we do not see any improvements, in s “ 1024. Table 5: Non-linear layer depth vstest error for the Pade neural model.´
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


Test
s = 2048 s = 4096 s = 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: Pade exponential (MWT Leg) model trained at lower resolutions can predict the output 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 _θ_ _, W, b_, the model can be trained and tested on different values of _p, q_ by just
t _L_ u p q
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.


-----

1 0.8 1 0.8

3 3

0.6 0.6

p5 0.4 p5 0.4

7 7

0.2 0.2

9 9

1 3 5 7 9 1 3 5 7 9

q q


Figure 5: Relative L2 error when Pade model trained and tested on different´ pp, qq. (Left) The KdV data using
pp, qq “ p5, 6q for training and other pp, qq for testing. (Right) The KdV data using pp, qq “ p4, 2q for training
and other pp, qq 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 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 KURAMOTO-SIVASHINSKY EQUATION

Comparison of relative error for different neural operators evaluated over KS equation is presented in
Table 10.

B.5 COVID-19 FORECASTING

Additional forecasting results for COVID-19 by considering the corner cases are presented in Figure 8.

B.6 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:

_wt_ _x, t_ _u_ _x, t_ _w_ _x, t_ _ν∆w_ _x, t_ _f_ _x_ _,_ _x_ 0, 1 _, t_ 0, T
p q ` p q ¨ ∇ p q ´ p q “ p q P p q[2] P p s

_∇_ ¨ upx, tq “ 0, _x P p0, 1q[2], t P r0, T_ s (12)

_w0_ _x_ _w_ _x, t_ 0 _,_ _x_ 0, 1
p q “ p “ q P p q[2]

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 vorticity w. The initial condition is generated in Gaussian
´ 3
random fields according to w0 0, 7 2 ∆ 7[2]I with periodic boundary conditions and
the forcing function is f _x_ „ 0 N.1 psin 2π p´x1 ` x2 q[´][2] cos[.][5]q 2π _x1_ _x2_ . We experiment with
different viscosities ν, final timep q “ _Tp, and the number of training pairsp_ p ` qq ` p p ` N : pqqqiq ν “ 1e ´ 3, T “ 50,


-----

_ν “ 1e ´ 3_ _ν “ 1e ´ 4_ _ν “ 1e´4_ _ν “ 1e ´ 5_
Networks _T “ 50_ _T “ 30_ _T “ 30_ _T “ 20_
_N “ 1000_ _N “ 1000_ _N “ 10000_ _N “ 1000_

Pad´e 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 at various viscosities ν and prediction horizon T .

Networks N = 200 N = 400 N = 600 N =800 N = 1000

Pad´e Exp **0.00864˘5.1e-4** **0.00439˘ 2.8e-4** **0.00365˘ 2.2e-4** **0.00322˘ 1.7e-4** **0.00295**
MWT Leg 0.00898˘16.1e-4 0.00641˘7.7e-4 0.00463˘3.5e-4 0.00420˘2.7e-4 0.00392
FNO 0.00970˘6.4e-4 0.00781˘3.3e-4 0.00706˘2.1e-4 0.00679˘1.2e-4 0.00672


Table 8: Korteweg-de Vries (KdV) equation benchmarks for different numbers of training samples N . Top: Our
method. Bottom: the state-of-the-art methods.

Networks N = 200 N = 400 N = 600 N =800 N = 1000

Pad´e Exp **0.00764˘2.0e-4** **0.00489˘ 2.2e-4** **0.00416˘ 1.8e-4** **0.00376˘ 1.0e-4** **0.00338**
MWT Leg 0.00849˘5.7e-4 0.00612˘5.5e-4 0.00496˘4.2e-4 0.00478˘3.3e-4 0.00445
FNO 0.01024˘1.1e-3 0.00771˘3.4e-4 0.00625˘1.7e-4 0.00508˘1.4e-4 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; 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 “ 4 × 10 1
3 are concatenated using ReLU nonlinearity. The results are reported
in Table 7, and we observe that Pade´ 3 × 10

recent state-of-the-art neural operator
approaches. For the less data setup of 2 × 10 1
_Nin Figure 6 that the Pad “ 1000 with ν “ 1ee Exp model´ ´ 4, we see_ Relative error
quickly converges to the lowest value


B.7 NUMERICAL VALUES

4 × 10 1

1 Padé Exp
3 × 10

MWT Leg
FNO

2 × 10 1

Relative error

1 100 200 300 400 500

epochs

Figure 6: Relative L2 error vs epochs for MWT Leg with different

The numerical values for the training number of OP basis k.
samples variation experiment for the
KdV and KS equation in the Figure 2
are shown in Table 8 and Table 9, respectively.


-----

Networks s = 64 s = 128 s = 256 s = 512 s = 1024

Pad´e 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 10: Kuramoto-Sivashinsky (KS) equation benchmarks for different input resolution s. The relative L2
errors are shown shown for each model.

C NOTATIONS

**Operator Learning**

_T, L_ Operators between function spaces

_H[s,p]_ Sobolev spaces such that constituent functions and their weak
derivatives upto order s have finite L[p] norms
**Multiwavelets**

Subspace addition

**Vàn[k]** _f_ _f are polynomials of degree_ _k defined over interval_
t | ă
p2[´][n]l, 2[´][n]pl ` 1qq for all l “ 0, 1, . . ., 2[n] ´ 1, and assumes
0 elsewhereu

**Wn[k]** Orthogonal space to Vn[k] [such that][ W]n[k] **Vnk** _n_ 1

[“][ V][k]`
_Pn_ Projection operator such that Pn : À **Vn[k]**
_H[s,][2]_ Ñ
_Qn_ Projection operator such that Qn : **Wn[k]**
_H[s,][2]_ Ñ
_L_ Coarsest scale of the multiwavelet transform
**Pad´e Approximation**

}.} Operator norm

_._ 2 Euclidean L-2 norm
} }

rp{qsf Pade approximation of function´ _f with numerator/denominator_
polynomial degree p/q, respectively


_θ_ Parameters used to represent the operator
_L_ _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 Θ

D MULTIWAVELET EXPONENTIAL OPERATOR ARCHITECTURE

The Pade neural operator based multiwavelet transform model is shown in Figure 7. The input to the´
model is s[p][n][`][1][q], where n log _M_, and the output is Us, where n is the finest scale (or the log of
“ p q [p][n][q]
input resolution). For a detailed description of multiwavelet transform, we refer the reader to (Gupta
et al., 2021).

E PROOF OF THEOREM 1

First note that, for the Pad´e approximation rp{qse[x], the polynomial coefficients in eq. (5), we have

_j_

_p_ _p_ 1 _p_ _j_ 1 _p_

_aj_ [1] p ´ q ¨ ¨ ¨ p ´ ` q _._ (13)
“ _j!_ _p_ _q_ _p_ _q_ 1 _p_ _q_ _j_ 1 _j!_ _p_ _q_

p ` qp ` ´ q ¨ ¨ ¨ p ` ´ ` q [ď][ 1] ˆ ` ˙


-----

## +


## +


Figure 7: Multiwavelet Pade Exponential Model´ . The Pade exponential neural operator from Figure 1 is used´
in the skeleton of the multiwavelet transform based neural operator model. The inputs and outputs are recursively
updated using the decomposition cell (Left) and reconstruction cell (Right).


and similarly,

_j_

_q_

_bj_ [1] _._
| | ď _j!_ _p_ _q_

ˆ ` ˙

For the simplicity of notation: let A _z_ _Apq_ _z_ _j_ 0 _[a][j][z][j][,][ A][1][p][z][q “][ ř]j[p]_ 1 _[ja][j][z][j][´][1][,][ C][p][z][q “]_
_q_ p q “ p q “ “ “
_j_ 0 _j_ 1 _[j][|][b][|][j][z][j][´][1][. Note][ |][B][pq][p][z][q| ď][ C][p|][z][|q][ and][ |][B]pq[1]_ [p][z][q| ď][ C] [1][p|][z][|q][. Using eq.]
(13), we can further write that“ [|][b][j][|][z][j][,][ C] [1][p][z][q “][ ř][q]“ [ř][p]
ř

_p_ _p_ _p_
_A_ _z_ exp _,_ _A[1]_ _z_ (14)
| p q| ď _p_ _q_ | p q| ď _p_ _q_ [exp] _p_ _q_
ˆ ` [|][z][|]˙ ` ˆ ` [|][z][|]˙

_q_ _q_ _q_
_C_ _z_ exp _,_ _C_ [1] _z_ (15)
| p q| ď _p_ _q_ | p q| ď _p_ _q_ [exp] _p_ _q_
ˆ ` [|][z][|]˙ ` ˆ ` [|][z][|]˙

Further let d and nθ denote the dimensions of x and θ respectively: x, y, u, v R[d] and θ : _θ_ R[n][θ] .
For the operation x _y, we denote y_ _F_ _x; θ_ _, W, b_ : _p_ _q_ _e[L] Px_ . Then y“ _L A P_ _v,_
ÞÑ “ p _L_ q “ r { s p q “ pLq
_v “ σpWu ` bq, u “ BpLqx, and_
_y_ 1
B max

BθL [“] _gPR[nθ]_ _,}g}2“1_ _tÑ0_ _t_ [p][F] [p][x][;][ θ][ `][ tg, W, b][q ´][ F] [p][x][;][ θ, W, b][qq] 2

max _θ_ _y, h_ _, g_ (16)
“ _g_ R[n][θ] _,_ _g_ 2[lim] 1 |x∇ x y y|
_h P_ R[d], }h }2 “ 1
P } } “

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,θ “ _j“1_ _[g][j]_ BBθj [, the]

differential operator induced by g. The quantity of interest is _g,θ_ _y, h_ . Next, by chain rule
|B x y|
_g,θ_ _y, h_ _g,θA_ _v, h_ _A_ _g,θv, h_ [ř][n][θ]
B x y “ xpB pLqq y ` x pLqB y
_g,θA_ _v, h_ _A_ _DW_ _g,θB_ _x, h_ (17)
“ xpB pLqq y ` x pLq B pLq y
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_ _B_ _x_ 2 _b_ 2.
Substituting back in eq. (17) we obtain } } ď } }} } `} } ď } }} pLq}} } `} }
_g,θ_ _y, h_ _g,θA_ _W_ _B_ _x_ 2 _g,θA_ _b_ 2
|B x y| ď }B pLq}} }} pLq}} } ` }B pLq}} }
_A_ _W_ _g,θB_ _x_ 2 (18)
` } pLq}} }}B pLq}} }
Using eq. (14)-(15), the spectral norms are bounded further by


}ApLq} ď Ap}L}q ď exp

}BpLq} ď Cp}L}q ď exp


_p_

_,_ (19)
_p_ _q_
` [}][L][}]˙

_q_

_._ (20)
_p_ _q_
` [}][L][}]˙


-----

|Confirmed|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||


6 3


|Recovered|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||


6 3


|Deaths|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||
|||


6 3

time (days)

Deaths


900

850

800

750

700

650

8000

7500

7000

6500

6000


1150

1100

1050

1000


25

20

15

10

1600

1500

1400


950

900

55000


14

|time (days) Confirmed|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||



14


14

|time (days) Recovered|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||



14


14

|actual Pade Exp seq2seq|Col2|
|---|---|



14


50000

45000


time (days)

Confirmed

|actual Pade Exp seq2seq|Col2|
|---|---|



6 3

time (days)

Confirmed


time (days)

Recovered

|actual Pade Exp seq2seq|Col2|
|---|---|



6 3

time (days)


time (days)

Deaths

|actual Pade Exp seq2seq|Col2|
|---|---|



6 3

time (days)


2500

2000


65000

60000

55000

50000

45000

40000

0

2000


100000

90000

80000

70000

60000

64000

62000

60000

58000


1500

1000


14

|actual Pade Exp seq2seq|Col2|
|---|---|



14


14

|Recovered|Col2|
|---|---|
|||
|||
|actual Pade Exp seq2seq||



14


14

|actual Pade Exp seq2seq|Col2|
|---|---|



14


Deaths

6 3


1500

1250

1000

750

500

250


4000


time (days)

Confirmed



150000

100000

50000

0

2.5

2.0

1.5

1.0

0.5

0.0


210000

207500

205000

202500

200000

2.75

2.70

2.65

Counts (Texas)

2.60


4000

3000

2000

1000

50000

|actual Pade Exp seq2seq|Col2|
|---|---|

|time (days) Recovered|Col2|
|---|---|
|||
|actual Pade Exp seq2seq||
|||

|time (days) Deaths|Col2|
|---|---|
|||
|actual Pade Exp seq2seq||
|||


14 9 6 3

time (days)


14 9 6 3

time (days)


14

|time (days) Deaths|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||



14


time (days)

Deaths

6 3

time (days)


50000


|time (days) 6 Confirmed|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||


14


|time (days) 6 Recovered|Col2|
|---|---|
|actual Pade Exp seq2seq||
|||


14


time (days)


time (days)


Figure 8: COVID-19 Forecasting. The Pade Exp prediction and the best non-neural operator scheme from´
Table 3 (seq2seq) is shown. (Top 3 rows) Prediction for the week 06/12/20 – 06/18/20 using previous 2 weeks
data. This week is the test sample with least MAE averaged across Confirmed, Recovered, and Death counts
for Pade Exp´ ; or the test sample with best averaged prediction. Next, row 1-3 are US states in the order
of best-middle-worst averaged prediction for this test sample. (Bottom 3 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 are US states in the order of best-middle-worst averaged prediction for this worst test sample. The
missing values are substituted with zeroes which is one of the reason for higher errors.


-----

Next, by using the product rule, we obtain

_p_

_g,θA_ _aj_ _g,θ_ _. . ._
}B pLq} ď | |}B pL ˝ L ˝ ˝ Lq}

_jÿ“1_ _j´times_

_p_

looooooomooooooon

_a_ _jaj_ _g,θ_
p q }L}[j][´][1]}B _L}_
ď _jÿ“1_

_A[1]_ _g,θ_
“ p}L}q}B _L}_


_p_

_g,θ_ _._ (21)
_p_ _q_ }B _L}_
` [}][L][}]˙


_p_ _q_ [exp]
`


Similarly,


_q_ _q_
_g,θB_ _C_ [1] _g,θ_ _g,θ_ _,_ (22)
}B pLq} ď p}L}q}B _L} ď_ _p_ _q_ [exp] _p_ _q_ }B _L}_

` ˆ ` [}][L][}]˙

_nθ_ _nθ_ 2[¸][1][{][2]

}Bg,θL} ď _j_ 1 |gj|} [B]θ[L]j } ď }g}2 ˜j 1 } [B]θ[L]j } _._ (23)

ÿ“ B ÿ“ B

Using eq. (18), (22), and (23), we obtain the bound for eq. (16) as

_nθ_ 2[¸][1][{][2]

_y_ _p_ _p_

BBθL [ď] ˆexpp}L}q}W }}x}2 ` _p ` q_ [exp] ˆ _p ` q_ [}][L][}]˙ }b}2˙ [˜]jÿ“1 } B[B]θ[L]j } _,_ (24)

which is used to write the eq. (7) in Theorem 1. Similarly, the analysis is extended as

_y_ 1
B max

BW [“] ΦPR[d][ˆ][d],}Φ}2“1 _tÑ0_ _t_ [p][F] [p][x][;][ θ, W][ `][ t][Φ][, b][q ´][ F] [p][x][;][ θ, W, b][qq] 2

max (25)
“ ΦPR[d][ˆ][d],}Φ}2“1 [}][A][lim][p][L][q][D][Φ][Φ][B][p][L][q][x][}][2][ ď }][A][p][L][q}}][B][p][L][q}}][x][}][2][,]

and

_y_ 1
B max
Bb [“] _hPR[d],}h}2“1_ _tÑ0_ _t_ [p][F] [p][x][;][ θ, W, b][ `][ th][q ´][ F] [p][x][;][ θ, W, b][qq] 2

max (26)
“ _hPR[d],}h}2“1_ [}][A][lim][p][L][q][D][h][h][}][2][ ď }][A][p][L][q}][,]

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.


-----

