Multiwavelet-based Operator Learning for Differential Equations
Authors: Gaurav Gupta
Xiongye Xiao
Paul Bogdan
Venue: NeurIPS 2021 Spotlight (2021)
Keywords: Neural operators
multiwavelet transform
partial differential equations

Abstract: The solution of a partial differential equation can be obtained by computing the inverse operator map between the input and the solution space. Towards this end, we introduce a $\textit{multiwavelet-based neural operator learning scheme}$ that compresses the associated operator's kernel using fine-grained wavelets. By explicitly embedding the inverse multiwavelet filters, we learn the projection of the kernel onto fixed multiwavelet polynomial bases. The projected kernel is trained at multiple scales derived from using repeated computation of multiwavelet transform. This allows learning the complex dependencies at various scales and results in a resolution-independent scheme. Compare to the prior works, we exploit the fundamental properties of the operator's kernel which enable numerically efficient representation. We perform experiments on the Korteweg-de Vries (KdV) equation, Burgers' equation, Darcy Flow, and Navier-Stokes equation. Compared with the existing neural operator approaches, our model shows significantly higher accuracy and achieves state-of-the-art in a range of datasets. For the time-varying equations, the proposed method exhibits a ($2X-10X$) improvement ($0.0018$ ($0.0033$) relative $L2$ error for Burgers' (KdV) equation). By learning the mappings between function spaces, the proposed method has the ability to find the solution of a high-resolution input after learning from lower-resolution data.

