MgNO: Efficient Parameterization of Linear Operators via Multigrid
Authors: Juncai He
Xinliang Liu
Jinchao Xu
Venue: ICLR 2024 poster (2024)
Keywords: neural operator; multigrid; universal approximation; boundary condition
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: In this work, we propose a concise neural operator architecture for operator learning. Drawing an analogy with a conventional fully connected neural network, we define the neural operator as follows: the output of the $i$-th neuron in a nonlinear operator layer is defined by $\mathcal O_i(u) =  \sigma\left( \sum_j \mathcal W_{ij} u + \mathcal B_{ij}\right)$. Here, $\mathcal W_{ij}$ denotes the bounded linear operator connecting $j$-th input neuron to $i$-th output neuron, and the bias $\mathcal B_{ij}$ takes the form of a function rather than a scalar. Given its new universal approximation property, the efficient parameterization of the bounded linear operators between two neurons (Banach spaces) plays a critical role. As a result, we introduce MgNO, utilizing multigrid structures to parameterize these linear operators between neurons.  This approach offers both mathematical rigor and practical expressivity. Additionally, MgNO obviates the need for conventional lifting and projecting operators typically required in previous neural operators. Moreover, it seamlessly accommodates diverse boundary conditions. Our empirical observations reveal that MgNO exhibits superior ease of training compared to CNN-based models, while also displaying a reduced susceptibility to overfitting when contrasted with spectral-type neural operators. We demonstrate the efficiency and accuracy of our method with consistently state-of-the-art performance on different types of partial differential equations (PDEs).

Published as a conference paper at ICLR 2024
MGNO: EFFICIENT PARAMETERIZATION OF LINEAR OPERATORS VIA MULTIGRID
Juncai He, Xinliang Liu JJinhhoo Xu King Abdullah University of Science and Technology (juncai he, xinliang. 1iu, jinchao www.ooust.kkk...o.o edu. sa
ABSTRACT
In this work, we propose concise neural operator architecture for operator learn- ing. Drawing an analogy with conventional fully connected neural network, we define the neural operator follows: the output of the i-th neuron in a nonlinear operator layer is defined by O,(n) or (Σ, Wij + Bi)). Here, Wij denotes the bounded linear operator connecting j-th input neuron i-th output neuron, and the bias Byy takes the form of a function rather than scalar. Given its new universal approximation property, the efficient parameterization of the bounded linear oper- ators between two neurons (Banach spaces) plays critical role. As result, we introduce MgNO, utilizing multigrid structures to parameterize these linear opera- tors between neurons. This approach offers both mathematical rigor and practical expressivity. Additionally, MgNO obviates the need for conventional lifting and projecting operators typically required in previous neural operators. Moreover, it seamlessly accommodates diverse boundary conditions. Our empirical obser- vations reveal that MgNO exhibits superior ease of training compared to CNN- based models, while also displaying a reduced susceptibility to overfitting when contrasted with spectral-type neural operators. We demonstrate the efficiency and accuracy of our method with consistently state-of-the-art performance on different types of partial differential equations (PDEs).
INTRODUCTION
Partial differential equation (PDE) models ubiquitous in physics, engineering, and other dis- ciplines. Many scientific and engineering fields rely on solving PDEs, such as optimizing airfoil shapes for better airflow, predicting weather patterns by simulating the atmosphere, and testing the strength of structures in civil engineering. Tremendous efforts have been made to solve various PDEs arising from different areas. Solving PDEs usually requires designing numerical methods that are tailored to the specific problem and depend on the insights of the model. However, deep learning models have shown great promise in solving PDEs in more general and efficient way. Recently, several novel methods, including Fourier neural operator (FNO) (Li et al. (020). Galerkin trans- former (GT) (Cao) 202T), deep operator network (DeepONet) (Lu et al. 2021), and convolutional neural operators Raonic et al. (2023), have been developed to directly learn the operator (mapping) between dimensoaail parameter and solution spaces of PDEs. These methods leverage ad- vanced architectures, such as Fourier convolution and self-attentions, to handle variety of input parameters and achieve high performance for forward and inverse PDE problems. The methodology has been applied to biomechanical engineering and weather forecasting (You et al. [022] Pathak et al. (022). Although deep neural operators have shown great potential to learn complex patterns and relationships from data, they still have some limitations compared to the classical numerical methods, such as lower accuracy and less flexibility to complex domain and boundary conditions. To address these challenges, we investigate the intrinsic properties of PDE-governed tasks and de- sign the neural operator architecture based on multigrid methods, one of the most commonly used and efficient classical numerical algorithms.
'These authors contributed equally to this work. "ooresponniig author.

---

Published as a conference paper at ICLR 2024
In this work, propose concise and elegant neural network architecture for operator learning. Our contributions are summarized as follows:
We introduce novel formulation for neural operators where the connections between neu- rons are characterized as bounded linear operators within function spaces. This approach yields new universal approximation result, making the commonly used lifting and pro- jecting operators in traditional constructs unnccessary.
Central this formulation is the efficient parameterization of linear operators. In response, we propose a distinet neural operator architecture, denoted MgNO. It employs the multi- grid, with standard convolutions to parameterize the linear operator between neurons, align- ing with our overarching goal of conciseness. Given the inherent ties between convolutions and multigrid methods in PDE contexts, MgNO naturally accommodates various boundary conditions.
Our method establishes its superiority not only in prediction accuracy but also in efficiency, both in terms of parameter count and runtime several PDEs, including Darcy, Helmholtz, and Navier-Stokes equations, with different boundary conditions. Such efficiency and pre- cision resonate with the core philosophy of the multigrid method. Despite its minimal parameter count, our model exhibits rapid convergence without susceptibility to overfitting compared with other models, indicating a problem worth further investigation from both theoretical and practical perspectives.
2 BACKGROUND AND RELATED WORK
Neural operators recent times, substantial efforts have been dedicated to the development of neural networks or operators tailored for solving partial differential equations (PDEs). Lu et al. (202T) introduced the DeepONet, employing a branch-trunk architecture grounded in the universal approximation theorem in Chen et al. (1996). FNO parameterizes the global convolutional opera- tors using fast Fourier transform. Furthermore, geo-FNO et al. (2022) addresses tasks involving intricate geometries, such as point clouds, by transforming data into latent uniform mesh and back. U-NO Rahman et al. (2022) improves FNO with a U-shaped architecture. The latent spec- tral model (LSM) (Wu al. (003) tried to identify the latent space and represent the operator in latent space. In addition, MWT Gupta et al. (202T) introduces a multiwavelet-baste
operator by parametrizing the integral operators using a fast wavelet transform. Then, Cao explored the self- attention mechanism Vaswani et al. (2017) and introduced a Galerkin-type attention mechanism with linear complexity for solving PDEs. More recently, a convolution-basee
architecture has been employed as a neural operator, as presented in Raonic et al.] (2023). Among the notable approaches, spectral-type neural operators have gained attention for their global operation characteristics. The non-local operators have shown very promising results. Mathematical viewpoints provide universal approximation property for such operators (Kovachki et al. [202]] Lanthaler et al. 2023). However, their inherent limitation lies in their incapability effectively handle boundary conditions. Meth- ods like FNO, MWT, UNO, and LSM, rely on encoding boundary information within the input data, making them reliant specific data representations. These methods tend to learn www-fequency component (Liu et al. 2022), which is not surprise since high-frequency modes are often truncated for the spectral type method. In these methods, convolutions in the neural operator construction parametrized by Fourier, or wavelet transforms. However, Fourier or wavelet- based methods are not always appropriate for solving PDEs because of boundary conditions and aliasing errors. In con- trast, local and geometrical neural networks, such as ResNet, DitResNet, UNet He et al. (2016a b); Ronneberger et al.] (2015), offer greater flexibility in managing diverse boundary conditions. Nev- ertheless, these approaches were not originally designed for solving PDEs and often necessitated substantial number of parameters to achieve high accuracy. Furthermore, they tend to act as high- pass filters, focusing primarily on local features like edges, surfaces, and textures, which limits their ability to capture long-diss -isttane relationships crucial various physical phenomena. Large kernel convolution could mitigate the issue, but the design of large kernels requires hand-craftes since it is often hard to train (Ding et al. 2022).
Multigrid Multigrid methods Hackbusch 2013) Xu 1989), Trottenberg al. (2000) stand out as
some of the most efficient numerical approaches in scientific computing, especially when tackling elliptic PDEs. Intriguingly, within the deep learning community, multigrid methods first received
2

---

Published as a conference paper at ICLR 2024
mention in the original ResNet paper He al. (016a), where the authors cited the methods core rationale behind the utility of residuals. Subsequently, the authors He & Xu 2019); He et al. (2023) established deeper and more extensive links between multigrid methods and ResNet, giving rise to what is now known as MgNet. The research rigorously demonstrated that the linear V-cycle multigrid structure for the Poisson equation, can be represented as convolutional neural network, despite their experimental focus on nonlinear multigrid with single-cycle structure for vision-related tasks. The original MgNet and its subsequent adaptations in the context of numerical PDEs Chen et al. (2022) and forecasting scenarios (Zhu et al. (2023), have not provided a simple but effective solution in the realm of operator learning. In this work, we aim to fundamentally integrate multigrid methodologies with operator learning based on very general framework.
3 A MORE GENERAL DEFINITION OF NEURAL OPERATORS
In this section, we present general deep artificial (abstract) fully connected (regarding neurons) neural operators maps between two Banach function spaces X H²( and y H* (Ω2) bounded domain Ω c Rdd
3.1 SHALLOW NEURAL OPERATORS FROM X TO y
To begin, define the following shallow neural operators with n neurons for operators from X to y as
O(u) Σ A (WW, B))
VueX
(I)
where W, E L,X.M), E V, and A, E LVVJ)). Here, L(X.)) denotes the set of all bounded (continuous) linear operators between and y, and or HR defines the nonlinear point-wise acti- vation.
Despite the absence of commonly used lifting and projection layers we still have the following universal approximation theorem based this unified definition of shallow networks.
Theorem 3.1 for any continuous operator O² X + y, compact set C c X and e > O, there is such that
H22)) and H² (12) for some s, s' 2 1, and o e C(R) is non-polyooial,
inf sup (ppll("()) O(u)lly Se.
OE uec
where E, denote the shallow networks defined in equation) with neurons.
Sketch of the proof: We begin by approximating (u) Σ f( where this is viewed as a piecewise-constant operator. For this approximation, o; E y and fi X - R are continuous functionals for all i m. Subsequently, we approximate fi(u)(); by the aforementioned shallow neural operator (₁()) by using the properties of X = H²) and the classical universal approximation results for shallow neuron networks on R detailed proof of this can be found in Appendix A
3.2 DEEP NEURAL OPERATORS FROM X To y
For clarity, let us define the product space y" := y a......) ay. Subsequently, the bounded
linear operator acting between these product spaces can be expressed W (VV" y"") To be more specific, the relation [Wh)], = Σド= Wi,h, holds for i = == 1: m, where h, E y and Wiy E (V. Now, let us denote E N+ for all l = 1 L as the number of neurons in (-th hidden layer with ni,1 1. Then, the deep neural operator with L. hidden layers and nN neurons in {-th layer is defined as
ho uEX h' h'(w) ( B¹) y (=1:L O(u) wiih y
(3)
where W ( , with W E L(X. and B' E y Unless otherwise specified, for simplicity, in the rest of the article, we assume n1

---

Published as a conference paper at ICLR 2024
We note that a standard hidden layer in GNO Kovachki et al. (2023 FNO Li al. (2020) given
by
h tu()) =o (wwn (w)+ B' +b) EY" y".
(4) B E R and b' E R.. The term Bhh
where (w"h (x) x,xx,x) in GNO FNO is typically referred to as local linear operator. In approach, the term is interpreted as a specific parameterization of 3' in equation 3 To parameterize B' E Y", the task
100) boils down to selecting appropriate bases for y. If we adopt the bases ( .... ,100) y
for each
y with i = I,....II, then we can express / 2 B B + b/1((x) Further, this can be compactly written vector form as B B'h -1 + b³, where B IR and b' E RR.
[x]
1
From equation a notable feature of our framework is the elimination of the commonly used lifting layer in the first layer and the projection layer the last layer, which are prevalent most prior neural operators Kovachki et al. (2023); Li et al. (2022). It's important to highlight that these two artificial layers are considered crucial for demonstrating the approximation capabilities across a wide range of neural operators, evidenced in Lanthaler et al. (2023).
However, in our proof of Theorem 3.1 a key insight that allows us to remove the lifting and projec- tion layers is the high expressive capability of the linear operator ", e (V, which maps from the i-th input neuron to the j-th output neuron for all i, 1 As result, primary contribution of this work is the introduction of a convolutional multigrid structure to parameterize W
4 PARAMETRIZATION OF W. AND MGNO
In this section, we begin by discussing the rationale behind using multigrid to parameterize W We then describe a standard multigrid process, framed in convolution language (single channel), for solving an elliptic PDE. Subsequently. we outline the global parameterization of W inspired by the multi-channel multigrid structure. Finally, we introduce our neural operator, MgNO, and touch upon its Houndary-preserving characteristics.
Motivation Given the Schwartz kernel theorem Schwartz (1950 and 1951); Duistermaat et al. 2010)) it is reasonable to select Wi as kernel version, i.e. W. "Wi,a(x) = K, x' )u(x' )dx' as equation
[ The primary challenge lies in appropriately parameterizing kernel functions K,, on Ω x Ω. We further narrow our focus to a specific type of kernel known as Green's functions, which correspond certain elliptic PDEs with specific boundary conditions. For practical purposes, we consider the discrete case for direct operator learning. Specifically, we have Q = [0, 11 for 1,2,3 and
(5) For instance, when d 2, the input is discretized u(x) = Ei uijhe(()) to Rdxd, where d A. Consequently, for any W VV, the dimensionality of the space L (V) becomes d² x d². This dimensionality is challenging to parameterize directly for large values of d. Most existing neural operators, such as those in Kovachki et al. [2023); Lanthaler et al. (203), primarily offer low-rank approximations of the kernel function in the spectral or wavelet domain. In this work, we advocate for using multigrid structures to directly parameterize within the spatial domain.
X y (V)) := linear finite element space on Q with mesh size h.
4.1 MULTIGRID METHODS FOR DISCRETE ELLIPTIC PDEs WITH BOUNDARY CONDITIONS
First, we offer concise and practical overview of the rationale and methodology behind using multi- grid techniques to parameterize Green's 's functions. Consider the elliptic PDEs given by Lu(x) = f() defined over the domain Ω (0. 1)2 and subject to Dirichlet, Neumann, or periodic boundary con- ditions. Employing linear FEM discretization with mesh size defined as h J. the discretized system can be expressed as:
A A+u=f.
(6)
where u.f f Rddd Here, represents the standard convolutional operation for a single channel, complemented by specific padding schemes determined by the boundary conditions. The kernl
4

---

Published as a conference paper at ICLR 2024
A, of dimensions X 3, is dictated by the elliptic operator L in conjunction with the linear FEM. Consequently, the inverse operation A* corresponds to the discrete Green's function associated with L' under lincar FEM framework. As demonstrated in [e & Xu (2019): He et al. (202T), the V-cycle multigrid approach for solving equation can be precisely represented as a conventional convolutional neural network with one-channel.
We provide a concise overview the essential components of multigrid structure the language of convolution an operator mapping from
1. Input (/) and Initialization:
Set f f and initialize with u = O.
2. Iteration (Smoothing) Process: The algorithm iteratively refines u based on the relation:
ut nd Bddtt (f' -A A 1fx
where = 1 and i ve.
3. Hierarchical Structure via Restriction and Prolongation: The superscript denotes the specific hierarchical level or grid within the multigrid structure. Specifically, using the residual, restrict the input f and the current state u' to coarser level through convolu- tion with a stride of 2:
= = RCt *2 (f' + u')
" E+10 O.
Subsequently, we apply the smoothing iteration as in equation 7 to derive the correction from the coarser level. The correction is then prolonged from the coarser the finer level using a de-convolution operation P/,, with stride of (acting as the transpose of the re- striction operation). After this prolongation, one either opt prolong further an even coarser level (known as the Backslash-cycle) or proceed with post-smoothing followed by prolongation (referred as the V-cycle).
Let's represent the linear operator defined by the aforementioned V-cycle multigrid operator as WMg. The convergence result presented in Xu Zikatanov (2002)
! ||u|A
||u WM (A
4) 40llA
demonstrates the uniform approximation capabilities of WMy relative to the inverse of A+, which corresponds to the Green's 's function associated with the elliptic operator L. Here, is constant that independent of the mesh size h, and ||ы][ = (u, A 4) denotes the energy norm. Please refer the Section C in Appendix for more details regarding the approximation property. We quantify the constant e equation numerically using a concrete example.
Invu if. wftoo
nestion
Imput:
resticc
wiltSleyer
Figure 1: Overview of WMs using multi-channel V-cycle multigrid framework.
4.2 ARCHITECTURE OF MoNO
Finally, we introduce a surrogate operator, denoted as MgNO. This operator maps from the linear finite element space of input functions, represented by Ruxdxe. X VV(2), to the linear
5

---

Published as a conference paper at ICLR 2024
finite element space of output functions, represented by Raxdxea =y (2( The mapping defined as:
= =EX. o h-) (w) +b10)+11eee EM, yr e=l:L.
(10)
Mş
0= Golu)= = Mg (W') y.
Here, W MM V"X, W Eyy" for t= J and W LV ,) are multi-channel
Mg
linear operators that process multi-channel input h- and output W (h(-). The function or is the point-wise GELU activation, with B' R=xa and b' R..
Mg
We observe that directly using a one-channel "Mr to parameterize each W in equation 3 indi- vidually lead to expressivity and efficiency limitations. It is advantageous restriction operators R and prolongation operators Pe,, for l = 1 J. as multi-channel convolu- tional kernels. detailed exposition of the multi-channel multigrid in the convolutional context is available in Algorithm|P Appendix and figure I Two primary motivations drive our choice of this strategy. Firstly, as established in He et al. (2022). channels in deep CNNs, for 2D inputs, function analogously to neurons in terms of their universal approximation capabilities. Consequently, aug- menting the number of channels multigrid enhances the expressiveness approach aligns with classical multigrid methods applied to solve elliptic partial differential equa- tion system, including linear elasticity problems in 2D or 3D scenarios. In such cases, the V-cycle multigrid methods can be naturally represented linear MgNO with multi-channel inputs. From practical standpoint, the multi-channe variant of WM possesses a parameter count on the order of O (log(d)n²) and complexity the order of Od?rr²). Ultimately, by substituting W in equation with WW Mg we derive our MgNO as presented in equation 10 It's important to note that W does Mg
design all A B
Me- Secondly, this
not include any nonlinear activation.
Goundary-preserving discretization in MgNO It is imperative to highlight that MgNO is adept at accommodating the boundary conditions of various PDEs. Leveraging the inherent relationship between convolutions and multigrid, the output of W Mg can seamlessly satisfy Dirichlet, Neumann,
or periodic boundary conditions. Consequently, Go(u) can maintain these boundary conditions with- out any training phase, provided by taking B' = O in equation 10 and ensuring 1(x) satisfies the same boundary conditions. One can refer Table Section for specific convolution configura- tion for specific boundary conditions, and also Section C in Appendix for an concrete example of integrating boundary conditions into neural operators.
EXPERIMENTS
General setting We consider pairs of functions (uj, D where uy is drawn from probability measure μ and ;j= G(u,). Given the data (u,, vi)a we approximate G, by solving the network parameter set 0 via an optimization problem:
min L(0) == min AED deo
N [IG.tu)) - o,
(1I)
Benchmarks To substantiate the superiority of our method, we conducted comprehensive method comparisons across multiple benchmark scenarios. These benchmarks were categorized based on their association with specific partial differential equations (PDEs, namely Darcy, Navier-Stokes, and Helmholtz. benchmarks. Furthermore, we incorporated tasks encompassing both regular and ir- regular domains, each subject to variety of boundary conditions. This comprehensive comparison allows us to thoroughly assess the robustness and adaptability of our method, ensuring its effective- ness in real-world scenarios. For clearness, we summarize the benchmarks in Table S Please also refer to Section D for more detailed descriptions.
We perform a comprehensive comparison with existing methods: (1) We include the original FNO (Li et al. 2020); (2) UNet (Ronneberger al. 2015) and U-NO (Rahman et al.
Baselines

---

Published as a conference paper at ICLR 2024
2022 (3) multiwavelet neural operator (MWT) (Gupta et al. 2021 (4) Galerkin transformer (GT)
Cao 2021); (5) latent spectral models (LSM) (Wu et al. 2023); (6) learned simulator for turbu-
lence (DilResNet) (Stachenfeld et al. 202T). (7) Fine-tuned ENO (sFNO-v2) (Benitez et al.] 2023).
The baseline models were implemented using their official implementations. For consistency, we executed the baselines with default hyperparameters unless the experiments' details specified oth- erwise. In cases where specific experiment details were lacking in the literature, we employed a reasonable network scale to ensure fair comparison. In our training setup, we implemented com- prehensive approach. We trained the baselines using multiple configurations, including their default training settings encompassing different loss functions, training algorithms, and schedulers. We then reported the best results from these trials to mitigate any training-related variations. Notably, our observations revealed that all models benefited from the scheduler utilizing cosine annealing learn- ing rates. Figure 3 will visually demonstrate that all models underwent sufficient training, ensuring that their performance accurately reflects their capabilities.
Darey The Darcy problems, widely recognized for their applicability to multigrid methods, serve as a key focus. We anticipate substantial performance improvements compared to the current state- of-the-art. To thoroughly assess our method, we incorporated three distinct Darey benchmarks: Darcy smooth, Darey rough in Li et al. (020)), and Darcy multiscale in [iu et al. (022).. Darcy smooth and Darcy rough both present two-phase media within the domain, featuring different in- terface roughness levels. In contrast, Darcy's multiscale presents media with multiple scales of permeability. Darcy rough and multiscale Darcy challenges our method's s ability to capture operator dependencies fine-scale features, task that poses difficulties for existing methods. Please refer D.I for detailed datasets description and experiments setup.
Table I: Performance comparison for Darcy benchmarks. Performance are measured with relative La errors (x10-3) and relative H' errors (x10-
Darcy smooth
Darey rough
Darcy multiscale
Model
time (s/iter)
params (m)
H
H
FNO2D
7.4
2.37
0.684
2.583
1.613
7,516
1,800
9.619
Du RESNET
14.9
1.04
4.104
5.815
7.347
12.44
1.417
3.528
UNET
9.1
17.27
2,169
4.885
3,519
5.795
1.425
5.012
U-NO
11.4
16.39
0.492
1.276
3.784
1.187
MWT
21.7
9.80
1.138
4.107
1,021
7,245
GT
38.2
2.22
0,945
3,365
1.790
6,269
1,052
8.207
LSM
18.2
4.81
0.601
2,610
2.658
4,446
1,050
4,226
MaNO
6.6
0.57
0.176
0.576
0.339
1,380
0.715
1.756
McNO-Baat-
0.103
0.275
0.799
0.504
.WT (Gupta et al. 2021 only supports resolution with powers of two.
-NN2D. U-NO'and and MWT's 's performance are further improved from originally reported because of the usage of H' loss and schedulor. The runtime number of parameters count are using Darcy rough as the example case.
Remark 5.1 Note that MGNO-HIGH- is not directly comparable with other models. The McNO- HIGH-IN model differs from others, utilizing a higher input resolution of x 512 the Darey rough case, compared to the standard 256 256 used in other models. Inspired by reduced-order modeling (Lucia et al.) 2004 Hou et al.) 1999- Engquist & Souganidis) 2008 Mālqvist & Peterseim 2204 Owhadi & Zhang 2007), it represents high-resolution input but targets a solution in lower dimen- sional space. This is achieved by adjusting the depthflevel of the W that condenses the input to the desired output resolution. Our findings suggest that high-resolution inputs enhance fine-scale predictions, implying the model effectively balances computational efficiency with high fidelity On the other hand, by incorporating appropriate finite element basis functions, MgNO, compatible with the finite element method, demonstrates the discretization invariant capability. MgNO can train on rower-resolution datasets and evaluate on higher resolutions without the necessity of risolutstttto training data, thus accomplishing zero-shot super-resolution. Please see E for the detailed experi- ments.
Table 5 and Figure 2 illustrate the precision with which our method predicts the fine-scale features of the solution. The qualitative comparisons of the contour lines in the top row indicate the model's aptitude in accurately capturing scall-scale variations, evidenced further by the low H' error re- ported in Table 5 Moreover, our method establishes its superiority not only in prediction accuracy
7

---

Published as a conference paper at ICLR 2024
Figure 2: Qualitative comparisons on Darcy rough benchmark. Top: coefficient a. ground truth u, and predictions; bottom: the corresponding prediction error map for each model in the color scale.
but also in efficiency, both in terms of parameter count and runtime. Such efficiency and precision resonate with the core philosophy of the multigrid method.
0.0
0.43
0.0
MgNOOO..51 57M perameters)
FNO12.37M parameters)
UNDI20 39M parameterss)
Figure 3: Comparison of training dynamics between MgNO, FNO and UNO. The x-axis represents the number of epochs, and the y-axis is the error in the log scale. We present both the L and H training and testing accuracy (errors). For full comparisons, please refer D.I
Training Dynamies Figure 3 we present the training dynamics different models on bench- mark Darey rough. We train all the models in the same setting. Please see D.I for further training set- tings). Despite its minimal parameter count, our model exhibits rapid convergence without suscep- tibility to overfitting. It is noteworthy that while larger models like MWT(. 80MM, UNO(16.39M), and LSM(4.81MM can achieve training errors comparable to, or even lower than, our model, their testing are considerably higher. This discrepancy underscores model's robust generaliza- tion capabilities on unseen data.
Helmholtz We conduct experiments the Helmholtz equation. The Helmholtz equation poses significant challenges for numerical methods, primarily due to the resonance phenomenon. Please refer D.3 for more details about the dataset. We present the comparison in Table and Figure
Navier Stokes Numerous studies have explored the use of neural operators for fluid dynamics simulation Kochkov et al. (2021): :: Li et al. (0020); Stachenfeld et al. (202I): Mi et al. (2023). Our focus is on the 2D Navier-Stokes equation its vorticity form defined over the unit torus, T (refer to Section D.2 for an in-depth discussion). The vorticity, represented by w(xx, f). is defined for x T), which projects
and T). We aim to learn the operator S w(. 0 < s 9) w(: 10 the vorticity up to time 9 onto the vorticity up to a later time T. following the benchmark set by Li al. 2020). Our experiments consider viscosities of le 5 and conclude at time 7 20.
8

---

Published as a conference paper at ICLR 2024
Furthermore, also test models on the same Navier-Stoker task with -e - 5 but introduce an alternative training configuration labeled as "NS2". It is imperative to note that in the Navier Stokes le- 5 configuration with training tricks, all models incorporate these specialized training techniques. These modifications consistently enhance performance when compared to the original setup. The inclusion of these training tricks contributes more stable generalization errors across various models, justifying their use in performance comparisons. Hence, both evaluation methods offer an equitable basis for contrasting the efficacy of our approach with existing baselines. A com- prehensive overview of the training configurations and specific tricks can be found in Section D.2
Table 2: Performance on Navier Stokes
Table 3: Performance Helmholtz
Model
time
paramiima
Pipe
NS
NS2
Model
time
params(i))
Pox10 3)
ff'000 3)
FNOID
2.0
0.94
0.67
15.56
4.57
ENO2D
5.1
1.33
1.68
9.92
UNET
7.6
17.27
0.65
19,82
10.08
UNe
7.1
17.26
3.81
23.31
DiRReNNer
5.9
1.04
0.31
23.28
10.91
DERISNET
4.34
34.21
U-DD
21.5
1.26
MWT
6.2
9.80
0.77
15,4]
3.09
sENO-v2
30.1
12.0
1.72
10.40
LSM
24.3
4.81
0.50
15.35
12.93
LSM
28.2
4.81
2.55
10.61
MONO
6.4
2,89
0.24
10.82
1.63
M&NO
7.58
0.71
4.02
Ablation & Hyperparameter Study The
ablation study uses NS2 as reference exam- plc. The baseline MgNO is configured with 5 levels, 32 channel dimensions, utilizes periodic boundary conditions, and comprises 5 architec- tural layers. One can find the detailed config- uration in Table 7 We conduct the ablation& hyperparameters study by changing one model configuration and fixing the others. In the first and second rows, our analysis evaluated the im- pact of the number of levels (depth) within the MgNO architecture. Given that structure simply merges convolution with a multigrid touch, this hyperparameter examination doubles as an ablation study for the influence of neural network depth. The outcomes highlight the performance sensitiv- ity of MgNO to its depth. This mirrors properties observed in traditional multigrid methods where deeper algorithms are requisite for effectively mitigating low-frequency component errors. The third row emphasizes the role of boundary conditions in convolution operations, underscoring the distinct advantages of our method. Fourth row indicate adding more layers improves performance, but also increases computational costs. The MGNO-FNO configuration offers higher-level ablation. Here, we substituted the multigrid parameterized Wi, with the spectral convolution detailed in Li et al. (020). The ablation study reaffirms the strength of the McNO approach, especially when integrat- ing boundary conditions. This method not only resonates with traditional multigrid propertics but also encapsulates the essence of convolution- based neural operators, delivering optimal performance in comparison to alternative convolution- -based configurations.
Table 4: Ablation & hyperparameter study
Model Configuration
L Emur
McNO. levels
2.10
McNO, levels
2.44
McNO, without houndary condition
3.18
McNO, layers
1 1.47
MGNO-PNO
3.91
Baseline MgNO
L.63
CONCLUSION
In this work, we introduced novel formulation for neural operators, where neuron connections are characterized as bounded linear operators within function spaces, eliminating the need for tra- ditional lifting and projecting operators. Central our approach the MgNO architecture, which leverages multigrid structure and its multi-channel convolutional form to efficiently parameterize linear operators, naturally accommodating various boundary conditions. Empirically, MgNO has demonstrated superior performance in both accuracy and efficiency across several PDEs, including Darcy, Helmholtz, and Navier-Stokes equations.
Moreover, there are still some questions worth future exploration. Our work does not sufficiently benchmark tasks on irregular domains. For irregular domains, our current implementation of MgNO requires deformation mapping to transform the data into regular format as shown in the pipe task. Therefore, incorporating algebraic multigrid methods Xu & Zikatanov (2017) will enhance the MgNO's applicability to inputs and outputs in general format such as arbitrary sampling data points and point clouds. Furthermore, MgNO achieves rapid convergence without overfitting, highlighting promising direction for future research.
9