Multiwayelet-based Operator Learning for Differential Equations
Gaurav Gupta, Xiongye Xiao, Paul Bogdan
Ming Hsieh Department of Electrical and Computer Engineering University of Southern California, Los Angeles, CA 90089 {ggaurav, xiongyex, pbogdan)(usee edu
Abstract
The solution of partial differential equation can be obtained by computing the inverse operator map between the input and the solution space. Towards this end, introduce mulriwavvlet-baaed neural operator learning scheme that com- presses the associated operator's 's kernel using fine- grained wavelets. By explicitly embedding the inverse multiwavelet filters, learn the projection of the kernel onto fixed multiwavelet polynomial bases. The projected kernel is trained at mul- tiple scales derived from using repeated computation of multi transform. This allows learning the complex dependencies at various scales and results resolution-independent scheme. Compare the prior works, we exploit the fundamental properties of the operator's kernel which enable numerically efficient representation. We perform experiments on the Korteweg-de Vries (KdV) equation, Burgers' equation, Darcy Flow, and Kavier-Stoker equation. Compared with the existing neural operator approaches, our model shows significantly higher accuracy and achieves state-of-ie titte-of-the-art in range of datasets. For the time-varying equations, the proposed method exhibits a (2X 10X) improvement (0.0015 (0.0033) relative 7.2 error for Burgers' (KdV) equation). By learning the mappings between function spaces, the proposed method has the ability to find the solution of a high-resolution input after learning from lower-resolution data.
Introduction
Many natural and human-built systems (e.g., aerospace, complex fluids, neuro-glia information processing) exhibit complex dynamics characterized by partial differential equations (PDEs) [52] 60] For example, the design of wings and airplanes robust to turbulence, requires to learn complex PDEs. Along the same lines, complex fluids (gels, emulsions) are multiphasic materials characterized by macroscopic behavior (55] modeled by non-linear PDEs. Understanding their variations in viscosity as a function of the shear rate is critical for many engineering projects. Moreover, modeling the dynamics of continuous and discrete cyber and physical processes in complex cyber-physich systems can be achieved through PDES [68].
Recent efforts on learning PDEs (i.e., mappings between fifnite-dimensionae spaces of functions), from trajectories of variables, focused on developing machine learning and particular deep neural networks (NNs) techniques. Towards this end, stream of work aims at parameterizing the solution map as deep NNs 2 13 40 71). One issue, however, is that the NNs are tied specific resolution during training, and therefore, may not generalize well to other resolutions, thus, requiring retraining (and possible modifications of the model) for every set of discretizations. In parallel, another stream of work focuses on constructing the PDE solution function
NN architecture
31 42 57 65 This approach, however, is designed to work with one instance of a PDE and, therefore, upon changing the coefficients associated with the PDE, the model has to be re-trained.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Additionally, the approach is not complete data-depenne detendeeenndntt one, and hence, cannot be made oblivious the knowledge of the underlying PDE structure. Finally, the closest stream of work to the problem
we investigate is represented by the "Neural Operators" [4 47
48
49
56
Being complete
data-driven approach, the neural operators method aims at learning the operator map without having knowledge of the underlying PDEs. The neural operators have also demonstrated the capability of Obtaining the data for learning the operator map could be prohibitively expensive or time consuming (e.g., aircraft performance to different initial conditions). To be able to better solve the problem of learning the PDE operators from scarce and noisy data, we would ideally explore fundamental properties of the operators that have implications in data-efficient representation.
Our intuition is to transform the problem of learning PDE to a domain where a compact representa- tion of the operator exists. With mild assumption regarding the smoothness of the operator's kernel, except finitely many singularities, the multiwavelets (5), with their vanishing moments property, sparsify the kernel in their projection with respect to (w.r.t.) measure. Therefore, learning an operator kernel in the multiwavelet domain is feasible and data efficient. The wavelets have a rich history in signal processing [24 225, and are popular in audio, image compression [8 61.. For multiwavelets, the orthogonal polynomial (OP) w.r.t. a measure emerges as a natural basis for the multiwavelet subspace, and appropriate scale / shift provides sequence of subspaces which captures the locality at various resolutions. We generalize and exploit the multiwavelets concept to work with arbitrary measures which upens-up new possibilities to design series of models for the operator learning from complex data streams.
We incorporate the multiwavelet filters derived using variety of the OP basis into our operator learning model, and show that the proposed architecture outperforms the existing neural operators. Our main contributions are as follows: (i) Based on some fundamental properties of the integral operator's kernel, we develop a multiwavelet-basee model which learns the operator map efficiently. (ii) For the 1-D dataset of non-linear Korteweg-de Vries and Burgers equations, we observe an order of magnitude improvement in the relative L2 error (Section 3.1 3.3). (iii) We demonstrate that the proposed model is in validation with the theoretical properties of the pseudo-differential operator (Section3.2 (iv) We show how the proposed multiwavelet based model is robust towards the fluctuation strength of the input signal (Section3..) (v) Next, we demonstrate the applicability on higher dimensions of 2-D Darcy flow equation (Section 3.4). and finally show that the proposed approach can learn at lower resolutions and generalize higher resolutions. The code for reproducing the experiments is available at:
hwtp://aaaaa//1111//ttttooo operator
2 Operator Learning using Multiwavelet Transform
We start by defining the problem of operator learning in 2.1 Section 2.2 defines the multi- wavelet transform for the proposed operator learning problem and derives the necessary transformation operations across different scales. Section outlines the proposed operator learning model. Fi- nally, Section 2.4 lists some of the useful properties of the operators which leads to an efficient implementation of multivavelet-basee models.
2.1 Problem Setup
Given two functions a(x) and u(x) with x E D, the operator is a map T such that Ta = u. Formally, let A and L/ be two Sobolev spaces HPP (s 0,p 1), then the operator T is such that T A - U. The Sobolev spaces are particularly useful in the analysis of partial differential equations (PDEs), and we restrict our attention to > 0 and p 2. Note that, for s O, the Hd coincides with LP. and, f E HP does not necessarily have derivatives in LP, We choose = 2 in order to be able to define projections with respect to (w.rt.)) measures j in a Hilbert space structure.
We take the operator T as an integral operator with the kernel K D x D + L² such that
Ta(x) - 1/ K(x, K(I.y)a(ydd)
(1)
For the case of inhomogeneous linear PDEs, Lu = S. with ʃ being the forcing function, L the differential operator, and the associated kernel is commonly termed as Green function. In our case, we do not put the restriction of linearity on the operator: From eq. ② is apparent that learning the complete kernel K (....) would essentially solve the operator map problem, but it is not necessarily

---

numerically feasible solution. Indeed, better approach would be to exploit possible useful properties (see Section such that a compact representation of the kernel can be made. For an efficient representation of the operator kernel, we need an appropriate subspace (or sequence of subspaces), and projection tools to map to such spaces.
Norm with respect to measures: Projecting given function onto fixed basis would require measure dependent distance. For two functions f and g. we take the inner product w.r.t measure μ
as (J, ub = f (())g(()) (x)))(x)dd and the associated norm as I| = [f 1/2 We now discuss the next ingredient, which refers the subspaces required to project the kernel.
2.2 Multiwavelet Transform
In this section, we briefly overview the concept of multiwavelets 4 and extend it to work with non- uniform measures at each scale. The multiwavelet transform synergizes the advantages of orthogonal polynomials (OPs) as well as the wavelets concepts, both of which have a rich history in the signal processing. 'he properties of wavelet bases like (7) vanishing moments, and (ii) orthogonality can effectively be used to create a system of coordinates in which a wide class of operators (see Section 2.4) have nice representation. Multiwavelets go few steps further, and provide fine fine-grained representation using OPs, but also act as basis on a finite interval. For the rest of this section, we restrict our attention the interval [O, I]; however, the transformation any finite interval [a, bj could be straightforwardly obtained by an appropriate shift and scale.
Multi Resolution Analysis: We begin by defining the space of piecewise polynomial functions, for k E and E Z+ U {0} as, V " [flde(()) < for x (2-"1, 21-"1,--((1- + 1)) O, elsewhere]. Clearly, dim(Vá)
k and for subsequent n, each subspace is contained in
another as shown by the following relation:
V c v V n-
V
Similarly, we define the sequence of measures μo: Pl,... such that f V is measurable w.r.t 44x and the norm of is taken as ||ʃl| (J, J'?." Next, since V n-1 C Vk. we define the multiwavelet subspace W for Z+ U {0), such that
V
W n
V T W
(3)
For given OP basis for Va as do, 01.- Φh-1 w.r.t. measure Mo, basis of the subsequent spaces V& n > 1 can be obtained by shift and scale (hence the name, multi-scale) operations of the original basis as follows:
(x)()) =""))="22,,""'> -)), j= l.......kk 1, 1= 0,1,....,22 w.r.t.
Jny
where, is obtained as the collections of shift and scale of μo. accordingly.
Multiwavelets: For the multiwavelet subspace W the orthonormal basis (of piecewise polynomials) are taken as vo, ..,,,..,, .h-1 such that (vh, W5) = 0 for & j and 1, otherwise. From eq. (3).
V 1 W n: and since Vk spans the polynomials of degree at most k therefore, we conclude that
{ x v (x)d O, AA j k. (vanishing moments)
(5)
Similarly to eq. a basis for multiwavelet subspace W are obtained by shift and scale of s/: as (z) 2n (2"z 1) and 2 are orthonormal w.r.t. measure nng i.e. (W dpr),. if jfj=j. 1 1', and 0 otherwise. Therefore, for a given OP basis for V (for example, Legendre, Chebyshev polynomials), we only require compute v., and complete basis set all the scales can be obtained using scale/shift of d;, vi.
Note: Since V v + W from eq. B therefore, for given basis 6; of Va w.r.t. measure #o and 03 basis for Vt w.r.t. μl, set of basis fia can be obtained by applying Gram-Schmidt
3

---

K[ v)
(ii)
(ii)
H
H G
coarse
MAg
G
reconstuct
Bo
Ao
(0) H
decompose
GO
G
fine
B1
Figure Multiwavelet representation of the Kernel. (i) Given kemee ((x,, y) of an integral operator T, (ii) the bases with different mcasures (fo, (f1) at two different scales (coarse=0, finc=l) projects the kemne into components A;, B:. C. (iii) The decomposition yields sparse structure, and the entries with absolute magnitude values exceeding le shown in black. Given projections any scale, the finer coarser scale projections can be obtained by reconstruction / decomposition using a fixed multiwavelct filters H and G 1.
Orthogonalization using appropriate measures. We refer the reader to supplementary materials for the detailed procedure.
Note: Since Va and W lives in Vt. therefore, 6;, v can be written as linear combina- tion of the basis of V/A. We term these linear coefficients as multiwavelet decomposition filters (H(0), H ( (())), since they are transforming a fine n = 1 to coarse scale n 0. A uniform measure (1t) version discussed in 4. and we extend any arbitrary measure by including the correction terms ()) and E(1) We refer supplementary materials for the complete details. The capability of using the non-uniform measures enables us to apply the same approach to any OP basis with finite domain, for example, Chebyshev, Gegenbauer,
For a given f(x), the multiscale, multiwavelet coefficients at the scale n are defined as s, [(,,, oi) kh1 df = [J, va) 1., k-1 respectively, w.r.t. measure Pa with s df E RAx2", The decomposition reconstruction across scales written
s = ( n-1 n-
21 + H (1) n+1
(6)
22 = [0) ( (H()) + G
n+I
(8)
2111'
dp (1 (0(0ssn-11 +
(4)
2/411
(7)
2/ (
(
The wavelet (and also multiwavelet) transformation can be straightforwardly extended to multiple dimensions using tensor product the bases. For our purpose, a function Rd has multiscale, multiwavelet coefficients s/ dl TR.xx..xXx2"" which are also recursively obtained replacing the filters in eq. - with their Kronecker product, specifically, (0)) with ()) & H (0), where the Kronecker product repeated d times. For eq. 8 - (0(0) (and similarly others) are replaced with their dltimes Kronecker product.
The multiwavelet representation of the operator kernel K y) can be ob- tained by an appropriate tensor product of the multiscale and multiwavelet basis. One issue, however, in this approach, is that the basis at various scales coupled because of the tensor product. To untangle the basis at various scales, we use a trick as proposed I called the non-standard wavelet representation. The extra mathematical price paid for the ntn-standard representation, actually serves as a ground for reducing the proposed model complexity (see Section(2.3) thus, providing data efficiency. For the operator under consideration T with integral kernel [(x, y). let denote T,, the projection of T on V which essentially is obtained by projecting the kernel K onto basis
Non-Standard Form:

---

n))
(t(n) d
U!!)
s ((++1)
s(x)
H,C
U(L)
H G
s(L)
T
B
U(n-1) US! U U!
d(n)
A
decomposition
reconstruction
Figure 2: MWT model architecture. (Left) Decomposition cell using 4 neural networks (NNs) A, B and C. and T (for the coarsest scale L) performs multiwavelet decomposition from scale n + (Right) Reconstruction module using pre-defined filters H( G performs inverse multiwavelet transform from scale 1 to
6) w.t.t measure fAn If Pn the projection operator such that Pnf = E5. Sj. then Tn = P.TPn. Using telescopic sum, T, is expanded as
Ta T Σ 1= +..T.++ +/T11)
PrTPr.
(10)
where, Q; = P: P.-1 and L is the coarsest scale under consideration (L > 0). From eg. B is apparent that Q; the multiwavelet operator. Next, we denote A, = Q,TQ,, B, Q,TP; C2 P:-17Q:. and T Pt. In Figure we show the non-standard multiwavelet transform for given kernel (xx, 3/. The transformation has a sparse banded structure due to smoothness property of the kernel (see Section For the operator T such that Ta = the map under multiwavelet
domain is written as
UR = A, B
= U T=
(11)
where, (UT) (S) d are the multiscale, multiwavelet coefficients of u/a, respectively, and L is the coarsest scale under consideration. With these mathematical concepts, now proceed to define our multiwavelet-based operator learning model in the Section 2.3
2.3 Multiwavelet-basee Model
Based on the discussion in Section 2.2 we propose a multiwavveet- model (MWT) as shown in Figure For a given input/output as a a/u, the goal of the MWT model is to map the multiwavelet- transform of the input (s/²) output (UN) at the finest scale N. The model consists of two parts: (i) Decomposition (dec), and (ii) Reconstruction (rec). The dec acts as a recurrent network, and at the input is used obtain multiscale and multiwavelet coefficients at a coarser level s and d respectively. Next, to compute the multiscalemultivavelee coefficients of the output u, we approximate the non-standard kernel decomposition from (II) using four neural networks (NNs) A, B, C and T such that U & Ae, (df') Bag (s)), UYr R Coe (d 0A vO < < L, and U. 2 To (sf). This is ladder-down approach, and the dec part performs the decimation of signal (factor 1/2), running for maximum of L cycles, L log, (M) for given input sequence of size M. Finally, the rec module collects the constituent terms U", Un. Un (obtained using the dec module) and performs a ladder-up operation to compute the multiscale finer scale
each iteration the input is s*+ Using (6)- 7
coefficients of the output scale N obtained for the output.
using ( The iterations continue until the finest
At each iteration, the filters in dec module downsample the input, but compared to popular techniques (e.g., maxpool), the input is only transformed to a coarser multiscalemulliwavelet space. By virtue of its design, since the standaadd wavelet representation does not have inter-scale interactions, it basically allows us to reuse the same kernel NNs A, B, C' at different scales. A follow-up advantage of this approach is that the model is resolution independent, since the recurrent structure of dec is input invariant, and for a different input size M, only the number of iterations would possibly change for maximum of log, M. The reuse of A, B, C by re-training various scales also enable to
5

---

learn an expressive model with fewer parameters (fA, fg. (0,00c.,))))))) 07). single-layered CNN for A, B, C' is sufficient for learning the operator.
We see in Section 3 that even
The dec / rec module uses the filter matrices which are fixed beforehand, therefore, this part does not require any training. The model does work for any arbitrary choice of fixed matrices H, G. We show in Dooctton[.44(bbr for randomly selected matrices, the model does not learn, which validates that careful construction of filter matrices necessary.
2.4 Operators Properties
This section outlines definition of the integral kernels that are typically useful in an efficient com- pression of the operators through multiwavelets. We then discuss a fundamental property of the pseudo-differential operator.
Definition (154)). Calderón- Zyggmund Operator. The integral operators that have kernel K (x, v) which is smooth away from the diagonal, and satisfy the following.
1
[x v)1
[ y
(12)
Co
10 K(x, 10 K K(x
S
x y MTT
The smooth functions with decaying derivatives are gold to the multiwavelet transform. Note that, smoothness implies Taylor series expansion, and the multiwavelet transform with sufficiently large k zeroes out the initial & terms of the expansion due to vanishing moments property (5). This is how multiwavelet sparsifies the kernel (see Figure where K (z, y) is smooth). Although, the definition of Landerrón-Zymmund is simple (singularities only at the diagonal), but the multiwavelets are capable to compresses the kernel long the number of singularities are finite.
The next property, from [1),, points out that with input/output being single- single-dimensional functions, for any pseudo-differential operator (with smooth coefficients), the singularity at the diagonal is also well-characterized.
Property 1. Smoothness of Pseudo Pseudo-Differental Operator. For the integral kernel K(x,y) of pseudo-differente operator, K(r,y) v) E C Vx # y, and for x = K(x, v) E C where T + the highest derivative order in the given pseudo-differental equation.
The property [ iimpiees that, for the class of -eeudo-differential operator, and any set of basis with the initial J vanishing moments, the projection kernel onto such bases will have the diagonal dominating the non-diagonal entries, exponentially, if J > T
[19]. For the case of multiwavelet basis with & OPs, / = k (from eq. (5)). Therefore, k > T 1 sparsifies the kernel projection onto multiwavelets, for a fixed number bits precision We the implication of the Property on our proposed model in the Section 3.2
3 Empirical Evaluation
In this section, we evaluate the multiwavelet-based model (MWT) on several PDE datasets. We show that the proposed MWT model not only exhibits orders of magnitude higher accuracy when compared against the state-of-the-art (Sota) approaches but also works consistently well under different input conditions without parameter tuning. From numerical perspective, we take the data point-wise evaluations of the input and output functions. Specifically, with a; = a(x;), u, u(x;) for Ii, I2, ........ N E D, where x; are M M-point discretization of the domain D. Unless stated otherwise, the training set of size 1000 while test is of size 200.
have the dataset (a,, u4)
Unless otherwise stated, the NNs A, B and C in the proposed model (Figure are chosen as a single-layered CNNs following a linear layer, while T is taken as single l X l linear layer. choose = in all our experiments, and the OP basis as Legendre (Leg), Chebyshev (Chb) with uniform, non-uniform measure μlo respectively. The model in Figure is treated as single layer, and for 1D equations, we cascade 2 multiwavelet layers, while for 2D dataset, we use total 4 layers with ReLU nin-linna
Model architectures:

---

Networks
s= 64
S 128
256
512
s= 1024
MWT
0.00338
0.00375
0.00418
0.00393
0.00389
MWT Chb
0.00715
0.00712
0.00604
0.00769
0.00675
FNO
0.0125
0.0124
0.0125
0.0122
0.0126
MGNO
0.1296
0.1515
0.1355
0.1345
0.1363
LNO 0.0429
0.0557
0.0414
0.0425
0.0447
GNO
0.0789
0.0760
0.0695
0.0699
0.0721
Table I: Korteweg-de de Vries (KdV) equation benchmarks for different input resolution Top: Our methods.
Bottom:
previous works Neural operator.
Leg
M
W
uplx)
uix,1) 1]
0.0
0.2
0.4
0.5
0.8
1.0
0.0
0.2
0.4
0.6
1.0
Figure 3: The output of the KdV equation. (Left) An input (x)) with À = 0.02. (Right) The predicted
output of the MWT Leg model learning the high fluctuations.
From a mathematical viewpoint, the dec and rec modules in Figure transform only the multiscale and multiwavelet coefficients. However, the input and output to the model are point-wise function samples, i.e., (a:; u.). remedy around this is to take the data sequence, and construct hypothetical functions fa Σi a.o);, and fu Σmm u.o"; Clearly, fas Su lives in V with n = log, N.
Now model can be used with g(m) = 0; and U!!) = u, Note that fax fu are not explicitly used, but only matter of convention.
We compare our MWT model using two different OP basis (Leg, Chb) with the most recent successful neural operators. Specifically, we consider the graph neural operator (GNO) 48) the multipole graph neural operator (MGNO) 49 the LNO which makes low-rank (r) representation of the operator kernel ((x, y) (also similar to unstacked DeepONet (50), and
Benchmark models:
the Fourier neural operator (FNO) 47 We experiment on three competent datasets setup by the work of FNO (Burgers' equation (I-D), Darcy Flow (2-D), and Navier-Stokes equation (time-varying 2-D)). addition, we also experiment with Korteweg-de Vries equation (I-D). For the 1-D cases, modified FNO with careful parameter selection and removal of Batch-normalization layers results in a better performance compared with the original FNO, and we use in our experiments. The MWT model demonstrates the highest accuracy in all the experiments. The MWT model also shows the ability to learn the function mapping through lower-resolution data, and able to generalize to higher resolutions.
All the models (including ours) are trained for a total of 500 epochs using Adam optimizer with an initial learning rate (LR) of 0.001. The LR decays after every 100 epochs with factor of = 0.5. The loss function taken as relative L.2 error (7]. All of the experiments are performed on single Nvidia V100 32 GB GPU, and the results are averaged over a total of 3 seeds.
3.1 Korteweg-de Vries (KdV) Equation
The Korteweg-de Vries (KdV) equation was first proposed by Boussinesq [T6] and rediscovered by Korteweg and de Vries (3). KdV is a 1-D non-linear PDE commonly used to describe the non-linear shallow water waves. For given field u(z, t), the dynamics takes the following form:
du = -0.5u
Ou
Id u
E (0. 1), E (0.
Øt
dr
dr
(13)
(x()) = u(ær t 0)

---

The task for the neural operator is to learn the mapping of the initial condition uz(z) to the solutions
u(x,t t = 1). We generate the initial condition in Gaussian random fields according to uo (O0 7' 7 255 with periodic boundary conditions. The equation is numerically solved using chebfun package [27 with a resolution 2 and datasets with lower resolutions are obtained by sub-sampling the highest resolution data set.
Varying resolution: The experimental results of the KdV equation for different input resolutions s are shown in T We see that, compared to any of the benchmarks, our proposed MWT Leg exhibits the lowest relative error and lowest nearly by an order of magnitude. Even the case of the resolution of 64, the relative eтrоr is low, which means that a sparse data set with a coarse resolution of 64 is sufficient for the neural operator to learn the function mapping between spaces.
Varying fluctuations: We now vary the smoothness of the input function (n 0) by controlling the pa- rameter A, where low values of À imply more fre- quent fluctuations and À + O reaches the Brownian motion limit [30] To isolate the importance of incor- porating the multiwavelet transformation, we use the same convolution operation as in FNO, i.e., Fourier transform- based convolution with different modes km (only single-layer) for A, B,C. We see in Fig- ure that MWT model consistently outperforms the recent baselines for all the values of X. sample in- put/output from test set is shown in the Figure The FNO model with higher values of km has better per- formance due to more Fourier bases for representing the high-frequency signal, while MWT does better even with low modes in its A, B, C CNNs, highlight- ing the importance of using wavelet-based filters in the signal processing.
Logt
FND+10
PNFF14
L2 emor
10
0.05
0.04
0.03
0.02
Figure 4: Comparing MWT by varying the degree of fluctuations À in the input with resolution 1024. For each convolution, we fix the number of Fourier hases as km. For FNO, the width is 64.
3.2 Theoretical Properties Validation
We test the ability of the proposed MWT model to capture the theoretical properties of the pseudo- differential operator in this Section. Towards that, we consider the Guler-Bernoull equation [62] that models the vertical displacement of finite length beam over time. A Fourier transform version of the beam equation with the constraint of both ends being clamped is as follows
8²u Ox
du
2 = fx),
dx
===
=
(14)
m(0) = (01)=) 0,
where u(x)) is the Fourier transform of the time-varying beam displacement, w is the frequency, f(x) is the applied force. The Buler-Bernoulli is a pseudo-differential equation with the maximum derivative order T + 4. We take the task of learning the map from I to u. Figure we see that for & 2 3, the models relative error across epochs is similar, however, they are different for & < 3, which is in accordance with the Property the diagonal of the kernel which is CT hennce, sparsification cannot occur, and the model learns slow.
[ For < 3, the multiwavelets will not be able to annihilate
3.3 Burgers' Equation
The 1-D Burgers' equation is a non-linear PDE occurring in various areas of applied mathematics. For given field u(x, 1) and diffusion coefficient v, the 1-D Burgers' equation reads:
du = -
Ou
y2
+ v
E (0,2 e (0.
Ōt uo(x) = u(x, t = 0).
dr
Oz
(15)
The task for the neural operator to learn the mapping of initial condition u(ar, 1. = 0) to the solutions at t = 1 u(x, = 1). To compare with many advanced neural operators under the same conditions,
8

---

10
10
Relative error
Relative error
10
10
MWT Clb
10
100
200
300
400
500
10
epochs
256
2048
4096 Resolution (8)
8192
Figure 5: Relative L2 error vs epochs for MWT Leg with different number of OP basis & L,.... 6.
Figure 6: Burgers' Equation validation various input resolution 8. Our methods: MWT Leg, Chh.
Networks
$
S 64
$ 128
$ 256
$ 512
MWT Leg
0.0152
0.00899
0.00747
0.00722
0.00654
MWT Chb
0.0174
0,0108
0.00872
0,00892
0,00891
MWT Rnd
0,2435
0.2434
0.2434
0.2431
0.2432
FNO
0.0177
0.0121
0.0111
0,0107
0.0106
MGNO
0.0501
0.0519
0.0547
0.0542
LNO
0.0524
0.0457
0.0453
0.0428
Table 2: Benchmarks on Darcy Flow equation at various input resolution 8. Top: Our methods. MWT Rnd
instantiate random entries of the filter matrices in 6 (.. Bottom: prior works on Neural operator.
we use the Burgers' data and the results that have been published in 47 and [49]. The initial condition is sampled as Gaussian random fields where uo ~ (0 5+ 52T)-2) with periodic boundary conditions.
is the Laplacian, meaning the initial conditions are sampled by sampling its first several coefficients from Gaussian distribution. In the Burgers' equation, is set to 0.1. The equation is solved with resolution 213, and the data with lower resolutions are obtained by sub-sampling the highest resolution data set.
The results of the experiments on Burgers' equation for different resolutions are shown in Figure Compared to any of the benchmarks, our MWT Leg obtains the lowest relative error, which is an order of magnitude lower than the state-of-the-an resolution, MWT Leg still maintains very low error rate, which shows its potential for learning the function mapping through low-resolution data, that is, the ability map between ddmiite-dimensionall spaces by learning limited finiteddimensional spaces mapping.
It's worth noting that even in the case of low
3.4 Darcy Flow
Darcy flow formulated by Darcy 22 one of the basic relationships of hydrogeology. describing the flow of a fluid through a porous medium. We experiment on the steady-state of the 2-d Darcy flow equation on the unit box, where it takes the following form:
V Vu(x)) f(x),
(0,
(16)
u(xx)
E 1)
We set the experiments to learn the operator mapping the coefficient a(x) to the solution u(x)) The coefficients are generated according to a ~ (0, 321) where is the Laplacian with zero Neumann boundary conditions. The threshold of a(x) is set to achieve ellipticity. The solutions u(z) are obtained by using a 2nd-order finite difference scheme on a 512 x 512 grid. Data sets of lower resolution sub-sampled from the original data set.
The results of the experiments on Darcy Flow for different resolutions are shown in Table MWT Leg again obtains the lowest relative error compared to other neural operators at various resolutions.

---

We also perform an additional experiment, in which the multiwavelet filters H GO, = 0, 1 are replaced with random values (properly normalized), We see in Table that MWT Rnd does not learn the operator map, in fact, its performance is worse than all the models. This signifies the importance of the careful choice of the filter matrices.
3.5 Additional Experiments
Full results for these experiments are provided in the supplementary materials.
The Navier-Stokes (NS) are 2d time-varying PDEs modeling the viscous, incompressible fluids. The proposed MWT model does a 2d multiwavelet transform for the velocity u, while uses single-layered 3d convolution for A, B and C' to learn dependencies across space-time. We have observed that the proposed MWT Leg is in par with the Sota on the NS equations.
Navier Stokes Equation:
Prediction at high resolution: We show that MWT model trained at lower resolutions for various datasets (for example, training with e = 256 for Burgers) can predict the output at finer resolutions 2048, with relative error of 0.0226, thus eliminating the need for expensive sampling. The training and testing with = 2048 yields a relative error of 0.00189.
We study the operator learning behavior when the training and evaluation datasets are obtained using the random function from different generating rules. The training is done with squared exponential kernel but evaluation is done on different generating rule with controllable parameter
Train/evaluation with different sampling rules:
Conclusion
We address the problem of data-driven learning of the operator that maps between two function spaces. Motivated from the fundamental propertics the integral kernel, we found that multiwavelets constitute natural basis to represent the kernel sparsely. After generalizing the multiwavelets to work with arbitrary measures, we proposed a series models to learn the integral operator. This work opens up research directions and possibilities toward designing efficient Neural operators utilizing properties of the kernels, and the suitable basis. We anticipate that the study of this problem will solve many engineering and biological problems such as aircraft wing design, complex fluids dynamics, metamaterials design, cyber-physical systems, neuron-neuron interactions that are modeled by complex PDEs.
Acknowledgement
We are thankful to the anonymous reviewers for providing their valuable feedback which improved the manuscript. We would also like thank Radu Balan for his valuable feedhack. We grate- fully 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-193261 the U.S. Army Research Office (ARO) under Grant No. 1091NNF-1-11-000 the Okawa Foundation award, and 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 and a Northrop Grumman grant. A part of this work used the Extreme Science and Engineering Discovery Environ- ment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. The views, opinions, and/or findings contained in this article are those of the authors and should not be interpreted representing the official views policies, either expressed implied by the Defense Advanced Research Projects Agency, the Army Research Office, the Department of Defense or the National Science Foundation.
10

---

References
[1] M. Abramowitz and L.A Stegun. Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables. Applied mathematies series. Dover Publications, 1965. ISBN 9780486612720.
[2] Jonas Adler and Ozan Öktem. Solving ill ill-posed inverse problems using iterative deep neural networks. Inverse Problems, 33(12), Nov 2017. ISSN 1361-6420, doi: 10.1088/1361-6420 aa9581. URL http://dx. doi. org/ 10. 0015//66616642//a22/222
[3] B. Alpert, G. Beylkin, R. Coifman, and V. Rokhlin. Wavelet-like bases for the fast solution of second-kind integral equations. SIAM Journal on Scientific Computing, 14(1):159-184, 1993.
doi: 10.1 (3370914405)
[4] B. Alpert, G. Beylkin, D. Gines, and L. Vozovoi. Adaptive solution of partial differential equationns in multiwavelet bases. Journal of Computational Physics, 182(1):149-199, 2002. ISSN 0021-9991.
[5] Bradley K. Alpert. class of bases in I for the sparse representation of integral operators. SIAM Journal on Mathematical Analysis, 24(1):246 262,, 1993. doi: 10. 10.1137//52401.
[6] 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
[7] Kevin Amaratunga and John Williams. Wavelet based green's 's function approach to 2d pdes. Engineering Computations, 10, 07 2001. doi: 10. 188/bb033333
[8] Radu Balan, Bernhard G. Bodmann, Peter G. Casazza, and Dan Edidin. Painless reconstruction from magnitudes of frame coefficients. Journal of Fourier Analysis and Applications, 16: 488-501, 2009.
[9] J.J. Benedetto. Wavelets: Mathematics and Applications.. Studies in Advanced Mathematics. Taylor & Francis, 1993. ISBN 9780849382710.
[10] G. Beylkin. On the representation of operators in bases of compactly supported wavelets. SIAM Journal on Numerical Analysis, 29(6): 176--774,, 1992. doi: 10. 13700729097
[11] G. Beylkin, R. Coifman, and V. Rokhlin. Fast wavelet transforms and numerical algorithms
i. Communications Pure and Applied Mathematics, 14(2):141-183, 1991. doi: https:
[12] Gregory BeylkIn and James M. Keiser. An adaptive pseudo-wavele approach for solving nonlinear partial differential equations. In Multiscale Wavelet Methods for Partial Differential Equations, volume 6 of Wavelet Analysis and Its Applications, pages 137-197. Academic Press, 1997.
[13] 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. 10007//00446-009-00114
[14] Kaushik Bhattacharya, Bamdad Hosseini, Nikola B. Kovachki, and Andrew M. Stuart. Model reduction and neural networks for parametric pdes, 2020.
[15] Nicolas Boulle, Christopher J. Earls, and Alex Townsend. Data-driven discovery of physical laws with unddrsstannee deep learning, 2021.
[16] Joseph Boussinesq. Essai sur la théorie des eaux courantes. Impr. nationale, 1877.
[17] R. Y. Chang and M. L. Wang. Shifted legendre direct method for variational problems. Journal of Optimization Theory and Applications, 39(2):299-307, Feb 1983. ISSN 1573-2878. doi: 10. 1000/BB009945355
[18] Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic recurrent neural networks. In International Conference on Learning Representations, 2020. URL https://openrevieet
[19] Kenneth C. Chou and Gary S. Guthart. Representation of green's 's function integral operators using wavelet transforms. Journal of Vibration and Control, 6(1):19-48, 2000. doi: 10.1 10775463000060010..
11

---

[20] N.M. Chuong. Psendodiffereniaal Operators and Wavelets over Real and p-adic Fields. Springer International Publishing, 2018. ISBN 9783319774732.
[21] M. Cotronei, L.B. Montefusco, and Puccio. Multiwavelet analysis and signal processing. IEEE Transactions Circuits and Systems II: Analog and Digital Signal Processing, 45(8): 970-987, 1998.
[22] Henry Darey. Les fontaines publiques de la ville de Dijon: exposition et application... Victor Dalmont, 1856.
[23] Olivier Darrigol. Worlds flow: history of hydrodynamics from the Bernoullis to Prandil. Oxford University Press, 2005.
[24] Ingrid Daubechies. Orthonormal bases of compactly supported wavelets. Communications on Pure and Applied Mathematics, 41(7):909 996, 1988. doi: http:://doi.or/10.1000022
3160410705.
[25] Ingrid Daubechies. Ten Lectures Wavelets. Society for Industrial and Applied Mathematics, 1992. doi:
[26] D. Deng, Y. Meyer, and Y. Han. Harmonic Analysis on Spaces of Homogeneous Type. Lecture Notes in Mathematics. Springer Berlin Heidelherg, 2008. ISBN 9783540887447.
[27] A Driscoll, N. Hale, and L. N. Trefethen. Chebfin Guide. Pafnuty Publications, 2014.
[28] Yuwei Fan, Jordi Feliu-Fabss Lin Lin, Lexing Ying, and Leonardo Zepeda- multiscale neural network based on hierarchical nested bases, 2019.
[29] Jordi Feliu-Fabà, Yuwei Fan, and Lexing Ying. Meta-learning pseudo-differente operators with deep neural networks. Journal of Computational Physics, 408:109309, May 2020. ISSN 0021-9991. doi: 10. 10.1016/j.jcp.2222 20w...090...00.,
[30] Silviu-loan Filip, Aurya Javeed, and Lloyd Trefethen. Smooth random functions, random odes, and gaussian processes. SIAM Review, 11:18550 01 2019. doi: 10.1 1
[31] Daniel Greenfeld, Meirav Galun, Ron Kimmel, Irad Yavneh, and Ronen Basri. Learning to optimize multigrid pde solvers, 2019.
[32] Albert Gu, Tri Dao, Stefano Ermon, Atri Rudra, and Christopher Ré. Hippo: Recurrent memory with optimal polynomial projections. In Advances in Neural Information Processing Systems, volume 33, pages 1474-1487, 2020.
[33] Xiaoxiao Guo, Wei Li, and Francesco Iorio. Convolutional neural networks for steady flow ap- proximation. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD '16, page 481-490. Association for Computing Machinery, 2016.
[34] Mustafa Hajij. Ghada Zamzmi, Matthew Dawson, and Greg Muller. Algebraicalll-informe deep networks (aidn): A deep learning approach to represent algebraic structures, 2021.
[35] M.H. Heydari, M.R. Hooshmandasl, and F. Mohammadi. Legendre wavelets method for solving fractional partial differential equations with dirichlet boundary conditions. Applied Mathematics and Computation, 234:267-276, 2014. ISSN 0096-3003, doi: https://doi.org// 10. 10166/..amm 2014.02..07.
[36] Chyi Hwang and Yen-Ping Shih. Solution of integral equations via laguerre polynomials. Computers & Electrical Engineering, 933)1123-12, 1982. ISSN 0045-7906. doi: https:
[37] Ameya D. Jagtap, Yeonjong Shin, Kenji Kawaguchi, and George Em Karniadakis. Deep kronecker neural networks: A general framework for neural networks with adaptive activation functions, 2021.
[38] M. Tavassoli Kajani, A. Hadi Vencheh, and M. Ghasemi. The chebyshev wavelets opera- tional matrix of integration and product operation matrix. International Journal of Computer Mathematics, 86(7):1118 (177::111-11111 1125, 2009. doi: 10. 00.08800000111600117666666
[39] F. Khellat and S.A. Yousefi. The linear legendre mother wavelets operational matrix of inte- gration and its application. Journal of the Franklin Institute, 348(2):181-199, 2006. ISSN 0016-0032. doi: https://doiorr///
12

---

[40] Yuehaw Khoo, Jianfeng Lu, and Lexing Ying, Solving parametric pde problems with artificial neural networks. European Journal of Applied Mathematics, 32(3):421-433 Jul 2020. ISSN 1469-4425.
[41] Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equations for irregular time series, 2020.
[42] Dmitrii Kochkov, Jamie A. Smith, Ayya Alieva, Qing Wang, Michael P. Brenner, and Stephan Hoyer. Machine learning- computational fluid dynamics. Proceedings of the Na- tional Academy of Sciences, 118(21), 2021. ISSN 0027-8424. doi: 10. 2..107/pnas.21017885 18.
[43] Vasil Kolev, Todor Cooklev, and Fritz Keinert. Design of simple orthogonal multiwavelet filter by matrix spectral factorization. Circuits, Systems, and Signal Processing, 39(4)2006--2411 Aug 2019.
[44] Samuel Lanthaler, Siddhartha Mishra, and George Em Karniadakis. Error estimates for deep- onets: A deep learning framework in infinite dimensions, 2021.
[45] Yuanlu LI. Solving a nonlinear fractional differential equation using chebyshev wavelets. Communicationss in Nonlinear Science and Numerical Simulation, 15(9):2284-222 2010.
ISSN 1007-5704. doi: http:://doi.or//I 10.1066/..cnsss.20 09.0000.
[46] Yuanlu Li and Weiwei Zhao. Haar wavelet operational matrix of fractional order integration and its applications in solving the fractional order differential equations. Applied Mathematics and Computation, 216(8):2276 2285,, 2010. ISSN 0096-3003. doi: https://doi.org/10.100
[47] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differen- tial equations, 2020.
[48] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations, 2020.
[49] 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, pages 6755-6766, 2020.
[50] 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.
[51] Lu Lu, Pengzhan Jin, Guofei Pang. Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation thcorem of operators. Nature Machine Intelligence, 3(3):218-22) Mar 2021. ISSN 2522-5839.
[52] Ryan McKeown, Rodolfo Ostilla-Mónico, Alain Pumir, Michael P. Brenner, and Shmuel M. Rubinstein. Turbulence generation through an iterative cascade of the elliptical instability. Science Advances, 6(9), 2020. doi: 10. 126/sciadv.aar2211..
[53] Y. Meyer and D.H. Salinger. Wavelets and Operators: Volume 1. Cambridge Studies in Advanced Mathematics. Cambridge University Press, 1992. ISBN 7780521458696
[54] Y. Meyer, R. Coifman, and D. Salinger. Wavelets: Calderón-Zygmann and Multilinear Opera- tors. Cambridge Studies in Advanced Mathematics. Cambridge University Press, 1997. ISBN 9780521420013.
[55] Guillaume Ovarlez, Anh Vu Nguyen Le, Wilbert J. Smit, Abdoulaye Fall, Romain Mari, Guillaume Chatté, and Annie Colin. Density waves in thear-thic suspensions. Science Advances, 6(16), 2020. doi: 10. 10.1126/sciddv.aa55588
[56] 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
13300.
[57] 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
13

---

differential equations. Journal of Computational Physics, 378:686- 707, 2019, ISSN 0021- 9991.
[58] C.E. Rasmussen and C.K.I. Williams. Gaussian Processes for Machine Learning. Adaptative computation and machine learning series. University Press Group Limited, 2006.
[59] Mohsen Razzggii and Samira Yousefi. The legendre wavelets operational matrix of integration. International Journal of Systems Science, 32:495-502, 04 2001. doi: 10. 10000000000022212227,
[60] M. F. Shlesinger, B. J. West, and J. Klafter. Lévy dynamics of enhanced diffusion: Application to turbulence. Phys. Rev. Lett., 58:1 100-1101: 103, Mar 1987. doi: 10. P..110//yyseee..et
100.
[61] B. W. Silverman, C. Vassilicos, and Nick Kingsbury. Image processing with complex wavelets. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 357 (57(1766))54332556 1999. doi: (0.1009//sta, 1999.0447.
[62] S. Timoshenko. History of Strength of Materials: With Brief Account of the History of Theory of Elasticity and Theory of Structures. Dover Civil and Mechanical Engineering Series. Dover Publications, 1983. ISBN 9780486611877.
[63] J. Towns, T. Cockerill, M. Dahan, I, Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, Peterson, R. Roskies, R. Scott, and N. Wilkins-Diehr. Xsede: Accelerating scientific discovery. Computing in Science & Engineering, 1665)x62-7, Septt-OO 2014. ISSN 1521-9615. doi: 20.110//CSE.2014..80 URL doi jew.ccccpptrrrocceety.ooo.ooooo. UCSE.2014.80 80
D. Lifka,
[64] Lifeng Wang, Yunpeng Ma, and Zhijun Meng. Haar wavelet method for solving fractional partial differential equations numerically. Applied Mathematics and Computation, 227:66-76,
2014. ISSN 0096-3003. doi: https://doi.or//1 1tt.:i/000/..11............ 101.1.06....mc2201 111.3.11.04...
[65] Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with physics-informed deeponets, 2021.
[66] Yanxin Wang and Qibin Fan. The second kind chebyshev wavelet method for solving fractional differentiaa equations. Applied Mathematics and Computation, 218(17):8922-601, 2012. ISSN 0096-3003. doi: https://doi.org/ ortp:/..//////0000000 www.///..m....12.......
[67] Shiying Xiong, Xingzhe He, Yunjin Tong, and Bo Zhu. Neural vortex method: from finite lagrangian particles to infinite dimensional eulerian dynamics, 2020.
[68] Yuankun Xue and Paul Bogdan. Constructing compact causal mathematical models for complex dynamics. In 2017 ACM/IEEE 8th International Conference on Cyber-Physicad Systems (ICCPS), pages 97-108, April 2017.
[69] Li Zhu and Yanxin Wang. Solving fractional partial differential equations by using the second chebyshev wavelet operational matrix method. Nonlinear Dynamics, (915):115-1997 Aug 2017. ISSN 1573-269X. doi: 10.1007/s1 107 1070077111331-11771
[70] Li Zhu and Yanxin Wang. Solving fractional partial differential equations by using the second chebyshev wavelet operational matrix method. Nonlinear Dynamics, 1915:1915-1922 Aug 2017. ISSN 1573-269X. doi: 10. 1007/s1 1071-017-33-111111
[71] Yinhao Zhu and Nicholas Zabaras. Bayesian deep convolutional secoder-decodde networks for surrogate modeling and uncertainty quantification. Journal of Computational Physics, 366: 415-447, Aug 2018 ISSN 0021-9991.
14

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] For example, see Table Figure for benchmarks on the datasets. Also, see Figure for robustness plot, and Figure for theoretical insights for pseudo-differential operators.
(b) Did you describe the limitations of your work? [Yes] We discuss in details the possible numerical issues that can occur in estimating the filter matrices H _I(1), ( G
for large values of k. The issue is not related to the mathematics involved but due to the nature of floating point precision. We discuss this in details the supplementary materials.
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] Please refer to Section where we list all the necessary derived results, while we have referred the reader (at appropriate places) to the supplementary materials for the complete details.
(b) Did you include complete proofs of all theoretical results? [Yes] We provide detailed derivations of all the measure dependent filter matrices H H()), GOO), G() and also the correction terms (()) [ in the supplementary materials.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed reproduce the main exper- imental results (either in the supplemental material or as a URL)? [Yes] The code is uploaded with the supplementary materials.
(b) Did you specify all the training details (e.g.. g. data splits, hyperparameters, how they were chosen)? [Yes] Please refer to Section 3 where we list all of the details regarding training and model architectures.
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes] All the results included in the paper are averaged over total of seeds. We have also mentioned the same in Section
(d) Did you include the total amount of compute and the type of resources used (e.g.. type of GPUs, internal cluster, or cloud provider)? [Yes] All of the experiments were performed on a single Nvidia V100 32 GB GPU, please refer to Section
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing asscts, did you cite the creators? [Yes] A part of the datasets are taken from the FNO work [47], while some are generated using the scripts provided by the same authors. We have properly cited the work in Section Benchmark models.
(b) Did you mention the license of the assets? [N/A] part of the code and datasets (47)) used by us are openly available with no license restriction, to the best of our knowledge.
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]
(d) Did you discuss whether and how consent was obtained from people whose data you'l 're using/curating? [Yes] The code and the dataset is openly available.
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] All the datasets are synthetically generated.
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
15


Reviews:
[{"review_id": "cR_541hzwQU", "reviewer": "Reviewer_RXFA", "summary": "This paper studies the problem of approximating the solution operator of a partial differential equation (PDE) that maps the inputs to the solution. The authors begin by reviewing the deep learning techniques that have been recently introduced and the multiple applications of this problem. The main contribution of the paper is to propose a new neural network architecture, based on the multiwavelet transform, to approximate the kernel of the solution operator. Finally, the authors compare their method against several Neural operator techniques (Fourier neural operator, graph neural operator, ...) on a number of numerical PDE examples and achieve lower relative error (2X - 4X).", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "## Originality\n\nThe problem of approximating the solution operator is not novel and many techniques have been recently proposed, which the authors reference properly in the introduction. A deep learning approach based on the wavelet transform has also been proposed by [1] and I believe that the authors essentially generalize this approach to employ an orthonormal polynomial basis (Legendre or Chebyshev basis) as basis for the wavelet construction. The paper would gain in quality by explaining the advantages and drawbacks of this approach compared to [1]. The dimension of the subspace generated is larger, which gives higher approximation power, but does it come with higher computational training time for the networks ?\n\n[1] Feliu-Faba, Fan, Ying, J. Comput. Phys., 2020.\n\n## Clarity\n\nThe introduction of the paper is well-written and summarizes the potential applications as well as the recent works on the topic. However, I found several places where the paper lacks clarity, beginning with section 2. In the problem setup, the problem is not clearly defined and I do not see what the function a(x) represents for the PDE: is it a forcing term, a coefficient function in the PDE ? \nSimilarly, the authors \"do not put the restriction of linearity\" on the operator T, but represent it as an integral kernel and try to approximate the kernel. This formulation won't hold for the nonlinear PDEs that the authors are considering later.\n\nSection 3.4 discusses standard properties of integral kernels but the definition and property cited are not mathematically rigorous as they omit the definition of the objects or the spatial domain. Furthermore, Property 1 does not hold for the Green's function of the Laplace operator in spatial dimension 3 (the kernel is singular and not continuous along the diagonal as stated by the property).\nThe experimental setup at the beginning of section 3 also lacks clarity: do the authors learn the operator for 1 input-output pair (a(x),u(x)), evaluated at N points of the domain (i.e. one-shot learning), or do they choose several input-output pairs, which they then sample on the domain ?\n\nMost of the section A of the supplementary material is very standard and the authors could just reference the relevant results later instead of redefining Legendre, Chebyshev polynomials or describing Gram-Schmidt orthogonalization.\n\n## Quality\n\nThe authors propose to use a compressed representation of the integral kernel (or Green's function) because learning a direct representation of the kernel is expensive in high dimensions (essentially because the kernel has twice the spatial dimension of the solutions to the PDE).\n\n However, it is feasible to learn the kernel for the problems that the authors are considering in the paper (spatial dimension 1,2, and simple square geometries) [2,3]. Moreover, even for simple problems (Poisson equation in 1D/2D), the kernel can be challenging to learn due to the singularity along the diagonal. I'm wondering whether the compressed approach discussed in the paper would be able to accurately learn the Green's function in this setting and outperform the methods proposed by [2,3] by comparing the L2 norm of the kernel function. \nIn  fact, using a testing set of functions to benchmark the method can be misleading, especially if the testing set is generated from the same distribution as the training set, because the method could simply overfit the training functions and not be accurate for any type of input functions.\n\nSimilarly, the choice of the Gaussian random fields for the input functions (e.g. $f \\sim \\mathcal{N}(0, 7^4(-\\Delta+7^2 I)^{-2.5}$)) seems a bit strange and varies across the different examples. Since the method is also assessed on functions generated by the same distribution, it is also a possibility that the method performs well for this particular choice but fails for any other choice. \nWhy did the authors not consider the standard squared-exponential covariance kernel ?\n\nFinally, the benchmarks are reported on an average of 3 experiments, but the authors did not report the individual errors, which questions the stability of the technique with respect to multiple runs and the standard deviation of the errors achieved.\n\n[2] Li et al., arxiv:2003.03485, 2020.\n\n[3] Boulle, Earls, Townsend, arxiv:2105.00266, 2021.\n\n## Significance\n\nThe gain in the relative error (2X-4X) achieved by the author is moderately significant but is compensated by the different flaws in the numerical experiments pointed out in the previous section. This neural network architecture could be of practical interest, provided it is either computationally faster than existing techniques or more accurate (but validated by more robust numerical experiments). The theoretical justification of the method is a bit vague due to the lack of mathematical rigour in sections 2.1-2.4 and probably does not apply for the range of examples that the authors are considering.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I do not see any negative societal impact in this study. One limitation that the authors could discuss is that the method currently works on simple tensor-products domain, on which simple polynomial bases exist. I do not know if the method naturally extends to more complex geometries or non-uniform spatial discretization.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4 hours"}, "extra_text": {}}, {"review_id": "QxqUtTC3i2i", "reviewer": "Reviewer_z5a1", "summary": "This paper proposes to use a multiwavelet representation in order to learn operators which compute non-linear PDEs. The input and output functions are represented in multiwavelet bases.  The key idea is that multiwavelets provide sparse representation of large class of Calderon-Zygmund operators with a \"non-standard form\" which keep both the wavelet and scaling function coefficients, and are thus used in numerical schemes to solve many PDEs such as Burger equations or KdV.  In these cases, the specific non-linear operator corresponding to a PDE is a non-linear map between two sparse representations and are thus expected to be easier to learn then from the original time or spatial functions. The operator learning is performed with neural networks applied at each scale. The algorithm is resolution independent because the resolution parameter is incorporated in the multiwavelet basis. Results are demonstrated over Burger, Korteweg-de Vries and Darcy's flow with better than state of the art results. ", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper is presenting an interesting approach to learn non-linear operators arising from PDEs. By using prior information on the nature of the solutions and the underlying operators, it  builds a sparse representations which simplifies the estimation problem. This type of work at the interface between harmonic analysis representations of PDEs solutions and neural networks is interesting and can potentially lead to real applications. At this point, the problems that is studied are very academic (one-dimensional Burger or KdV or Darcy flows) but they are none the less interesting. Whereas the paper spends a lot of time on the multiwavelet representation and its properties to represent Calderon-Zygmund integral operators, which are known in harmonic analysis, they spend little space on the description of the neural networks that perform the learning, the choice of network and an analysis of its properties. It would be very useful to see how the complexity of the network relates to the complexity of the underlying learned operator once represented in a sparse non-standard form. Overall this is a good paper that makes an interesting link between PDE and learning. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4h"}, "extra_text": {}}, {"review_id": "HPwOK3aWHqu", "reviewer": "Reviewer_rY49", "summary": "The work develops Multiwavelet-based neural operator using multiwavelet Transform and orthogonal polynomials. The model applies multiscale transfer and kernel integration. Empirically the model is more accurate compare to Fourier neural operator on many PDEs.", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "In my opinion, this work is a natural and novel extension to neural operators. The architecture is very solid and the performance is quite good. I would vote for strong acceptance.\n\nConcerns:\n(1) It's understandable that most of the benchmark numbers are copied from the previous papers. For the Burgers' equation, it was reported as 0.0149 on FNO but if the batch norm layers are taken out, the error on FNO will be ~ 0.002. This can also be true for the Korteweg-de Vries equation. This may explain why the gap of MWT and FNO is so larger on Burgers and Korteweg-de Vries, but not so much on Darcy and Navier-stokes. \nTo have a fair and solid comparison, I strongly suggest the author use the tuned FNO model and re-run for Burgers and Korteweg-de Vries equation.\n(2) MWT has a multi-scale structure, which may imply the overall runtime of MWT can be slower. It will be great if the author could report the runtime of MWT and how it compares to the baseline models.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "It will be better if the author can discuss when wavelet does not apply. To my knowledge, while there are many works on wavelets, but in general people still prefer FEM/FDM solvers, etc. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "BDkp6iOoqEn", "reviewer": "Reviewer_aK5Q", "summary": "The manuscript entitled “Multiwavelet-based Operator Learning for Differential Equations” exploits the fundamental properties of the operator’s kernel which enable numerically efficient representation. The proposed method exhibits a 10X improvement on relative L2 error for some of the PDEs compared with the state-of-the-art neural operators, which is very impressive. The proposed method also shows the ability to find the solution of a high-resolution input after learning from lower-resolution data.", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper will have a significant impact in the field of operator learning for differential equations. I am looking forward to its publication with some modifications to address the below issues:\n\n1.\tOn Line 115 Page 3, should it be V(0, k) instead of V(0, n)?\n\n2.\tIn Figure 1 on Page 4, the description for this figure is not clear. It is hard to understand (ii) and (iii) by the caption. It doesn’t even mention (iii) explicitly. What is the relation between H(0), G(0), H(1), G(1) the 3d matrix displayed in the figure? On the last line of the caption, it should be i = 0, 1 instead of i = 1, 2.\n\n3.\tOn Line 178 Page 5, it is not clear what the predefined filter matrices look like. Only looking at Figure 2 and the contexts in the 2.3 Multiwavelet-based Model, it is very hard for the readers to build the same model on their side, hence it is hard to reproduce the experimental results in the paper. The authors should plot another graph for the neural networks so readers can understand how to build the neural network in practice. The author should show at least one detailed neural network for one of its PDE experiments for readers to reproduce the results.\n\n4.\tI don’t see the supplementary materials.\n\nIn summary, I recommend its publication in NeurIPS after those commends are addressed.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Authors should address the limitation of the filter matrices used in this paper. Should they apply to all PDE problems? Why does this approach result in 10x accuracy improvement for some PDEs while only gives about 1x accuracy improvement (no improvement) for the Navier-Stokes (NS) are 2d time-varying PDEs?", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "19Br4TjPhIi", "reviewer": "Reviewer_SQX4", "summary": "This paper presents a multiwavelet-based neural operator learning scheme that compresses the associated operator’s kernel using fine-grained wavelets. Using the proposed model, the paper performs experiments on the Korteweg-de Vries (KdV) equation, Burger’s equation, Darcy Flow, and Navier-Stokes equation. Compared with the existing neural operator approaches, this model claims to show significantly higher accuracy and can find the solution of a high-resolution input after learning from lower-resolution data. This paper appears to be a novel proposal, and certainly of interest. However, the method only applies to some toy examples, its application prospect seems narrow.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths\n\n\nThis paper operates in a general area (neural PDE) of great interest to the machine learning community. The method demonstrates the applicability on higher dimensions of the 2-D Darcy flow equation in Section 3.4, which is a good point. Also, I appreciate the detailed supplementary material provided by the authors.\n \nWeaknesses\n\n\n- The introduction of the theory seems a little indirect and not concise, such as Section 2.1 is a bit redundant.\n- Although the paper shows a lot of error curves, the results are rarely shown. To enhance readers' intuition of the predictive power of the presented model, the temporal evolution of the PDE solutions should be visualized.\n- The result comparison is not convincing in its current form. It's described in a fairly vague way, and there's a lack of results showing that it works. This part would need significant extensions, or potentially could be left out.\n- Table 3 in supplementary material looks like the simulation error is too large to be practical.\n\n\nCorrectness\n\n\nThe technical content of this paper appears to be correct.\n \n\nClarity\n\n\nOverall, the exposition of the paper is fairly clear. But some details are confusing:\n- MWT model architecture in Figure 2 is confusing. I can't follow the data flow of the network.\n- What are the two subfigures in Figure 3?  Please specify the meaning of these two figures.\n \n\n\nRelation to prior art\n\n\nThis paper compares its MWT model using two different OP bases (Leg, Chb) with the most recent successful neural operators, such as the graph neural operator (GNO), the multipole graph neural operator (MGNO), the LNO which makes a low-rank representation of the operator kernel K(x; y), and the Fourier neural operator (FNO ). The MWT model claims to show the ability to learn the function mapping through lower-resolution data, and be able to generalize to higher resolutions. These methods are compared in the results section, but the demonstration seems to need to be strengthened. Also, it would be helpful to add some citations of recent Neural ODE/PDE solvers, such as papers related to conservation laws, symplectic neural networks, and neural vortex methods.\n \nReproducibility\n\n\nThe method description is clear in general and reasonably easy to follow.\n\n\n \nAdditional Feedback/Questions\n\n\nHow to handle the modal decomposition of a problem with nontrivial boundary conditions?\n \n\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I would love to see some limitations of the proposed method. In addition, the method should discuss more robustness to real noise data.", "review_text": "", "extra_scores": {"ethical_concerns": "None", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper investigates the problem of approximating the solution operator of a partial differential equation (PDE) and proposes a novel neural architecture based on the multi-wavelet transform to approximate the kernel of the solution operator. The key innovation of this paper is a clever application of the multi-wavelet transform that allows for a simpler solution. There is a general consensus among the expert reviewers that this paper is of high quality and offers a new solution to the important problem that would have a high impact on the PDE community. Some reviewers raise a criticism that the paper didn't provide the results on complex real-world scenarios but focus instead on the simple settings. I consider this weakness to be minor for the NeurIPS audience, especially after the authors have satisfactorily described potential directions to improve this work in their responses.

As a result, I recommend this paper for publication at NeurIPS2021.
Decision: Accept (Spotlight)