---

Published as a conference paper at ICLR 2024
Reproducibility Statement The complete code used in experiments is available at https::/ github. www./ooloooooii..oce /MgNO/ This repository includes all scripts, functions, and necessary files for reproducing our results. Datasets employed in experiments provided in the repository. We have disclosed all configurations, encompassing hyperparameters, model architectures, and optimization strategies. This ensures that experimental setup is fully transparent and can be replicated by others. To ensure the robustness of our results, each experiment was conducted at least three times with different runs initiated using unique random seeds. This approach is aimed at accounting for variability and guaranteeing the reliability of our findings.
also be accessed via URLs
Acknowledgement
This work was partially supported by the KAUST Baseline Research Fund.
REFERENCES
Jose Antonio Lara Benitez, Takashi Furuya, Florian Faucher, Xavier Tricoche, and Maarten de Hoop. Fine-tuning neural-operator architectures for training and generalization. arXiv preprint arXiv:2301. 11509, 2023.
Dietrich Braess. Finite elements: Theory, fast solvers, and applications in solid mechanics. Cam- bridge University Press, 2007.
Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural pde solvers. arXiv preprint arXiv:2202.03376 2022.
William L Briggs, Van Emden Henson, and Steve F McCormick. A multigrid tutorial. SIAM, 2000.
Shuhao Cao. Choose transformer: Fourier galerkin. Advances in Neural Information Processing Systems, 34, 2021.
Jiun-Shyan Chen, Chunhui Pan, Cheng-Tang Wu, and Wing Kam Liu. Reproducing kernel particle methods for large deformation analysis of non-linear structures. Computer Methods in Applied Mechanics and Engineering, 139(1-4):195-227, 1996.
Yuyan Chen, Bin Dong. and Jinchao Xu. Meta-mgnet: Meta multigrid networks for solving param- eterized partial differential equations. Journal of computational physics, 455: 10996, 2022.
Maarten De Hoop, Daniel Zhengyu Huang, Elizabeth Qian, and Andrew Stuart. The cost- accuracy trade-off in operator learning with neural networks. arXiv preprint 13X11. 2022.
Xiaohan Ding, Xiangyu Zhang, Jungong Han, and Guiguang Ding. Scaling up your kernels to 31x31: Revisiting large kernel design cnns. In Proceedings of the IEEEJCVF conference on computer vision and pattern recognition, pp. 11933-11975, 2022.
Johannes Jisse Duistermaat, Johan AC Kolk, JJ Duistermaat, and JAC Kolk. Distribution Kernels. Springer, 2010.
B. Engquist and P. Souganidis. Asymptotic and numerical homogenization. Acta Numerica, 17: 147-190, 2008.
Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavelet-basee operator learning for differen- tial equations. Advances in Neural Information Processing Systems, 34:24048-2444 24062, 2021.
Wolfgang Hackbusch. Multi-grid methods and applicatiinss, volume 4. Springer Science & Business Media, 2013.
Juncai He and Jinchao Xu. Mgnet: A unified framework of multigrid and convolutional neural network. Science China Mathematics, pp. 1-24, 2019.
Juncai He, Jinchao Xu, Lian Zhang, and Jianqing Zhu. An interpretive constrained linear model for resnet and mgnet. arXiv preprint arXiv:2112.07441, 2021.
Juncai He, Lin Li, and Jinchao Xu. Approximation properties of deep relu cnns. Research in the mathematical sciences, 9(3)338, 2022.
10

---

Published as a conference paper at ICLR 2024
Juncai He, Jinchao Xu, Lian Zhang, and Jianqing Zhu, An interpretive constrained linear model for resnet and mgnet. Neural Networks, 162:384-399 -928, 2023.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Proceedings of the IEEE Conference Computer Vision and Pattern Recognition, pp- 770-778, 2016a.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Identity mappings in deep residual networks. In European Conference Computer Vision, pp. 630-645. Springer, 2016b.
T. Y. Hou, X.-H. Wu, and Z. Cai. Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients. Math. Comp., 68(227): 13--993,, 1999.
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine learning-accelenee accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences, :11021):e2101184111 2021.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces. arXiv preprint ar arv: 2tt.::220.0000. 088808888888.8888811111881111181 2021.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew M Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. J. Mach. Learn. Res., 24(89):1-97, 2023.
Samuel Lanthaler, Zongyi Li, and Andrew M Stuart. The nonlocal neural operator: Universal ap- proximation. arXiv preprint arXiv:2304. 13221, 2023.
Moshe Leshno, Vladimir Ya Lin, Allan Pinkus, and Shimon Schocken. Multilayer feedforward net- works with nonpolynomial activation function can approximate any function. Neural networks, 6(6):811-867 1993.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier NeuralOperator for Parametric Partial Differential Equa- tions. In International Conference on Learning Representations, 2020.
Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for pdes on general geometries. arXiv preprint arXiv:2207.05209, 2022.
Xinliang Liu, Bo Xu, and Lei Zhang. Ht-net: Hierarchical transformer based operator learning model for multiscale pdes. arXiv preprint arXiv:2210. 10890, 2022.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3(3):218-229, 2021.
David Lucia, Philip S Beran, and Walter Silva. Reduced-order modeling: new approaches for computational physics. Progress in aerospace sciences, 40(1-2):51-117, 2004.
A. Mâlqvist and D. Peterseim, Localization of elliptic multiscale problems. Math. Comp., 83(290): 2583-2603, 2014.
Junyi Mi, Xiaowei Jin, and Hui Li. Cascade-net net for predicting cylinder wake at reynolds numbers ranging from suberitical to supercritical regime. Physics of Fluids, 35(7), 2023.
H. Owhadi. Multigrid with Rough Coefficients and Multiresolution Operator Decomposition from Hierarchical Information Games. SIAM Rev., 59(1):99 149999-19,,90 2017.
H. Owhadi and L. Zhang. Metric-based upscaling. Comm. Pure Appl. Math., 60(5):675-723, 2007.
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay. Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al. Four- castnet: A global data- driven high-resolution weather model using adaptive fourier neural opera- tors. arXiv preprint arXiv:2202. 12 2022.

---

Published as a conference paper at ICLR 2024
Md Ashiqur Rahman, Zachary E Ross, and Kamyar Azizzadeneshelii U-no: U-shaped neural oper- ators. arXiv e-prints, pp. arXiv-2204, 2022.
Bogdan Raonic, Roberto Molinaro, Tobias Rohner, Siddhartha Mishra, and Emmanuel de Bezenac. Convolutional neural operators. In ICLR 2023 Workshop Physics for Machine Learning, 2023.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedi- cal image segmentation. In International Conference on Medical image computing and computer- assisted intervention, pp. 234 241. Springer, 2015.
Laurent Schwartz. Théorie des distributions. and II. Paris, 1950 and 1951.
Kim Stachenfeld, Drummond Buschman Fielding, Dmitrii Kochkov, Miles Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho, Peter Battaglia, and Alvaro Sanzhez- Learned simulators for turbulence. In International conference on learning representations, 2021.
Ulrich Trottenberg, Comelius W Oosterlee, and Anton Schuller. Multigrid. Elsevier, 2000.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention all you need. Advances in neural informa- fion processing systems, 30, 2017.
Haixu Wu, Tengge Hu, Huakun Luo, Jianmin Wang, and Mingsheng Long. Solving high- dimensional pdes with latent spectral models. arXiv preprint arXiv:2301. 2023.
Jinchao Xu. Theory of multilevel methods, volume 8924558. Cornell University, 1989.
Jinchao Xu and Ludmil Zikatanov. The method of alternating projections and the method of sub- space corrections hilbert space. Journal of the American Mathematical Sociery 15(3)55735977, 2002.
Jinchao Xu and Ludmil Zikatanov. Algebraic multigrid methods. Acta Numerica, 26:591-721, 2017.
Huaiqian You, Quinn Zhang, Colton J Ross, Chung-Hao Lee, Ming-Chen Hsu, and Yue Yu. A physics-guided neural operator learning approach to model biological tissues from digital image correlation measurements. Journal of Biomechanical Engineering, accepted for publication, pp- arXiv preprint arXiw:2204.00205, 2022.
Jianqing Zhu, Juncai He, and Qiumei Huang. An enhanced v-cycle mgnet model for operator learn- ing in numerical partial differential equations. arXiv preprint arXiv:2302.0093: 00938, 2023.
12

---

Published as a conference paper at ICLR 2024
A PROOF OF THEOREM 3.1
Here we present the detailed proof of Theorem with the following steps.
1. Piecewise constant approximation of O": Since C c X compact set, for any > 0, there are (61. mml y and continuous functionals : R. for m such that
sup MEC
O'(u)
Mi 3
(12)
1m1
Thus, we only need to prove that there is O, 3n, such that
sup ( (zu)ly $ 3m uaC
(13)
2. Parameterization (approximation) of with finite dimensions discretize fi: Since X (I2( and C is compact, we can find k E N such that
sup
f( usC
Σ u. 1))
Vi = 1 m,
(14)
3m
Ly where fi are the orthogonal basis (2)) Then, for specific fi X + R, let us define the following fiittediieensionl continuous function F; R - R as
F2(x)
Σ 1=1
xxpi
Vce
(15)
where M := sup, sup (us s.).
3. Universal approximation of Fr on [-M, MJ using classical shallow neural networks Fx(x)
for x E [-M, Mj o R + R is not polynomial, given by the results in Leshno et al.
(199)), there is N+ and
F(()= Σ aioo(wiy x bi)).
(16)
where wy B and aij, bij R. such that
sup
[F(()) F.(x))
(17)
3mlldblly
4. Representation of (x) using Em, Now, let us define O, Em, as
m
O Σ An (Wip + Bij)
(18)
where
Wiju wi ((u,))), (u, 1())
(19)
and B,, b,,](x) E V, and
(v) = aj (v. 1)
|| 6
(20)
Here, 1(x) denotes the constant function Ω whose function value is 1. This leads to
O ((u,,)))) (u, (x)) d:
= (((, (u, (P))) 4;
(21)
= f Σ u. si) Φ di.
5. Triangle inequalities finalize the proof:
sup
O (w)
Σ
0,(0) SE
(22)
uec
l=1
where E Of(u) Ex with N Σ Bij
13

---

Published as a conference paper at ICLR 2024
B MULTIGRID IN CONVOLUTION LANGUAGE
Algorithm 1 WW WWMgSSuee u. J. Ye, vt,n)
1: Input: discretized input function f RRaddee, discretized initial output function u Rdixdixn number of grids J, number of smoothing iterations for = J, number of channels on each grid.
2: Initialization: f' K f Rdixdd u E Rdxidi u = 0 if n is not given as input, and
discretized spatial size d, z for 1=1:=
3: for (= 1 do
4: Feature extraction (smoothing):
5: for = 1 do
6:
ut w + Bld (f" A x ( Rdindon
(23)
7: end for
8: Note: uf u
9: if l then
Interpolation and restriction:
A² u uu)E Rde-ixdes u ++.0 and getl = -R(t) 42 *2((" -f---'
11: end if
12: end for
13: for -- do
14:
Prolongation:
d ul +
+P P/a x² ut)
15:
i- 1: Yt do Skip if Backslash-cycle OR do the following post post-smoothing if V-cycle: u u fi + a xf" - A (f'
ul u((1)) E Boxdon
17: end for
18:
end for
19: Output:
ut = lon E didden
C A NUMERICAL EXAMPLE ILLUSTRATING THE APPROXIMATION PROPERTY OF W Mg WITH SEPCIFIC BOUNDARY CONDITION IN EQUATION
We consider elliptic PDEs described by -A((x) f(x) over the domain = (0. 1)2, subject to Dirichlet boundary conditions where 0 on a.. Using linear Finite Element Method (FEM) discretization with mesh size (d+T) the discretized equation is represented as:
I
A A+u=f.
(24)
-1 0 where f E Rdsdd and = -1
0
-1
Here, " denotes the standard convolution operation
0
-1
for a single channel, incorporating a zero-padding scheme with padding size corresponding to the
Dirichlet boundary condition.
In the one-channel 'Wa, "WMg configuration, as detailed in Algorithm T we need further configurations to incorporating the Dirichlet boundary conditions. The convolutional layers A and B utilize zero- padding of size 1. In contrast, the layers RE(+l and pf Ptup for & 1 J, up without any
14

---

Published as a conference paper at ICLR 2024
10-1
10-3
Error
10-3
74
Number of layers
Figure 4: The error, quantified on logarithmic scale, numerically demonstrates the approximation rate, which is approximately 1 c 0.1, outlined in equation
padding. The convolution weights are configured as follows:
0
-1
0
1/64
A
-1
-1
B
1/64
12/64
1/64
0
-1
0
1/64
0
1/2
1/2
0
1/2
1/2
RCt1
1/2
1/2
and pt
1/2
1/2
1/2
1/2
1/2
1/2
0
Implementing "Wur My with these specific weights, apply it iteratively to the right-hand side f for / iterations, analogous to layers in the 'WM network. We define (e(l-1) := WM al)) with v := O. The convergence rate observed is illustrated in Figure
The configurations of these weights are based fundamental facts of piecewise linear finite ele- ments and multigrid algorithms. For further details on these configurations, one can refer to Briggs et al. 22000;; Braess (2007). For Neumann boundary conditions and periodic boundary conditions, the implementation of reflect padding and periodic padding, respectively, each with a padding size of 1, aligns seamlessly with the discretization of the finite element method, facilitating straightforward application without essential difficulty. The convergence behaviors are the same.
D FURTHER DETAILS ON BENCHMARKS
Benchmarks
Time dependent
Regular domain
# Dimension
BD
High resolution
Multi-channel input
SMOOTH
No
Dirichlet
No
DANCY RO0121
No
Yes
2D
Dirichlet
Yes
No
NAVIER STOKES
Yes
Yes
2D
Periedic
No
Yos
PIPE
No
Mixed
No
Yes
HELMHOUTZ
Yes
2D
Neumann
No
No
Table 5: Summary of benchmarks
D.I DARCY
The Darcy equation writes
(d(x)v((x)) = f(x) xE
(25)
u(x) =
x aD
where the coefficient Omin S a(x) S amax, Yx D, and the foreing term f HD)DR The coefficient to solution map is (D; (.D.,.) (D; R), such that S(a) is the target operator.
Two-Phase Coefficient (Darcy smooth and Darcy rough) The two-phase coefficients and solu- tions (referred Darcy smooth and Darey rough in Table generated according to htthtts///
15

---

Published as a conference paper at ICLR 2024
gi github.. cww. cww./oongyi-li/fourrr/oooo petuurrrrrrruraallle operaator/tree/nnastee//aaa geriiiatangeeatiinnnoooo and used as an operator learning benchmark in (Li et al. 2020; Gupta et al 2021 Cao 2021). The coef- ficients a(x) are generated according to := N (- with zero Neumann boundary conditions on the Laplacian. The mapping d/ R - R takes the value amax on the positive part of the real line and amin on the negative part. The push-forwart is defined in a pointwise manner. The forcing term is fixed as f(x) 1. Solutions are obtained by using a second-order finite difference scheme on suitable grid. The parameters Omax and dmin can control the contrast of the coeffi- cient. The parameter controls the roughness (oscillation) of the coefficient; a larger results in coefficient with rougher two-phase interfaces, as shown in Figure 2
Darey multiscale coefficient We examine equation 25 that features multiscale trigonometric coefficient adapted from (Owhadi, 2017). This corresponds to the Darey multiscale benchmark highlighted in Table formulated as
The domain, denoted D), spans the area 11 The coefficient a(x) is
[ 1+ cos(a;a(x) + x2))
sin(a - 3xnn)
a(x)
where each ag is uniformly distributed between 2 and 1.5 2 The forcing term is consistently set to f(x) = 1. Reference solutions are derived using Pi FEM on grid resolution of 1023 x 1023.
Refer Figure for visual representations of both the coefficient and the solution.
coefficient in log scale (b) reference solution 2D
Figure (a) multiscale trigonometric coefficient, (b) reference solution.
Training settings In the Darcy rough scenario, our data comprises 1280 training, 112 validation, and 112 testing samples. For the Darcy smooth and multiscale trigonometric cases, the split is 1000 training. 100 validation, and 100 testing samples. Training limited to 500 epochs for Darcy smooth and rough, and 300 epochs for Darcy multiscale.
We trained various baselines using multiple configurations, focusing primarily default training settings with varying loss functions and schedulers. The optimal results were selected to ensure consistency. The consistent boost in performance across all models observed using the H! loss and the OneCycleLR scheduler with cosine annealing. Specifically, optimal models are trained with batch size of 8, the Adam optimizer, and OneCycleLR with cosine annealing. Although learning rates varied to optimize training across models, MgNO started with rate of 5 x 10-4 decreasing to 2.5 10-6 FNO. UNO, and LSM default to 10 with weight decay models, like UNet and DilResNet, have a maximum rate of 5 x 10 As indicated in Figure thorough training process ensured equitable model comparison.
10 CNN-based
All experiments were executed
NVIDIA A100 GPU.
Training dynamics We present the training dynamics of different models on benchmark Darcy rough. We train all the models in the same setting (we use the same batch size 8, Adam optimizer and OneCycleLR scheduler with cosine annealing. OneCycleLR scheduler with a small learning rate at the end allows models to be trained sufficiently. Note that the learning rate varies for different models to allow models to be trained sufficiently) for 500 epochs and record the history of training and testing errors during the process.
16

---

Published as a conference paper at ICLR 2024
Irain ecc L2 Train 80C ML
train acc brain DcC
Tran DCC F2 - train acc
tedt ace 12 Ted NCC
tedt W1 1047
Text aCC In test aot
100
200
100
100
5
100
200
300
400
120
200
300
6OC
MgNOOO.557 57M perameters)
FNO12 37M parameters)
UNO016 39M parameteers)
Train
brann
Uan UOC
Scan acc H1 ce4t
buin train
bast 12
0.
0.0
0.0
0.03
100
100
100
300
400
1:00
200
LSM444 91M parameters)
UNet117.277 2.M parameters)
MWT(9.80M parameters!
Figure 6: Comparison of training dynamics between MgNO, FNO and UNO. The x-axis represents the number of epoch and the y-axis the error in log scale. We present both the 1.2 and H1 training and testing accuracy (errors).
D.2 NAVIER STOKES
We take the dataset to simulate incompressible and viscous flow on the unit torus, where the density of the fluid is unchangeable (p in Eq. In this situation, energy conservation is independent of mass and momentum conservation. Hence, the fluid dynamics can be deduced with:
U
Ow
Ōt + U Vw =WW
w/.00 Wo.
is the vorticity, Wo E the initial vorticity at O. In this dataset, viscosity is set 10-3 and the resolution of the 2D field is 64 64. Each generated sample contains 20 successive time steps, and the task is to predict the future 10-time steps based the past 10 time steps. The training dataset contains 1000 samples, while the testing dataset contains 100. In the following, we describe two training setups.
where U = (u, v) a velocity vector in 2D field, w x UI= -3 åx
Two training setups
Original training setup: Samples consist of 20 sequential time steps, with the goal of predicting the subsequent 10 time steps from the preceding 10. Using the roll-out prediction approach as described by Li et al. 2220)), the neural operator uses the initial 10-time steps to forecast the immediate next time step. This predicted time step is then merged with the prior 9 time steps to prediet the ensuing time step. This iterative process continues to forecast the remaining 10 time steps.
New training setup: Our findings indicate that an amalgamation of deep learning strate- gies is pivotal for the optimal performance of time-dependent tasks. While Li et al. (2020) employed the previous 10-time steps as inputs for the neural operator, our approach simpli- fies this. We find that leveraging just the current step's data, similar traditional numerical solvers, is sufficient. During training, we avoid model unrolling. Instead, from the initial 1000 samples consisting of 20 sequential time steps, we derive 19,000 samples with pairs
17

---

Published as a conference paper at ICLR 2024
of sequential time steps. These are then shuffled and used to train the neural operator to predict the subsequent time step based on the current one. For testing, we revert to the roll-out prediction method only the initial time step's ground truth available. These methods align with some techniques presented in Brandstetter et al. (2022). As shown in T the second approach provides better performance.
For both setups, with a cosine annealing strategy. The learning rate starts with x 10 and decays to I x 10 All models were trained for 500 epochs.
the 1.2 loss function, the Adam optimizer, and the OneCycleI.R scheduler
Pipe This dataset in Li et al. 2022 is devoted to simulate incompressible flow through a pipe.
The governing equations are:
U-U=0
aU
+ VU =f- 1 D-Vp+wV P
The dataset generated in the pipe-shaped structured mesh with the resolution of 129 x 129, For experiments, we adopt the mesh structure the input data, and the output is the horizontal fluid velocity within the pipe.
W
2000
111
30000
118
506
Figure 7: The pipe task: On the left, the ground truth of the pipe flow is displayed, whereas the right illustrates the prediction error of the MgNO.
D.3 HELMHOLTZ EQUATION
Helmholtz equation in highly heterogeneous media is an example of multiscale wave phenomena, whose solution is considerably expensive for complicated and large geological models. We adopt the setup from (De Hoop et al. 2022), for the Helmholtz equation on the domain D = [0, IJF. Given frequency 10 and wave speed field Ω R, the excitation field Ω - R solves the
equation
2
0 in Ω
(xxx
du 0 on 022, 0Q
ön
du
aQ2,
On
where 02; is the top side of the boundary, and 8Q other sides. The speed field is where
c(x) = 20 + tanh(č(x)), where ē is sampled from the Gaussian field ē (- and d are chosen to control the roughness. The Helmholtz equation is solved 100x 100 grid by finite element methods. We aim to learn the mapping from c 100x to u E R 100x 100 shown in Figure
as
18

---

Published as a conference paper at ICLR 2024
(a) Prediction
MANNO
(b) Error
Figure 8: The Helmholtz benchmark.
The Helmholtz equation is notoriously difficult to solve numerically. One reason is the so-called resonance phenomenon when the frequency w is close to an eigenfrequency ator for some particular wave speed We found that it necessary to large training dataset of size 4000 examples. The test dataset contained 400 examples. All models were trained for 100 epochs.
Helmholtz oper-
E DISCRETIZATION INVARIANCE
FNO intrinsically leverages the Fourier basis, allowing models trained on lower-resolution datasets to seamlessly process high-resolution inputs. By incorporating appropriate finite element basis func- tions, MgNO, compatible with the finite element method, demonstrates an equivalent capability. MgNO can train on lower-resolution datasets and evaluate higher resolutions without the neces- sity of high-resolution training data, thus accomplishing zero-shot super- resoutiino.u
FNO
MgNO
Test
128
256
512
128
256
512
Train
5.2808
7,9260
1054
1.032
1.178
1.205
Table 6: Comparison of discretization invariance property for MgNO and FNO for the Darcy multi- scale benchmark. The relative L,2 error (x10-2) with respect to the reference solution on the testing resolution is measured.
MODEL CONFIGURATIONS
We detail the primary configurations of the models in this section. Taking the NUMBER oF IteRATIOns LEVEL an example, the format [[..I/D..III.,,,II,,,I...CC the first level. the representation [1, 1] implies one a priori iteration followed by one post iteration (as indicated by the two values separated by comma).
is adopted. In the context of
19

---

Published as a conference paper at ICLR 2024
Modules
Darcy smmotthnougy/mmlliaaall
Navier Stokes
Helmhulte
NEHHER CHANNELS
12424,24212424
922.22.22.22.221 (1000011
28266122225.5121
ITERATIONS FER LEVEL;
Wate"
ConvouTTINN A- CONMELUTION
3 3x3, zmas FADDING
3.3 CIRECLAR PARDING
3 3. KIPLECT PADRING
333,
3
ConvoutTION
BARDING 1,
CIRCULAR 1,
3, STRIDE
CONTRANNO P.
X No FADDING, STRIDE No
4x 4, NO TARDENG, STRIDE YEs
4 4, NO PAROONN, STRIEE Yes
LAVER NORUALIAAAINN:
4/4/5
ACTIVaTIOO FUNCTIN:
Table 7: Model Configurations
20


Reviews:
[{"review_id": "qWMNqrQqS0", "reviewer": "Reviewer_5jaY", "summary": "This paper proposes a new neural network architecture based on analogies between convnets and the multigrid method (these analogies have been previously identified with an architecture called MgNet). The authors prove a universality theorem for a network whose neurons are operators and use this theorem to motivate a neural operator which is a deep network with neurons being (elliptic) operators implemented via multigrid. Numerical experiments show competitive performance against existing neural operator architectures.", "questions": "- Why is there no comparison with FNO in Figure 7 (Helmholtz)? It might be nice to compare / combine with https://arxiv.org/pdf/2301.11509\n- The sentence \"Regarding our training setup, we adopted a meticulous approach.\" is strange.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- The paper is very clearly written, the ideas are solid and the numerics are strong. I like the spirit of an architecture where neurons become operators: this appears very natural and the right thing to do.\n\n- Using a multigrid (or in general multiscale) representation of operators is also the right thing to do. Similar ideas go back to efficient wavelet approximations of operators of Beylkin, Coifman, and Rokhlin. It is nice to see this done explicitly and clearly it yields strong performance.\n\n- The experiments are very well executed.", "weaknesses": "- Many of the ideas might have already been present in MgNet; it would be nice to see a comment.\n\n- It is not clear to me what (9) is stating since you don't state the typical values of c. How should we use / interpret (9)?\n\n- The motivation via exact encoding of boundary conditions is sound but I would also say quite easy to implement with other neural operators. For example the FNO by default works on the torus but via zero padding could be easily adapted to Neumann or Dirichlet boundaries (though not as easily as MgNO).", "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": "", "extra_scores": {}}, {"review_id": "JK2mVhiFel", "reviewer": "Reviewer_H6UK", "summary": "The authors present a concise neural operator architecture designed for operator learning. This architecture draws inspiration from a conventional fully connected neural network, with each neuron in a nonlinear operator layer being determined by a bounded linear operator connecting input and output neurons, in addition to a function-based bias term. To address the efficient parameterization of these bounded linear operators within the architecture, the authors introduce MgNO (Multigrid Neural Operator). MgNO utilizes multigrid structures to effectively model and approximate these linear operators, eliminating the need for conventional lifting and projecting operators while accommodating diverse boundary conditions. Numerical experiments highlight the advantages of MgNO over CNN-based models and spectral-type neural operators, demonstrating improved ease of training and reduced susceptibility to overfitting.", "questions": "1. The proposed method may not scale well to larger datasets or more complex problems. Any remark for this?\n\n2. Is it possible to extend the universal approximation theorem to encompass general Banach spaces for X and Y?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The authors introduced a novel formulation for neural operators that characterizes neuron connections as bounded linear operators within function spaces. This eliminates the need for traditional lifting and projecting operators, simplifying the architecture.\n\n2. The MgNO architecture, which leverages multigrid structure and multi-channel convolutional form, efficiently parameterizes linear operators while accommodating various boundary conditions. This approach enhances both accuracy and efficiency.\n\n3. Empirical evaluations demonstrate superior performance of MgNO across multiple partial differential equations (PDEs), including Darcy, Helmholtz, and Navier-Stokes equations. This indicates the effectiveness and versatility of the proposed approach.", "weaknesses": "The limitations and drawbacks of the proposed methods are not explicitly mentioned", "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": "", "extra_scores": {}}, {"review_id": "dZsnNG6n5G", "reviewer": "Reviewer_RPM7", "summary": "This paper proposes a new neural operator architecture that does not require lifting and projecting operators. It uses multigrid structures (V-cycle multigrid) to parameterize the linear operators between neurons. Then this paper proves the universal approximation of their proposed parameterization and uses experiments to show the efficiency and accuracy of the method.", "questions": "1. What is the dependence of $n$ with $\\epsilon$ in Theorem 3.1?\n2. Why are the approximation capabilities so important? It might be a necessary condition for neural operators to learn but is definitely not sufficient. In my understanding, FNO is good at working with different resolutions which probably is not caused by its approximation capability.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The neural operator proposed in this paper uses the multi-scale method to get a better parameterization of the weights. It is applied to the spatial domain hence it seems like it won't significantly increase the complexity. From the experimental results, MGNO works quite well compared to other models, which provides convincing evidence that MGNO can be useful.", "weaknesses": "1. I am confused about the discretization of the input. It seems like the input still needs to be discretized (Eq.(6)) and the multigrid is applied to the spatial domain. Then I suppose the performance and the efficiency of MGNO are affected by the mesh size $h=1/d$. The performance of MGNO with respect to $d$ is missing in the paper.\n\n2. What about the parameterization with multi-channel inputs? Can this V-cycle be applied to multi-channel inputs? If there is no easy extension, I believe the significance of this work is compromised.", "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": "", "extra_scores": {}}, {"review_id": "uHBVFh2Vor", "reviewer": "Reviewer_4o6c", "summary": "In this paper, the authors introduce a novel formulation for neural operators and establish (in Theorem 3.1) the corresponding universal approximation property.\n\nThe authors then propose a multi-grid-based parameterization of linear operators, called MgNO, that can be efficiently parameterized and naturally accommodates various boundary conditions.\n\nNumerical results are provided on several popular PDEs, including Darcy, Helmholtz, and Navier-Stokes equations, with different boundary conditions, showing the superiority in prediction accuracy and efficiency of the proposed approach.", "questions": "I do not have specific questions but the following general comments for the authors:\n\n1. just being curious, can something similar to Theorem 3.1 be said about deep networks as defined in Section 3.2? How the network depth may play a role here?\n2. it would be helpful to discuss also the limitation of the proposed approach.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper provides novel insights and novel methodology of neural operators.\nThe numerical results on different popular PDEs with various boundary conditions look compelling.", "weaknesses": "I do not see particular weakness for this paper. See below for some comments.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)