Understanding the Expressivity and Trainability of Fourier Neural Operator: A Mean-Field Perspective
Authors: Takeshi Koshizuka
Masahiro Fujisawa
Yusuke Tanaka
Issei Sato
Venue: NeurIPS 2024 poster (2024)
Keywords: Fourier Neural Operator
PDE
Mean-Field Theory
Primary area: machine_learning_for_physical_sciences

Abstract: In this paper, we explores the expressivity and trainability of the Fourier Neural Operator (FNO). We establish a mean-field theory for the FNO, analyzing the behavior of the random FNO from an \emph{edge of chaos} perspective. Our investigation into the expressivity of a random FNO involves examining the ordered-chaos phase transition of the network based on the weight distribution. This phase transition demonstrates characteristics unique to the FNO, induced by mode truncation, while also showcasing similarities to those of densely connected networks. Furthermore, we identify a connection between expressivity and trainability: the ordered and chaotic phases correspond to regions of vanishing and exploding gradients, respectively. This finding provides a practical prerequisite for the stable training of the FNO. Our experimental results corroborate our theoretical findings.

Understanding the Expressivity and Trainability of Fourier Neural Operator: A Mean-Field Perspective
Takeshi Koshizuka
Masahiro Fujisawa RIKEN AIP masahiro. fujahiro.fujlssxx@rikenn jp
Department of Computer Science
The University of Tokyo
toshizuka-ttkeshi99888 9334444499984440gg..mm ecc. ecc. u-tokyo. jc.
Yusuke Tanaka
Issei Sato
NTT Communication Science Laboratories
Department Computer Science
ysk. tanaka@ntt.com
The University of Tokyo satołg. ecc. u-tokyo. ac jp
Abstract
In this paper, we explores the expressivity and trainability of the Fourier Neural Operator (FNO). We establish a mean-field theory for the FNO, analyzing the behavior of the random FNO from an edge of chaos perspective. Our investigation into the expressivity a random FNO involves examining the grdered-chaos phase transition of the network based on the weight distribution. This phase transition demonstrates characteristics unique to the FNO, induced by mode truncation, while also showcasing similarities those of densely connected networks. Furthermore, we identify connection between expressivity and trainability: the ordered and chaotic phases correspond regions of vanishing and exploding gradients, respec- tively. This finding provides practical prerequisite for the stable training of the FNO. Our experimental results corroborate our theoretical findings.
Introduction
The recent surge in interest in solving partial differential equations (PDEs) has led to the of neural network bNN)-based surrogate models. One promising line of work is the neural operator (NO), which learns the solution operator of PDEs, thereby bypassing the need for mesh dependency. Among the variants of NO, the Fourier neural operator (FNO) (Li et al, 2020c) has gained popularity because of its ad advantggcous cost/accuracy trade-off. The FNO can capture long-distance distance spatial interactions using the Fourier transform, whereas convolutional neural networks (CNNs) (Wen et al. 2019 Jiang al. 2021b) and message-passing graph neural networks (GNNs) (Li et al. 2020a b) are limited to operating solely on local variables. From a computational cost perspective, the Fourier transform is performed in quasi-linear time by the fast Fourier transform (FFT), making it significantly faster than the Transformer (Li et al. 2023).
Despite the widespread use of FNO as architecture, there is a lack of comprehensive theoretical analysis on its expressivity and trainability. The universal approximation property (Kovachki et al. 2021), recognized the basic expressivity of the FNO, is well-known however, the exponential expressivity depending on the weight distribution, which are known for the densely connected network (DCN) (Schoenholz et al.) 2016), a.k.a. fully connected network, and CNN (Xiao al.) 2018), remains unexplored. Regarding the trainability of FNO, the training instability in deep FNO has been experimentally reported by Tran et al. (2022), but the causes and conditions of the difficulty have not been clarified either theoretically or experimentally.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

homogeneous
Ordered phase
o²
Backpropagation
Vanishing gradients
H
H(4)
the edge of chaos
Exploding grudients
heterogencous
H
Chaotic phase
H
H
Figure 1: Illustration of ordered-chaos phase transition for the weight initialization parameter o². In the ordered phase, the spatial hidden representations H on the grid converge uniform state during forward propagation and the gradient vanishes during backpropagation. In the chaotic phase, the representations either converge a distinct state or diverge and the gradient explodes.
We analyze the exponential expressivity (how far can two different input vectors be pulled apart) and trainability (how much gradient explosion average) of the random FNO from the perspective of whether the network is ordered or chaotic. This viewpoint is grounded mean-field -eead theory, an analytical framework for NN established by Poole et al. (2016); [Schoenholz al. (016): Yang & Schoenholz (2017); Xiao et al.] (2018). A network is considered ordered when brings all representations of different spatial positions closer together, and chaotic when it drives them apart during forward propagation. Furthermore, a network can only be stably trained when initialized close the edge of chaos, which the transition point between the ordered phase and the chaotic phase. In fact, He initialization (He et al. (015)) is an example of a commonly used edge of chaos initialization for the DCN with ReLU activation (Burkholz & Dubatovkal 2019).
In this study, we establish mean-field theory to analyze the expressivity and trainability of the FNO. Our investigation reveals the expressivity of random FNO at initialization by examining the transition point between ordered and chaos phases. The phase transition exhibits -speccccciic characteristics induced by mode truncation, well as similarities with the characteristics of DCN and CNN. We also find a link between expressivity and trainability: the ordered and chaotic phases correspond to regions of vanishing and exploding gradient, respectively. This discovery offers practical initialization prerequisite for the stable training of the FNO.
Background
2.1 Fourier Neural Operators
The FNO (Li et al. scientific problems Yang al. 1 Wen et al.] 2022b Hwang al. 2022 Jiang al.) 2021a Pathak et al. 2022). An M -ddimensional FNO with the number hidden features /) and the spatial size N for learning the operators between scalar-valued functions is defined as follows.
2020c) is one of the well-established methods for solving PDEs across many
X (e+1) & (D (x() + (()) (x()))
where (0 R N is the l-th hidden representation, M is the number of spatial dimensions, and ø is activation. The hidden representations X and (T) are the output of the lifting operator and the input of the projection operator, respectively. The architecture of these operators does not affect analysis long the network stays shallow, implemented in (Li et al. (020cc). The E-th densely connected (DC) module ()() is an affine point-wise map in the physical space and the L-th

---

Fourier convolution module ((()) is a parameterized kernel integral operator using Fourier transform. The bias term is considered to be included in either or both modules K(()) and D(
Several FNO variants have been developed to address specific challenges, such geo- FO- (Li ct al. 2022) for irregular regions and group equivariant FNO (G-FNO) (Helwig et al. 2023). which
maintains equivariance to rotation, reflection, and translation. U-NO (Rahman al. 2022 and
U-FNO (Wen et al., 2022a)) integrate FNO with U-Net for multiscale modeling. Additionally,
WNO TTipura & Chakraborty 2022 utilizes wavelet bases, while CFNO (Brandstetter al.
2023) enhances the use of geometric relations between different fields and field components through Clifford algebras. Adaptive FNO (Zhao et al. 2022 Guibas et al. 2021) and F-FNO (Tran et al. 2022) have improved computational and memory efficiency through incremental learning and architectural modifications. Other approaches for improving performance include methods with increasing physical inductive bias (Li et al. 2024), data augmentation (Brandstetter et al. 2022), and a criiacc--rreerving weight initialization scheme (Poli al. 2022).
While numerous new models and learning methods have been proposed, relatively little research has been conducted to understand the intrinsic nature of these methods. Issues such spectral bias (Zhao et al. 2022) and training instability (Tran et al.) (022) have been reported. [Tran et al. (2022) observed that training did not converge even at 24 layers. They successfully addressed the stability and accuracy degradation issues associated with an increase in the number of layers by implementing skip connections behind activation and introducing various training techniques. However, it is still unknown that the theoretical basis for why the original architecture of the FNO has problems with training instability and accuracy degradation.
2.2 Mean-field Theory for Neural Networks
The mean-field theory has been used to provide a mathematical framework for understanding the expressivity and trainability of neural networks (Poole et al., 2016; Schoenholz et al. 2016 Yang Schoenholz] 2017 Hayou et al. 20188) Xiao et al. 2018)) A series of papers (Poole et al. 2016 Schoenholz al. 2016) delved into the average behavior of infinite- random deep DCN, with weights and biases initialized by zero-mean Gaussian distribution. The formulation is given below.
10
(0)
(0
(0 x (e-1) +b
w rid. N
0,
b rid. N (0,0%)
where xe (*) E RR. the f-th hidden representation, (()) E R DxD ((9) E RR." are the f-th learnable parameters, and the width is assumed to be sufficiently large D > 1.
Poole et al. (2016) and [ccoennolzz et al. (2016) explored the exponential expressivity of random DCN determined by two phases depending on the initial variance parameters o² and oū. as shown in Fig. 5a| Poole et al. (2016) first examined the forward propagation of random DCN with Tanh activation. They demonstrated that the covariance X of the ε-th pre-activation representations (()
and ( corresponding two different inputs xr and (0) are obtained by
Ad E [D], { = o²R
[ ( (n( & (NC (R
+ +
N((-1)) ~ AO 1)) The where the expectation is taken over the pre-activations [hd (e-1) [h(e-1) hid covariance converges exponentially to a fixed point Σ determined by parameters o² and of
A network is considered ordered when it brings two distinct representations closer together, which implics a state of small expressivity. Conversely, network is chaotic when it drives them apart during forward propagation, implying state of large expressivity Networks with either excessively small or large expressivity can disrupt the structure of the input: the difference between two distinct inputs quickly becomes indistinguishable in networks with small expressivity, while similarities between inputs no longer recognized in networks with large expressivity. The network is ordered the initial variance of the weights is small. For larger values, and beyond a certain threshold, the phase shifts, and the network behaves chaotically. This phase shift point is termed the edge of chaos.
Subsequently, [choenholz et al. (2016) discovered the connection between expressivity and trainabil- ity of DCN through analysis of the backpropagation behavior. In an ordered network, the expected
3

---

value of the gradient norm becomes exponentially small during backpropagation, while it becomes exponentially large in a chaotic network. This implies that the gradient vanishes/explodes in ordered or chaotic networks, respectively. These findings suggest that deep DCN can be stably trained only near the edge of chaos. [cchoenhll et al. (2016) also provided an estimate of the maximum depth at which network can be trained when initialized away from the edge of chaos. These insights are not limited to DCN and similar results have been observed for residual networks (Yang & Schoenholz 2017) and CNN (Xiao et al. 2018).
3 A mean-field theory for FNO
In this section, establish mean-field theory of FNO. We demonstrate the exponential expressivity of random FNO by examining the ordered-chaos phase transition during the forward propagation. Furthermore, we identify the connection between expressivity and trainability by concentrating on backward propagation behaviors. Our analysis is an advanced version of the approach developed
by Poole al. 6) Schoenholz et al. (2016); Yang & Schoenholz (2017); Xiao et al. (2018). In Section 3.1 we outline the problem setup. In Section 3.2 we analyze the forward and backward
propagation behavior of random FNO at initialization. In Section [3.3] we discuss the practical
prerequisites for initialization to stabilize the training of FNO, leveraging the similarities between
FNO and DCN. The proofs for the lemmas and theorems are provided in Appendices A and
3.1 Problem setting
Here, we consider a simplified one-dimensional (ID) FNO. Note that our theory is extensively
applicable to the original FNO, as discussed in Section 3.3 The simplified ID FNO, with depth
of L, is defined by the number of hidden features D, spatial size N 2 (where m is an integer), the number of Fourier modes K S N + 1, two learnable weights (t.A)) E IR DxD and ((t,)) DxD, and bias ( RD. Denote Φb: R IR by the non-decreasing activation function. Let (() NXD and H E NXD be the post and pre-activation representations defined by
K-1
( = Φ (()) H
Σ
Ck
H
(E,R)
(((b))
+ 1
2
(3)
((M))) = F (1k) FX ((-1)
(0)( + V=
where dabb is the Kronecker-delta, Cs 2 5x,0 0g,N/2 /2 is constant, 1 is all-ones column
vector with the size N, H ((,k) is the conjugate of (t.k) corresponding the (N k-)-th frequency components, it is the transpose conjugate, F E NNN is the Discrete Fourier Transform (DFT) matrix defined by Fk.n N/2 exp( 2
N n), and (()) (k) a diagonal matrix with a at position ((A) k,k-
There are two differences from the original FNO proposed by Li et al. 2020aa: (1) the DC module is dropped for the simplicity. and (2) (f. is multiplied by v2 with respect to k 0, N appropriate normalization. We assume that the weights of FNO are initialized by i.i.d. samples from
z for
Gaussian distribution, i.e. (f.k) i.dd. NO $)) E.
((.kk i.idd NN(0, ei) ) b(0) i.dd N(0, For 2 the parameter = (E,k) is set to zero exceptionally. For all d E [D] = ...... D - 1}. the pre-activations (9)) E TR are i.i.d. random variables. When > 1, by the central limit
i.j
k = 0, N
..d
theorem, the variables ( follow Gaussian distribution with mean 0 and covariance matrix x)((
(f) Eo 14 a'.d where the expectation is taken over all random variables [0 E
( M) (((,k') )e Our theory can be easily extended to 2D and 3D FNOs.
3.2 Expressivity and trainability of FNO
Firstly, the forward propagation of single input signal with spatial features is described as follows.

---

Lemma 3.1 (Iterated map). For all [D), the covariance E Ed NiguitE.1 gad
H
H
(6)
is
.
obtained recursively by the iterated map C defined by
K
(0
Ca
|IF
(H. allel
cos (0ae)) + o.
(4)
0,00
k==0
((E((11) Ja,a'
where the expectation is taken over the pre-activations H N(,, ((-))) gl0) a.a' = represents the scaled positional difference.
2 N (a a!)
The indices a and a' correspond to different spatial locations with the mean-field theory for CNN (Xiao et al. 2018). Note that F6 (H..a) ..]] is the k-th Fourier modes of the post-activation representation. When applying DCN (Poole et al. 2016 Schoenholz et al., 2016) CNN (Xiao et al., 2018) to the spatial signal, the iterated map depends only on local spatial locations, while in the case of FNO, the iterated map depends all spatial locations because of the global Fourier convolution. In addition, only periodic spatial correlations with iift--ivariiat are propagated, and high-freeueen components exceeding mode K are eliminated.
Next, we explore the fixed point E' of the iterated map C satisfying E' = C(E*). By linearizing the dynamics of signal propagation around this fixed point and analyzing the stability and rate of convergence to the fixed point, we can determine the depth to which each component of the input can propagate. [Schoenholz al. (2016) showed that the iterated map of DCN defined in Eq. 2 has fixed point of the form:
E' -"= - (INN1 - IN),
(5)
where q²' ce are the fixed points of variance and correlation, and IN is the identity matrix. Meanwhile, Xiao et al. (2018) showed that any fixed point for the iterated map the DCN is also a fixed point for that of CNN. We show that random FNO has the same fixed points of the form of Eq. (5) with
e" = in the following lemma.
Lemma 3.2 (Exsistance of fixed points). When random DCN defined by Eq. 2 has the fixed point (g" e* ===11 1) for the initial parameters (o²,al) o%), then random simplified FNO defined by Eq. (3) has
a fixed point E of the form
E' ='=q'IN + (1N1 IN) 1"1N1N
Lemma 3.2 indicates that the fixed point for the iterated map Σ* of the DCN serves as a fixed point for the iterated map of the simplified FNO (as well as CNN). To analyze the stability and convergence rate, we linearly approximate the C-map around the fixed point Σ*, i.c., C(E) se E' + + (E - E"), where Jz- the Jacobian linear map of the iterated map defined Eq. 15 We then derive the eigenvalues and eigenvectors for the Jacobian linear map Js as follows.
Definition 3.3.
Xy o'' ("H
(6)
Xe (HH (Har.a),
(7)
0² Xx = of E. B
+ ol a) (H (Ha),)) .)) + o R (c (o (Ha.a) (H
(8)
where the expectation is taken over the pre-activations H. ~ N(0, E*), and Φ, o" are the first- and second-order derivatives of the activation Ø.
The bases sb, (1),, (k1) E NxN using the above quantities are defined below.
20, =1
Xx + Xe
Xg:
Σ C2 cos
(g(+)
N
Xx
3,9'
(9)
K
(K)
(0(k)
(0($)
8.9'
cos
8,3'
Σ
ΣØ
k
C₁ ==0
Co cos
3,9'
5

---

From Lemma A.4, K - 1 matrices in {25(4) [eE are eigenbases with the eigenvalue Xe- of the Jacobian linear map. From Lemma A.5, the matrix & the eigenbases with the eigen- value X of the Jacobian linear map. Since the rank of the Jacobian linear map is at most K (Lemma A.3), the deviation from the fixed point ) - E is spanned by K-dimensional eigenspace span (fab(")) (EE[A] U {œ))) Then, the fixed point stability and the convergence rate are shown in the following theorem.
Theorem 3.4 (Exponential expressivity). Let E = D E' be the deviation from the fixed point at the f-th layer. Suppose that the deviation at the first layer is decomposed as E = a EK=1 exap(h) ee. The scalars e, Ej represent the scale of the perturbation for each eigencomponent of the linearty approximated map E(d) + E((+1), The component e E RNNN XN belongs to the
orthogonal complements of the space span (tab, sb(1), ..... 11k-1))).
Then, the deviation at the E-th layer is obtained by
E
=xEe
Σ
X Ex2b(1))
(10)
K-1
K-1
CxX +
1 (1-
Gs
(Xr + Xc+).
N
8=0
In particular, when the Fourier mode is K
+ 1, Eqs. and 10
99
reduce to the following.
K-1
E ex Σ
X c₁ (b())
k=1
(11)
v3, [N], 63, 1, 4 cos (0%/3) 6g.
Theorem B.4 shows the expressivity of the FNO, which is characterized by the ordered-chaos phases and varies exponentially with respect to the number of layers. Theorem 3.4 indicates that the asymptotic behavior of the zero-frequency deviation is mostly determined by X and the periodic deviation is determined by Xc. If < 1 and Xe- < 1, the fixed point is stable as the deviation from the fixed point converges exponentially to zero. When the fixed point remains stable at c* 1, random network exists in ordered phase, where all spatial representations are correlated in an asymptotic manner. Conversely, when the fixed point with e* = becomes unstable, the network transitions into a chaotic phase, exhibiting behavior dependent on the activation function Φ. The boundary between these two phases is referred to as the edge of chaos.
The convergence rates X4- and Xe. are the same as the convergence rates of the variance and correlation to the fixed point for DCN (Schoenholz al., 2016) and CNN (Xiao et al. 2018). However, only periodic spatial correlations are propagated in the FNO, resulting in different eigenspace of the map E(() + E({+1) compared to the DCN and CNN. In DCN, the deviation belongs to vector space with dimension N(N-1) in DCN, whereas in FNO, the dimension is K, or most NN + 1. CNN possess diagonal eigenspaces associated with eigenvalues Xo- and non-diagonal eigenspaces associated with eigenvalues Xe In contrast, FNOs without mode truncation exhibit similarity, possessing eigenspaces Xg- for pero-freyynccy and eigenspaces Xe- for k-frequencies with diagonal components removed. Furthermore, mode truncation increases the convergence rate of cero-frequency deviation from Xg- to X and affects all eigenbases as well. For further discussions on the similarities between CNN and FNO, please refer to Appendix C A visualization of the covariance of the FNO with Tanh and ReLU activations shown in Appendix F
Finally, we demonstrate the connection between expressivity and trainability. By examining the covariance of the gradient each layer during backpropagation, we investigate the conditions under which training is stable without gradient vanishing or exploding.
Theorem 3.5 (Trainability). Let ÊO RNNN be the gradient covariance with respect to some loss C. e.g. mean squared error, at the b-th layer. Suppose that the gradient covariance at the L-th layer decomposed as ((L) a,a' = Σk= EK cos (o(t) + ē, where E is the coefficient of each basis and

---

(a) Simplified FNO with Tanh (b) Simplified FNO with ReLU
(c) Original FNO with ReLU
Tr (r(D during the backpropagation of several FNOs plotted as a function of layer &. Each line corresponds to the result of different initial values of x² from 0.5 to 4.0 in increments of 0.5. The x-axis the layer and the y-axis the log-scale of the gradient norm. Depending on the value of o², the gradient norm increases or decreases consistently as the gradient propagates shallower layers.
Figure 2: Average gradient
belongs to the orthogonal complements of span((cos (Da.d) ).oe)> 3 Then, the gradient covariance at the E-th layer is obtained by
K-1
550 a.a'
Σ
Xe E cos (satdr)
--t
km0
Theorem 3.5 shows that gradient vanishing occurs when Xe < 1 (ordered phase) and gradient explosion occurs when Xe- (chaos phase). Thus, stable training of the FNO can be achieved close to the edge of chaos by setting the initial parameter o² to satisfy Xe se 1. We specifically present the initial parameter choices that achieve the edge of chaos for several FNOs in Section
When
there is no change in the dynamics during backpropagation due mode truncation.
When using the full mode K -/2 + 1, the condition X.- = always achieves the edge of chaos, which is consistent with the results for differences among FNO, DCN, and CNN, Theorem 3.4 and Theorem 3.5 demonsttatee the similarities the random behavior of FNO, DCN, and CNN. This allows existing results based on mean-field theory to be applied to the FNO.
DCN and CNN. Despite architectural and iterative map
3.3 Initialization requirements for stable training
For stable training, Theorem 3.5 suggests the necessity of initializing FNO near the edge of chaos, i.e., initializing FNO so that X the edge of chaos for several FNOs, each with slightly different architectures such as activation functions. Furthermore, the behavior of the gradient norm Tr (((())) D function of layer l are visualized Fig. 2 for several variants random FNO with different initialization parameters o². We used FNO with a width of D 32 and a number of layers I 64, and for simplicity, we computed the absolute value the output as the loss with respect to the input sampled from the standard normal distribution.
this section, we present the initial parameter choices that achieve
Simplified FNO with Tanh activation.
The behavior of Xe for the parameters o² and of of the Tanh network has been extensively studied by Poole al.) (2016): Schoenholz et al. (2016). The phase diagram drawn by [Pennington al.) (2017) is shown Fig. 5a By using parameters (o² ( around the two phase boundaries of ordered and chaotic that achieve Xe* = the training of the simplified FNO with Tanh activation can be stabilized. Figure 2a depicts the behavior of the gradient backpropagation in the simplified FNO with Tanh activation and the bias parameter being of 0.1. Figure 2ashows that when o² x 2, the gradient diminishes exponentially; otherwise, it explodes exponentially.
7

---

Simplified FNO with ReLU activation.
The iterated map C of the DCN with ReLU activation is given by Cho & Saul (2099) follows.
q(tt1)
1 + +ob.
02)(0)
ob.
(12)
2
(
1
0²
cl0)g()) g(0)
+ to
(13)
2
(i(c) = (V1 - d² (क arccos (c))e).
The edge of chaos initialization for the DCN with ReLU activation is known as He initialization (He al. 2015), which sets the initial variance parameter as x² 2 and initial bias as b(d) for Eq. (2) From the similarity of the DCN and the FNO, we can derive the FNO version of the He initialization that achieves Xc. 1 by setting o² 2, 6(()
for Eq. ( The He initialization scheme for the
simplified FNO with the activation Φ = ReLU is derived as follows.
(2.R) Lid. NNOO ( i.dd N D ) E
Figure 2h demonstrates that the choice of a² 2 preserves the magnitude of the gradient norm during backpropagation of deep simplified FNO with ReLU activation.
Original FNO.
In the original architecture of the FNO proposed by Li et al. (2020c), the DC module used together with the Fourier convolution module shown in Eq. Qpp. We initialize the weights of both layers consistently as follows. For all E [L]. l E [K], and i.3 E [D].
( iid. N (0.
0²
(2,K) Lid. N 0, 0²
4D
4D
W Lid. N
0, 2D
b. N (0. o
From the similarity of the initial network behavior of the FNO and the DCN, the fixed point with c" ===1 1 of the simplified FNO is also fixed point of the original FNO. In the neighborhood of the fixed point Σ', the eigencomponents spanned by {2)(() =k=1 will decay increase at the rate of Xe* Following the derivation of Theorem [.3] the eigenvalues of the cos function eigencomponents of the gradient covariance are also Xe.. These results show that the edge of chaos initialization scheme can be used for the original FNO with each activation function. Figure (2c shows that the original FNO with ReLU activation and the parameter fixed as 6(() = 0 exhibits similar backpropagation behavior as the simplified FNO, i.e,, σ² = is an appropriate choice.
Experiments
In this section, experimentally demonstrate that deep FNO training requires appropriate initializa- tion settings on a variety of datasets, consistent with the theory discussed Section
4.1 Datasets
We evaluated three models on commonly used PDEs: the advection equation, Burgers' equation, Darcy Flow equation, and incompressible Navier-Stoke (NS) equation. All datasets were generated by numerical simulations used in (Takamoto et al. 2022 Li et al. 2020c) and are publicly available.
summary of the dataset provided Table and more details are given in Appendix [D]
Advection equation and Burgers' equation. The linear advection equation for the function u(x; t) given by
0,u(x, 1) + 30,, (u(x, ()2) 0, u(x, 0) = uo(x),
where un is the initial condition and 8 = 2.0 is an advection speed. The non-linear Burgers' equation for the function u(x, 1) is given by
0,((,,)) t) + O2 (u (r, 2)) viluu(x, f), u(z, 0) = ux)
8

---

Table 1: Overview of dataset with number of spatial dimensions M, time dependence, spatial resolution N₂ = x ... . x N, temporal resolution N₁, and number of samples for training, validation, and testing.
PDE
M
time
N₂
N₂
Number of samples (Train / Val / Test)
advection
1
64
8000 1000 1000
Burgers'
1
64
8000 1000 / 1000
Darcy flow
2
64 64
900 100 / 100
Navier-Stokes (v le-3)
2
V
64 x 64
25
1000 / 100 100
Navier-Stokes (v/ le-4)
2
64 64
20
8000 1000 / 1000
Navier-Stokes (v le-5)
2
64 x 64
20
1000 / 100 100
where Up is the initial condition and V = 4.0 the diffusion coefficient.
Darey Flow equation. The Darey Flow equation for the function u(z) with Dirichlet boundary is given by
V (a(x)))u((z)) = f(x) (x E (0,1)
u(x) =0 (ar E 0(0,1)3)
where a the diffusion coefficient and f(x) is the forcing function.
Incompressible Navier-Stokes equation. The 2D NS equation on the unit torus is defined by
a,w(r, t) + u(x, t) Vw(x, t) = wV, t) + f(z), t)) 0, wxx 0) wu(x),
where w(x,t) t) is the vorticity, wo the initial vorticity, u(x,t) the velocity field for any > 0, V R4 is the viscosity, and f is the external forcing function defined by f(x) = 0.1 (2in(2 + x22 + cos(2m(11 + x2)). We experimented with the viscosities v =
le 3, 1,l.,,1e---- For the data with v 1c-4, the time resolution also downsampled by half.
4.2 Experimental Settings
In the experiments the 1D advection and Burgers' equation, compared the results of simplified FNOs defined in Eq. (3) with Tanh and ReLU activations for varying number of layer and initial parameters o², In the experiments on the 2D Darcy Flow and NS equations, we compared the results of the original FNO et al. 2020c) as shown in Eq. (T) for varying number of layer L and initial parameters 2²_ All models have width of D 32 and the Fourier mode of K 12. When using Tanh as the activation, we fixed the initial bias parameter of = 0.1, and when using ReLU activation,
We used the Adam W optimizer and cosine annealing scheduler. Training was stopped early at the epoch of minimal normalized mean squared error (nMSE) on the validation data. Details of the experimental setup are given in Appendix [D] the experiments on NS equation, we trained the original FNO with the autoregressive scheme using a teacher-forcing technique, input normalization, and regularization that adds small-Gaussial noise to the input (Tran et al. 2022). During the evaluation phase, only the initial 10 steps are provided as input, and the rollout results of all subsequent steps are evaluated. For all tasks, the mean squared error (MSE) is used the training loss and nMSE as the validation and testing (Li ct al. 2020a Tran et al. 2022
we fixed the initial bias (/")
4.3 Results
Heatmaps of the training loss measured by MSE at the last epoch and the test performance measured by nMSE on six different datasets, for different depths L and initial parameters a², are shown respectively. Despite the differences in architectures and datasets, Figure shows the same trend supporting the theory in all experiments. As the number of layers I. increases, the range of acceptable initial o² value settings becomes narrower, and initialization near the edge of chaos (o2 x 2) is required for stable training of deep FNO. Detailed analyses on training loss and test performance are presented in Appendix E. and Appendix E.2 respectively.
in Figs. 3] and 4
5 Conclusion
In this study, we developed mean-field theory for FNO. We showed the expressivity and trainability of the FNO, which is characterized by the ordered-chaos phases. Furthermore, we observed both

---

16-
201
109
10
12
(a) Tanh advection eq. (b) ReLU advection eq. (c) Tanh on Burgers' eq.
(d) ReLU Burgers' eq.
7
ReL.U Darey Flow (f) ReLU NS eq. (Ic-3) (g) ReLU NS eq. (1c-4) (h) ReLU on NS eq. (1c-5)
Figure 3: Training loss of FNOs at last epoch for four distinct PDEs. (a, b): the advection equation, (c, d): the Burgers' equation, (e): Darey Flow, (f-h): the NS equation. The heatmaps represents the training loss values for varying depth L {4,8, 8, 16, 32) and initial weight parameter o² {0.1, 0.5, 1.0, 2.0, 3.0, 4.0), with lighter colors signifying lower training loss. The presented results the mean training loss the last epoch over three different seeds.
12
(a) Tanh on advection eq. (b) ReL.U on advection eq. (c) Tanh on Burgers' eq.
(d) ReLU on Burgers' eq.
16
40
ar
166
(e) ReLU on Darcy Flow (f) ReLU on NS eq. (Ie-3) (g) ReLU on NS eq. (le-4) (h) ReLU on NS eq. (le-5)
Figure 4: nMSE of FNOs on test datasets for four distinct PDEs. (a, b): the advection equation (c, d): the Burgers' equation, (e): Darcy Flow, (f-h): the NS equation. The heatmaps for each nMSE correspond to the results of each heatmap of training loss in Fig.3 The lighter colors, the better the test performance. The presented results are the mean nMSE calculated over three different seeds.
unique FNO-speciff behaviors caused by mode truncation, well as common behaviors akin to those DCN. With our analysis as a basis, we identified the necessity of initializing FNO near the edge of chaos for stable training of the FNO. Experimental results supported our theoretical results.
limitation of our analysis is that limited the network at initialization and does not address the stability of the entire optimization process. While we do not provide sufficient conditions for stable training, we offer one necessary condition for achieving stable training. Future work may consider a mean-field analysis of the FNO when using skip-connection (Tran et al. (022), Dropout and batch normalization, as well as initialization methods that ensure the input-output Jacobian of the FNO satisfies dynamical isometry (Pennington et al. 2017 2018).
10

---

Acknowledgement
This work was supported by JSPS KAKENHI Grant Number 24H00709 Japan, JST ASPIRE Grant Number JPMJAP2329, RIKEN Special Postdoctoral Researcher Program, JST ACT-X Grant Number JPMJAX21OK, and JST ACT-X Grant Number JPMIAX2IOD.
References
Brandstetter, J.. Welling, M., and Worrall, E. Lie point symmetry data augmentation for neural pde solvers. In International Conference on Machine Learning, pp. 2241-2256. PMLR, 2022.
Brandstetter, J., van den Berg, R., Welling, M., and Gupta, J. K. Clifford neural layers for pde modeling. In The Eleventh International Conference on Learning Representations, 2023.
Burkholz, R. and Dubatovka, A. Initialization of relus for dynamical isometry. Advances in Neural
Information Processing Systems, 32, 2019.
Cho, Y. and Saul, L. Kernel methods for deep learning. Advances in neural information processing systems, 22, 2009.
Guibas, J., Mardani, M., Li, Z., Tao, A., Anandkumar, A., and Catanzaro, B. Adaptive fourier neural operators: Efficient token mixers for transformers. arXiv preprint arXiv:2111.13587, 2021.
Hayou, S., Doucet, A., and Rousseau, J. On the selection of initialization and activation function for deep neural networks. arXiv preprint 1eX.::8805.08266 2018.
He, K., Zhang, X., Ren, S., and Sun, Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In Proceedings of the IEEE international conference on computer vision, pp. 1026-1034, 2015.
Helwig, J., Zhang, X., Fu, C., Kurtin, J., Wojtowytsch, S., and Ji, S. Group equivariant fourier neural operators for partial differential equations. arXiv preprint arXiv:2306.05697, 2023.
Hwang, R., Lee, J. Y., Shin, Y., and Hwang, H. J. Solving pde-constrained control problems using operator learning. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, No. 4, pp. 4504-4512, 2022.
Jiang. P., Meinert, N., Jordão, H., Weisser, C., Holgate, S., Lavin, A., Lutjens, B., Newman, D., Wainright, H., Walker, C., et al. Digital twin earth-coasts Developing a fast and physics- informee surrogate model for coastal floods via neural operators. In Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021), 2021a.
Jiang, Z., Tahmasebi, P., and Mao, Z. Deep residual u-net convolution neural networks with autoregressive strategy for fluid flow predictions in large-scale geosystems. Advances in Water Resources, 150:70 2021b.
Kovachki, N., Lanthaler, S., and Mishra, S. On universal approximation and error bounds for fourier neural operators. Journal of Machine Learning Research, 22((90):1-76, 2021.
Li, H., Xu, Z., Taylor, G., Studer, C., and Goldstein, T. Visualizing the loss landscape of neural nets. Advances in neural information processing systems, 31, 2018.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar,
A. Neural operator: Graph kernel network for partial differential equations. arXiv preprint ar arXiv:2003.03485, 2020a.
Li, 7.,, Kovachki, N., Azizzadenesheli, K., Liu, B., Stuart, A., Bhattacharya, K., and Anandkumar, A. Multipole graph neural operator for parametric partial differential equations. Advances Neural Information Processing Systems, 33:6755-6766, 2020b.
Li, Z., Kovachki, N. B., Azizzadeneshelii K., Bhattacharya, K., Stuart, A., Anandkumar, A., et al. Fourier neural operator for parametric partial differential equations, International Conference Learning Representations, 2020c.

---

Li, Z., Huang, D. Z., Liu, B., and Anandkumar, A. Fourier neural operator with learned deformations for pdes on general geometries. arXiv preprint mxXiv2200.005209, 2022.
Li Z., Meidani, K., and Farimani, A. B. Transformer for partial differential equations" operator learning. Transactions on Machine Learning Research, 2023.
Li, Z., Zheng, H., Kovachki, N., Jin, D., Chen, H., Liu, B., Azizzadenesheli, K., and Anandkumar, A. informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 1(3):127,, 2024.
Pathak, J., Subramanian, S., Harrington, P., Raja, S., Chattopadhyay, A., Mardani, M., Kurth, T., Hall,
D., Li, Z., Azizzadenesheli, K., al. Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv preprint arXiv:2202. 12214, 2022.
Pennington, J., Schoenholz, S., and Ganguli, Resurrecting the sigmoid in deep learning through dynamical isometry: theory and practice. Advances in neural information processing systems, 30, 2017.
Pennington, J., Schoenholz, S., and Ganguli, S. The emergence of spectral universality in deep networks. In International Conference Artificial Intelligence and Statistics, pp. 1924-1932.
PMLR, 2018.
Poli, M., Massaroli, S., Berto, F., Park, J., Dao, T., Ré, C., and Ermon, S. Transform once: Efficient operator learning in frequency domain. Advances in Neural Information Processing Systems, 35: 7947-7959, 2022.
Poole, B., Lahiri, S., Raghu, M., Sohl-Dickstein, J., and Ganguli, S. Exponential expressivity in deep neural networks through transient chaos. Advances in neural information processing systems, 29, 2016.
Rahman, M. A., Ross, Z. E., and Azizzadenesheli, K. U-no: U-shaped neural operators. arXiv preprint arXiv:2204.1112 2022.
Schoenholz, S. S., Gilmer, J., Ganguli, S., and Sohl-Dickstein, Deep information propagation. In International Conference Learning Representations, 2016.
Takamoto, M., Praditia, T., Leiteritz, R., MacKinlay, D., Alesiani, F., Pflüger, D., and Niepert, M. Pdebench: An extensive benchmark for scientific machine learning. Advances in Neural Information Processing Systems, 35: 15:1516-1,,,, 2022.
Tran, A., Mathews, A., Xie, L., and Ong, C. S. Factorized fourier neural operators. In The Eleventh International Conference on Learning Representations, 2022.
Tripura, T. and Chakraborty, S. Wavelet neural operator: a neural operator for parametric partial differential equations. arXiv preprint arXiv:2205.02/91, 2022.
Wen, G., Tang, M., and Benson, S. M. Multiphase flow prediction with deep neural networks. arXiv preprint arXiv:1910.09660 2019.
Wen, G., Li, Z., Azizzadenesheli, K., Anandkumar, A., and Benson, M. U-fno-an enhanced fourier neural operator-basee deep-learning model for multiphase flow. Advances in Water Resources, 163:104180, 2022a.
Wen, G., Li, 7.,, Long, Q., Azizzadenesheli, K., Anandkumar, A., and Benson, S. M. Accelerating car- bon capture and storage modeling using fourier neural operators. arXiv preprint arX1v:2210.17051. 3, 2022b.
Xiao, L., Bahri, Y., Sohl-Dickstein, J., Schoenholz. S., and Pennington, J. Dynamical isometry and mean field theory of cnns: How to train 10,000-layer vanilla convolutional neural networks. In International Conference Machine Learning, pp. 5393-5402. PMLR, 2018
Yang, G. and Schoenholz, S. Mean field residual networks: On the edge of chaos. Advances neural information processing systems, 30, 2017.
12

---

Yang, Y., Gao, A. F., Castellanos, J. C., Ross, Z. E., Azizzadenesheli, K., and Clayton. R. W. Seismic wave propagation and inversion with neural operators. The Seismic Record, 1(3):126-134, 2021.
Zhao, J., George, R. J., Zhang, Y., Li, Z., and Anandkumar, A. Incremental fourier neural operator. arXiv preprint arXiv:2211.55188 2022.
13

---

0.25
0.25
q
Ordered
1.5
0.20
c==1 Vanishing gradicnts
0.15
Ordered
Chaotic
1.0
8
6
bounded, =1
unbounded, c* =
Vanishing gradients
Exploding gradients
0.5
0.10
Chaotic
0.05
0.05
c' 1
0.0
Exploding gradients
0.00
0.5
2.0
2.5
3.0
0.00
0.5
1.0
1.5
2.0
2.5
3.0
0
Tanh activation (Pennington et
2017
ReLU activation
Figure 5: Ordered-chaos phase transition diagram for the DCN
A Proof of Theorem 3.4
E* be the deviation from the fixed point the l-th layer. Suppose that the deviation at the first layer is decomposed as E c ΣK=1 <xap()) The scalars e, 6k represent the scale of the perturbation for each eigencomponent of the linearly approximated map (
Theorem 3.4 (Exponential expressivity). Let E(A) = Σ(O)
E The component N belongs to the
orthogonal complements of the space span (fab.ebl) ..... (.(.,))))
Then, the deviation at the 8-th layer is obtained
E
x'ab
Σ
X (26((0)
(10)
K- N
K
x
CxXa"
+
1
Σ
(xx + X
N
8=0
==0
In particular, when the Fourier mode is K
+ 1, Eqs.
9 and 10
reduce to the following.
K-1
E = X=- Eap + Σ
Xe (b)
k=1
(11)
V3,8' E (N), U3 =
(k) 3/ = cos
(0(0)
03,31
63.9
Proof. The theorem is obtained by the eigenvalue analysis on the the first-order Taylor approximation of the iterated map C at the fixed point Σ'. The Jacobian matrix J(Σ*) ) E N x of the iterated
map at the fixed point and the Jacobian linear map Jzz() (-): NNN - R N are defined as follows.
(J(E)) (aaa))))) = |C(()))
(14)
σΣg,8'
E E-
VΣ E Js ( mat (J(E") vec(E))
(15)
where vee performs the vectorization, i.e. transforming an x matrix into vector of size N and mat performs the inverse operation of vec.
From Lemma A.3 of the Jacobian linear map. From Lemma(A.4 the matrix &/ is the eigenbases with the eigenvalue of the Jacobian linear map. Since the sets of K matrices in ) ]RE U {2b are lincarly
( matrices in { (vec(1b(0) (4)( } (0) are eigenbases with the eigenvalue Xe-
14

---

independent (yet non non-orthose -othoooggonnl in RN*,, the subspace span (fab()) (x) U (1b)) is the K- dimensional eigenspace the Jacobian linear map. From Lemma the rank of the Jacobian matrix J(E) at the fixed point E" is at most K, thereby the rank of the Jacobian linear map in Eq. (15 is most K. Therefore, we have
Ve E span (12)(*) FRELKT u(bb))
JD (e)=
C( + E E" + Jx- (B(e-1))
K-1
Σ" + (x-)es--()) (bb + Σ
t- ExJx (B5(())
X
K-1
E" + x + Σ
Xe- cxb((1)
Ye
k=1
= E(
Lemma 3.1 (Iterated map). For all d [D], the covariance ( = Bigu f obtained recursively by the iterated map C defined by
H
H (8)
is
K-1
(0
FF (H. a)ll? cos (oat)r + o
C
(4)
0,0
k=0
CC(S((-1))
11 ))aaa'
where the expectation is taken over the pre- pre-activations H.d ((e-))))))). g(0) N (a reprexents the scaled positional difference.
o¹)
Proof. For simplicity, introduce Ya,,3,, as follows.
H ((.,))
(e-1) (0(()) aa Fu S ldd
+ (1 (1-66.. 6k.N/22)) V-I=((1))
Ya,k,8.
Since the weights are sampled independently, for different k o k',
E
H (E.k) a.d
H
(H( a.d
H
E
(E,K)
H((kk)
(e(k')]
0.
and
From Eq.
3). we have
K-1
Eg si 11
Σ
Ck
(m + (((
or + ((EAN) ] +E [6]0)2]
d'
2
k=0
15

---

First, we calculate the term k + 0,
N 2 where C₁ = 2.
E
H
+
H.dd)
(HU + H()))
a.d
D
N-1
Σ
Σ
Yakk (o² - (o²- 02) /2D
4 8,8'=0
D
N
Σ
Yar a' (o² - -o²) a2D
f,'=1 8,9=1
D
N-1
YaN.B.Y a (o² + to) D2D
3,9'0
D
N-1
Σ
Σ
Yak.s,.aa (o² + o 1
8.9'
Ya Y, a. + Y F.
3,3
0²
N-1
Σ
F x (e-1)
g ((e-)) Bi,
(It (Fr + Fk.e
=X,, ((e-1)
20
D
2/(21112
(b) (.a))
D
his
cos
where ((-1) is the k-th Fourier mode of the representaion (C-1)
Second, we calculate the terms k = 0, N. where es = and H a.d
(((.))
= H (f,b)
a.
Cl: (H((Ak) 2 a.d + -F (HL + =F
=
Σ
D
Σ
(Ya, v)) 2D
6is
o²
1,8'=1
8,0'=0
0²
{
(0(h) a,a"
D
bu
cos
Since the Fourier modes (f-1))
12 are i.i.d. for each hidden dimension i E [D], we have
K
Egur- Geun---.gw-1,.0rr-1BB Bor0
H (
a. a'dd
(f)
0²
C Be.rrri,guu
x(1 (((11)) 1/9/ cos ((aaa) ) +o²
k=
K
1) cos (o(k) ) + 0z.
Cx E H
NN (
(16)
To obtain more tractable expression in the subsequent proofs of theorems, we express the iterated map without using Fourier modes as follows.
[1
(9)
K-1
Σ
cos (0(k)
( (0(H), a)] of. (17)
Ck
3,0"
8,
16

---

Lemma 3.2 (Exsistance of fixed points), When random DCN defined by Eq. 2 has the fixed point 'q"'''== for the initial parameters (o² of), then random simplified FNO defined by Eq. 3 has
fixed point "" of the form
E' =-= + q C
= =
Proof. Using the properties of cosine functions, the following holds.
K
N N-1
b
2
N2
cos
(o(k)
Cos
(k) B,9'
og.
(18)
k=
B, =0
=5..000
Then, the following holds for all a, [M]
K-1
N-1
[C(2*) )a.a"
Σ cos
Σ
cos
(k()) 3,g'
(O a- (oo))//HHH ( +o))
N2
A=
8,3'=0
When the DCN is applied pointwise to each spatial representation Hs,,, [N], the iterated map of the random DCN (Poole et al. 2016) is given by CDCN (E) == O (0, (, (Hy a)l
oz. Since the covariance Σ is fixed point with respect to the iterative map of the DCN, i.e. CDON(E* = Σ", the following holds.
(C(E") Jan''
Ck cos
(0(4) a,a''
Σ 3,
cos
q (6 + +(1 53 B,P
(err)e" )c*) o¹
N2
=q
Thus, confirm that the covariance E satisfies the definition of fixed point in the map i.e. E' = C(E"). This means that the fixed point for the iterated map E of the DCN also serves as fixed point for the iterated map of the simplified FNO.
Lemma A.I. Let E" be the fixed point of the form in Eq. (5) Suppose that the symmetric perturbation E E NNN, where Eg,3 Eg 3 and gg non-zero for some B' [N], 3 + B, and all other elements are zero. Then, we have
E a E ( + of = + Ea,aXa + +0 ((
(19)
8² O HH N E) (o(He a of =q¹ é- + Ep,0Xn + Eg + O (1 13).
(20)
where Xy's x and Xv. are the constants defined by
Xa = o'HHH [ + (H
Xe o² O -i 18 16' (H),,)) (Hgr
0² 2
Xx =
6 HH
+ H(Ha,) (Hara) a) + 2c (Ha,a)o' (Ha
2016 and Section Proof. Equation 119 is obviously shown by the result of Section 2 in (Poole et al. 3 in (Schoenholzeed et 2016). We prove Eq. (20) with reference the results of (Schoenholz et al. 2016).
a)) + of =²² ( + +
EH a~ +E) 10 (H (HHH)))(HHH)
where f Dd V dze - is the measure for a standard Gaussian distribution, 81 = /q21.
U2 = V (08,8121 + 1
22 = 8² + Eg,g and 6g.8 c* + Eg.a
17

---

We consider the case where c* 1 and c* = separately. Later, we will show that the two results agree with each other. First, we consider the case where approximate ((u1) and ((u2) as follows. For the simplicity, we assume pe))) O (F
Using Taylor expansion. we can
1 Es,s
((u1) = ((1))
+
2 (w)) + O(E
2 Vq
((u2) = o (us) + Eg
c²
E
21
Ø (uz) +
(c 21 + - (u²) ((E,,))))
v1 (ct)
2
where uj V 21 and uz = V + V1 (22))
Thus, we have
d² / D21 D2)))))()))((u))+++++ oo
o² / D2, D2jD2sb(u()((uu) + oî +o² /
Egg
c¹
Dzi
V1 (cm)2 22
)(ui)o'(us)
/q-²
=--0"
EB, (c (1- o²
1 Es
o² DarDee
DinDess 2 /
n (ui))(u)) +
q
The results of the second term are obtained from the transformation of equations 36 to 39 in Appendix 7.2 (Schoenholz et al., 2016).
Eg.9' Vq"
/ Des 21 /c+)2 22 ((ut)d'(uz) -Ea.g o² f (D,, (212)) (ui)
=x=-
Utilizing the identity, (2) (2 we obtain the third term as follows.
1 EEg
2²
D21 Des
2c
(e
2 44-
o 2
/
d² EBg
The last term is calculated as follows.
0² / D2,Dis 2 o
1 Es. 2nd (((()()((z)) = 0² Eg, / DD (co()))
Summing the last two terms, we obtain the term E3,8Xx-
Next, we consider the case where c* = 1. As with the discussion of (Schoenholz et al. 2016). the perturbed correlation is defined by Cg,9' c² Ea.s. where Eg,y and then (0u;) is expanded follows.
((u2) p(uz) + (1/
(V2
Eg.9r Eg.0r 9"
o'(u)) + 2300 + Ea.3 d (w)) 0(
22
q
21
2/
18

---

Thus, we have
o² / D2, D2Dizo((u)((ss) + to
= o² / D2y (isDezdlui)((u)) + +o² / DesDess /2 FBg
220(1]) (ug))
96
-g'c'
EBS
o²
DaiD22
z +o? Dr 220
VT
+ +0² /
DaziDa2
Ea,
21
1 Egg + 02 DaiiDa22 vp
21 (ui) + +
2/4""
(21)
Using the fact that uz ui and un is independent of 22.
o² f
Ea, moll)) = o² /Day/
(ʃ (Dazaa)
Dzi De2
Eg D21
21 ( = o² / (Da
E
0²
(10(2)))')
Vp
Vq
= =o²Ep,8 /
(zi(o((u)) ((ui))"(u)))
0² /
Dz D2i DrID22B 22 230(1))0"
o² 23 / DaIIlui)O""M) (ui) Dz223
(u))) = o² /
E
d2
D2, 2 2
21
Dz1
21 (
2
0²
Dz
D2i(d'( +oujjo""
2
1 Eg d² Dzi D2iDzzs Dz2 v
Eg
21 26 (ui))(uz)
= 0²
D=
21
Vq"
Es,
0²
2
Substituting these facts into Eq. (1I). we obtain
d² / Dzj D21De)9(u1)0(u)) + of
c² - o2 / / of
+ oo
Dzi(o'(ut)2 + o(ut)o"(ui) (ui))
The above result agrees with that obtained by substituting c* = 1 for the result obtained when the case c* < 1.
Lemma A.2. The Jacobian matrix J(E") of the iterated map C defined in Eg. 14 is obtained as follows.
K
cos (o(2)) cos (0952) (6 Xx + 6x + (1 -5..)
[J ) =
N2
k'=0
Furthermore, the rank of the Jacobian matrix J(E") E N2 N is at most K.
19

---

Proof. Let some semi-positive definite matrix E E R.Vx be deviation from the fixed point E".
From Lemma 3.1 and Lemma A.I
we have
[C(2" + E)lara'
K-
N-1
Σ
(B()) a,a'
(o1k))
3,0r) (O
oî.)
ak cos
cos
N
8,0'
K-
N-
N2
Ck cos
(g(k')
a,a!'
cos
(8 B,8
(6 + Eg (1-11111))) (a + EB + Ex.)
9,6'=0
(22)
Note that Eq. (22 is obtained by neglecting higher order terms in Eqs. 19 and 20
By the definition of the fixed point,
K
N Σ 8,81=0
(O(K) 0,0'
(0/4)) q¹ (ó9,00r + (1 )c))
Ck cos
cos
N
By substituting the fact into Eq.
22
we obtain
ICE'* TM
+
E E'lara' laya'
((k')
Σ
(0tk') 8.801 (6g. Ex,axx-) + 6B, + Ex.))
N2
Cke cos
3.8'=U
(23)
N-1
K
(0xg) a.a'
(00) (0g.0r (Xa Xx + Nx, + (1 03 (x.,))
N2
Cp' cos
cos
E
8,3'=0
=/=0
=|J(E>) (1.0),100))
The last equation can be rewritten using the matrix calculation as follows.
C(E* + E) E' mat(J(E") vec(E)).
Since J(E') is the first-order coefficient to the deviaction E, J(E*) is exactly the Jacobian matrix of the itcrated map C.
Furthermore, the Jacobian matrix can be decomposed to two matricies A E Vx and B e x as follows.
K --1
[J (o..)))
Σ
A Bw A(a = cos (olt9)
K=O
B8, = N2 ct cos (0!5))
1
(5 (x
Xx Xx+ ox + +(1 5g.0.
Therefore, the rank of the Jacobian matrix most
Lemma A.3. Let Ex be the fixed point of the form in Eq. (5) and E(k) the perturbation expressed
as
K-1
[.(k] = 12
(k)
(k)
Σ
(0/5)/)
cos
2.9'
Cs
Cs cos
where en denotes the scale of the perturbation, assumed be sufficientll small. Then, we have, for all l E [K]][0}.
C(E¹ + +E()) x E" +
((k,
20

---

Proof. From Eq.
23
we have
K-1
N-1
[C(E" + (E(R)) Ha,a' 2
Ck cos
(g(h) a,ar
(a"
E(R)
xx)
N
k'=0
B=0
K-
N W-1
(0(k)
CDS (0 (9" + - + F -3g
N
OK
008
a,a'
9=0
From the definition of the fixed point, obtain
K
D
cos
(B(k) a.a'
Σ 8,8"=0
cos
(D
9¹ (6 + (1 - 6e,0 )c
N
K==0
By combining the above two results, we have
1 K-1
N
[C(2' E())) - E'Jaa' 2
Σ
(oCk')
E(R
N2
Ck cos
k==0
==0
K-
N
Cke cos
(((*))
a.aa
(o(k')
8,8'
(L (k)
B. X.)))
(k)
By
Since cos (0,)) for all k [K], we have (k)
and E 0. Thus, only the with
(k)
ggg remains.
K K-1
[C(2" + {()))
V"aao' Exxe
N2
1
Ck. cos
(O(k)
a,a!
Σ
cos
(0(6)) 1,87
cos
9,60
(24)
K
(0(h9)
Σ
(o(K')
Σ
(0l.))
Cs
Ck cos
aa,a'
cos
Gs cos
3.31
(25)
Given the orthogonality of the cosine functions, we have
Σ
Cge
(k')
cos
3,00
(K)
=
N2 C
N
(K (otherwise).
(26)
8.0' BA
Utilizing Eq. (26 leads to the following two facts:
K
K
Ckr cos
(e(K') a,
Σ
cos
(((k) 3,9'
cos
3,9'
= N2 cos
(0())
aaaa
Σ
Cpr cos
(o(K')
a,
==
k==
K-
K
K
K
Ck cos
(0(²)
a,a'
cos
9,8'
Σ
Es 008
(0(a)
3,83
N2
cos
(((h))
a.a'
N
Σ
Co
Σ
Ope cos
(0(h))
M=
k==0
s=0fl
A==0
21

---

By substituting these facts into Eq. 24 we obtain
[C'' + E(AA)
"0 Ja,a'
K
1
N² 2º
(0 (,00
( 0,00
~ Xo EK
N2
-
Ck- cos
k'=D
K-1
-1
K-1
1-
Σ
N
COS
(B(K') a.a'
+
Ch--cos
( a.a
$=0
k
K-1
K
(0 a,a''
(o(t?))
Xx=Ek
COs
Co
COs
(4)
This completes the proof.
Lemma A.4. Let E" be the fixed point of the form in Eq. (5)
and E be the perturbation expressed as
K-1
Eg. 61b,,' 33.
N
Xx
Xc
Xa
E
Cs
COs
(0!5)
3,88
Xx
where denotes the scale of the perturbation, assumed to be sufficiently small, 1N 1 is all-one matrix with size N x N.
Then, we have
Σ K-1
Σ
K- Ce
CEE'+E)
==
Xa
1
e=0
(Xx + xo-) E.
N
N
Proof. For simplicity, we introduce x. = N
DK-1
2==0 Ox) as follows.
CE + E
2x=D Co
Xa
+
===0 Cs
)
(Krr + xor) E
Σ.
1
Xx Xo Xo*
Σ s=0
Cs
Xx + Xct
E. (27)
N
Xx
From Eq.
23
we have
CC
E)
"aaa'
K-1
N-1
cos (o(k) a,a''
Σ
Ea.9Xq +
Σ cos (0))2)
(Eg,8Xn + Eo.))
(28)
N2
k==0
A=0
By
Our goal to derive Eq. 27 from Eq.
28
The following results are useful for the computation of
each term of Eq. (28).
K-1
K-1
V3, 8 [N], # Br,
Eg c(1 - x),
Eg,0r =e
1 1-x x
Σ
Co
ca
cos
8,31
==0
8=0
N
(o!52)
=
N2
N
(K' 0)
cos
N
(otherwise)
8=0 B300
N-1
Σ
(g(k')
cos (0())
{
N2
Ck N C
(K' = k) (otherwise)
CK cos
R,3'
3,8'
=
8.9'=0
22

---

Using the above results, the first and second terms of Eq.
28
are calculated as follows.
K
N-1
K
Ck cos
(B(K) 0,a0
Σ
E = c(1 =xx
Ck cos (B
N
N
k'=0
B=
k'=0
K-
N
N
CR oos
(o()))
a,a"
Cos
(0/11) (.90) EB = E(1 x)xx - =(1
x))x
N
Σ
ck oos
(D a,
==0
k'=0
The third term of Eq.
28
obtained follows.
N-1
(o(k')
Σ
cos (o E
cos
N
3=0
Byg
K
N-1
(o(k')
cos (o!k))
N
Cke cos
a.a'
3,8' EXc
k Mm
8m0 By n
(O(K')
(0l)) B,3r cos (09,5) EXc
x
C2
N
Ck cos
a..0
C2 cos
(o(k))
=ee'- Exe
CK COS
N
k'=0
K-1
K-
K-1
Σ Ck cos
(olh'))
CK cos (Ql.)
X
C2
xxxxx
N
$=0
k==0
CXe Xx
Σ
C2
Ck COS
(o(t2)
c(1
x)
Σ
cke cos
(B(k')
N
k==0
k-=0
By substituting these facts into Eq.
28
we have
(CEE* E) Ed aa."
K-1
K-1
2 =e((1 - x)xx + xxx) ExxEx
Σ
Co
Cp. cos
(oct')
s=1l
===0
1
Ck' (ork')
-c11 x) N (xx Xc- Xx
k==0
Exmx E0 co)
K K
-1 K-1
- (11=x)xx + x
1-
Σ
Cs
Cee cos
(o(k') a,a')
0=0
Finally, we obtain Eq. (7
B Proof of Theorem 3.5
Theorem 3.5 (Trainability). Let E(0) NxN be the gradient covariance with respect to loss
L. e.g. mean squared error, at the E-th layer. Suppose that the gradient covariance at the L-th layer decomposed as S(I.) - K-1, & cos (0tt) èē, where Et is the coefficient of each basis and belongs to the orthogonal complements of sssssan span[(cos] (o(t)) ) ) Then, the gradient covariance
23

---

at the E-th layer is obtained by
K Σ k=0
{
L-t Ex
cos (oltd)
a,a'
Xo
Proof. Recall the definition of the gradient covariance. We first demonstrate the iterated map of the gradient covariance, starting from this definition.
500) 0,a' B
(6) (8) = [as
Oc
N-1
Σ
Oc
(H!!
8=0
The Jacobian
is calculated as follows. a a.
(2+1) OH
3,
0H
K-1
Cl
(0 J + +V=I=ftt1a) + Faa (0 (0/(t1,)) - of (u
--T
Fkio
2
k=
K-1
Ck
(e)
COS (08))
1,) sin (0%/)) ø1 o((H
N
2
k=1
am aH
The covariance of Jacobian
is follows.
(H(( (H+1)))
Ba
Ey ELL
OH
B ((N
N D
Ck coos
(0(4) 9,00 cos (0%% sin (0(k) sin (0%).r)
2²
EH HH ((..00)))))
(H
Σ Ck cos (0 (0f
(022
N
Then, the covariance of the gradient is given by the following recurrence relation for a, al [N],
(0
(8) (9)
Egii gaL
0,a'
N-1
OH (1-1)
aH!!+1)
(F+1) (e+1)
B
I3.i 93
0 oH
8.9'=0
N
D
aH (e-1)
aH (1+1)
Σ
Σ
E
111 1411.g
93
(e+1) (e+1)
E
04 (
OH!0
8,9' 13'=0i,' i
N-1
K-1
(F+1) B.8'
EH
Σ
CK cos
(0
(0(A)
As with (Schoenholz et al. 2016). we approximate (+++) se Σ' since the number of layer is assumed be sufficiently large. Then, the linear iterated map of the gradient covariance ((+1)) (0) is given as follows.
N-1
K-1
500 ".a
Σ
((++)) "B,0'
Xe:
Σ
ck cos (k( (03-8r,a-)) (0
(29)
N²
3,0'=0
24

---

The rank of the linear iterated map E(++1) 1) H E is less than K since the matrix representation of the linear map can be decomposed into two matrices A N=N and B E IR K follows.
1 K-1
5 a.a
Σ Cs exp
( (-V-10
3,0' ) exp (V-10(g) a',a Xe 333"
5(et1)
A(a.a)), Bu
E
N2
0.3"
A(aaa),
N2 Che exp (V-10 ) Xo BR.c,,3) = exp (-V-100h)
Next, we show that the subspace span (fcos eigenspace with cigenvalue Xe. of the linear iterated map ((++))
(00) aaa'
RNNN is the K-dimensional
1A=0
E
By substituting -/
Ek=o
Xet
cos I.-(e+))
Ex (e()) B,8
into Eq. 29 obtain
N
K-1
K
50
Σ k=
(olt)
Σ
(ol4)
-- (3-P(cc-0) )
0,
N2
Xe
Es cos
0,000)
cos
9,6'=0
k'=0
From the orthogonality of the cosine and sine function, we obtain
1
Σ
(005) (0- (02-8,a-0')) cos
Cp cos
N2
8,9'
N-1
Σ
Σ
Ci. cos (0f/)) cos (0452) cos (A(A)
N
(30)
3
(o(k3) 9,00
sin ( cos 3,8'
+ sin
cos 0
(N)
Hence, we have
K-1
(5()) a.a'
Σ
cos (o(t)
Xe
k=0
This completes the proof.
C Discussions on the similarity of DCN and CNN
While CNNs perform local convolutions in the spatial domain, FNOs execute convolutions in the frequency domain, thereby achieving global convolutions in the spatial domain. Our theory for the FNO and the mean-field theory for CNNs (Xiao et al. 2018) share common focus the correlation dynamics ()) [1) (()) of the spatial representations ())), H HL) E Nx The
iterated maps for CNNs and FNOs are obtained as follows.
(++1))
Σ En (( (0(Ho + +o =: CCCNN(E
0,00
2r
+1
Beker
K-1
(1)
a,a" 0²
Σ
Og HH
[[ (M) [ cos (0(h) a.
+ + Coxo 1 (1)),
where 2r + 1 is the number of filter width and ker = {B E | S r} is set indices referring to the elements of the filter.
When we consider the FNO without mode truncation (K N/2+1 1), the propagation of the diagonal components (a, a) for any E [N] is equivalent to a CNN with filter size N performing global
25

---

convolution.
(())
a,a
o²
cn En ~
[|L (H (Hll
cos
(ock) a,a
to
k=0
r² N-1
H a(0.E1)) [1o(Ha()))
tob
E -(0,000) (.) (. o
The first equality follows from Perseval's equality and the second from the periodic boundary condition. In contrast, the propagation of the off-diagonal components differs between CNN and FNO, and the iterated map is different even in the presence of mode truncation.
Regarding the fixed point, Xiao al. (2018) demonstrated that any fixed point for the iterated map of the DCN is also a fixed point for that of the CNN. Consequently, Lemma A. indicates that the fixed points for CNN and FNO are consistent.
The behaviour of the iterated map around fixed points reflects the nature of each architecture. CNN possess diagonal eigenspaces associated with eigenvalues Xy and non-diagonal eigenspaces associated with eigenvalues Xe*. FNOs without mode truncation exhibit similarity, possessing eigenspaces Xa* for zero-frequency and eigenspaces Xe* for k-frequencies with diagonal components removed.
Finally, the iterated map during the backpropagation for CNN and FNO given by 5,()) 0,00"
Σ BEker
U3 U3
S(++1)
B (H, (d
a-B,a'
(0 0,00'
o²
(1+ B,8
(0(H (H
Σ ck cos (0 (022)a-)))
N2
=
where Va is the variance weight parameter dependent of the filter position 3, i.c. wi NO0, ((g)))) Esekor Ug
Using the approximation ((++)) se E' as with (Schoenholz et al. 2016), the backpropagation the diagonal components (a..0) of the FNO without mode truncation is equivalent that of the global CNN (see Eq. (2.16) in (Xiao et al.
2018) with Uy N. as shown below:
N-1
( 0,0
N2
Σ
(e+1)
3,9
Xe
CK cos (0!0) (09-8ra-a))
$ = Xe
((+1) 3,9
3,8'=0
NN
This equivalence suggests that the edge of chaos initialization (e.g., He initialization) is also valid for the FNO since the problems of gradient vanishing and explosion are determined by the diagonal components of E.
D Details of Experimental Setup
In this section, we summarize the detailed setup of the all experiments, including the experiments in Section[4]
D.1 Datasets
D.1.1 Advection equation
We used the advection equation data published by Takamoto al. (022).. The advection equation for the function u(x, 4) E (2, 1) x (0, 2]: TR) is given by
0,u(x, t) + 30,, (u(2, 2)/2) 0, u(r, 0) = uo(x),
26

---

where uo E ((0, R) is the initial condition and / E R is an advection speed set to 2.0. The exact solution is given as u(x,t) uo( 5t) for any initial condition un-
Only periodic boundary conditions were used in this dataset. The initial conditions are the super- position of sinusoidal wave given by
Eman
Σ
sin (ke + 0.),
(31)
where h;; = 2T Ey= R Ly are wave numbers whose 74. are integer numbers randomly chosen in [1, kmax]. Li = is the calculation domain size, N = 2 is the number of wave to be added, and kimax 8 is the maximum wave number. The amplitude A; is uniformly chosen in [0, 1], and the phase Φ; is the randomly chosen in (0, 2T). The 2nd-order temporal and spatial upwind finite difference scheme was used for generating the data. Settings are described in Appendix D of (Takamoto et al. 2022).
D.1.2 Burgers' equation
We used the Burgers' equation data published by Takamoto et al. (2022. The Burgers' equation for the function u(x,t) t) E ()) x (0, 2|: is given by
0zu(x, t) + az (u² (z, 2)) = vd u(x t),
u(x, 0) un ux()),
where Un E ((0, 1);R) RR) is the initial condition and the diffusion coefficient set to 4.0. The periodic boundary conditions and Equation (31) are used as the initial conditions. The 2nd-order temporal and spatial upwind finite difference scheme is used for generating the data. Settings are described in Appendix D of (Takamoto et al. 2022).
D.1.3 Darcy Flow equation
We used the data of 2D Darcy Flow equation on a regular grid published by Li et al. (2020a). The Darcy Flow equation for the function u E H) ((0, RR); with a Dirichlet boundary is given by
(a(x))u(x)) = f(x),
x E 1)?,
uxx)= = 0,
200,1)²
where E ((0, is the diffusion coefficient and f (,, 1)2; IR) is the forcing function.
The coefficients was generated by measure μ vo, (-A 97) 22) using the Laplacian with zero Neumann boundary and the binary point-wise mapping vixx) 12 (z 2 0). 3 (x 0). The forcing function is fixed fxx
Our aim is to predict the operator mapping the diffusion coefficient to the solution a - u. The solution function u was generated by using the second-order
finite difference scheme on a 421 x 421 grid. Settings are described in Appendix A.3.2 of LL et 2020ch.
D.1.4 Incompressible Navier-Stokes equation
We used the 2D NS equation on the unit torus defined by
a,l(I, 1) + u(x, t) Vw(x,
V vuu(z,t) 0, 0)=
where t(, (() is the vorticity, wp H", R)2 is the initial vorticity, u(I,t) E C(O0,T] 1);; is the velocity field for any r > 0, v E R,+ (2) is the external forcing function defined by f(x) 0.1 (sin(2 + x2) + c00((2 + +x2)). The initial vorticity wp was generated by wp ~ l where (3 + 49.) with periodic boundary conditions. The viscosity was set to 1e-3, 1e-4, or le-5. Our aim is to predict the operator that maps a solution u up to time 10 to solution up to some later time T > 10. The data was generated by the pseudo-spectral Crak-Nicholson second-order method on 64 x 64 grid. For the data with v 1e-4, the time resolution was also downsampled by half. Settings are described in Section 5.3 of (Li ct al. 2020c).
is the viscosity, and
27

---

Table 2: Training settings
PDE
architecture
batch size
initial Ir
max max_eppohs
advection
Simplified FNO with Tanh
40
5.0 10-3
200
advection
Simplified FNO with ReI.U
40
5.0 x 10-3
200
Burgers'
Simplified FNO with Tanh
40
5.0 x 10-5
200
Burgers'
Simplified FNO with ReLU
40
1.0 10-3
200
Darcy flow
2D FNO with ReLU
20
1.0 x 10-3
500
Navier-Stokes (v = le-3)
2D FNO with ReL.U
20
x.0 x 10-1
500
Navier- (v le-4)
2D FNO with ReLU
50
2.5 10-3
400
Navier-Stokes (v Le-5)
2D FNO with ReLU
20
2.5 x 10
500
D.2
Training settings
Our detailed training settings of the experiments in Section are provided Table Our experimental environment consists of an Intel Xeon Plantinum 8360Y (36-core) CPU and single NVIDIA A A100 GPU. Most of our code for experiments are based on the code PDEBench (https://ga thub http//eebnncc/PDEB...
(Takamoto et al. 2022). The only modifications the model are to multiply the outputs (variable out_ft in code of class FNOId and FNO2d) corresponding to mode 0, N/2 by /2 and to initialize the weights by Gaussian distribution, as described Section 3.1
E
Detailed Experimental Analysis
E.1 Analysis of Training Loss
Figure 6a shows the training loss for each epoch of the 32-layer FNOs with parameters o² e (0.1, 0.5, 1.0, 2.0, 3.0, 4.0} on the NS dataset with 1e-3. When the initial parameter o² is too small, the training loss is not well reduced due to gradient vanishing. On the other hand, when the initial parameter o² proposed edge of chaos initialization smoothly reduces the training loss in the initial epoch and enables stable training.
large, the initial training loss blows up due to gradient exploding. The
F.2 Analysis of Test Performance
The nMSE of the FNOs on test datasets for six distinct PDEs is presented in Tables 3 and Results are shown only for the FNOs with initial parameters where training was successful in many cases. For the NS equation with viscosity values of 1e-3, 1c-4, where sufficient data is available, Table shows that best performance is achieved with 8 or 16 layers. This suggests that while shallow FNOs are currently prevalent, deep FNOs could be advantageous in certain tasks, underscoring the significance of our analysis the bias in deep FNOs. Conversely, for other equations, the 4-layer FNO performs best and deeper FNOs result in drop in performance even with the edge of chaos initialization. We will discuss this test performance deterioration in detail.
The over- over-fitting phenomenon is observed in the Darcy Flow and NS equation datasets with 1e-5. where only limited training data available. The training loss for each epoch on the equation depicted in Fig. (bb] As demonstrated in Fig. [b]] the 16 and 32-layer FNOs yield lower training loss than the 4-layer FNO, but exhibit poorer performance on the test dataset as shown in Table 4 These results suggest over-fitting the training data, necessitating either abundant training data or appropriate regularization.
Conversely, the under-fitting phenomenon is apparent in the ID advection and Burgers' equation datasets in T The training loss of the FNO with ReLU activation for each epoch on the Burgers' equation is presented in Fig. 6c] Figure 6c indicates that the larger the number of layers, the higher training loss in the final epoch, and worse the test performance. This under-fitting to training data could be attributed to the escalating complexity of the loss landscape as the layer count increases, a known issue for DCN and CNN (Li et al.] 2018). This may be due to the emergence of local minima corresponding to operators that generate too complex functions, preventing the
28

---

100
10
1011
(a) NS eq. (Ic-3)
(b) NS eq. (Ic-5)
(c) Burgers' eq.
Figure 6: Training Loss Curve. (a): training loss curve of the 32-layer original FNOs with varying
initial parameters o² {0.1, 0.5, 2.0, 3.0, 4.0). on the NS equation with v = 1e-3. (b): training loss curve of the original FNOs with an initial parameter o² 2.0 with varying number of layers L E {4, 8, 16, 32} on the NS equation with y = le-5. (c): training loss curve of the simplified FNOs with ReLU activation and the initial parameter
2.0 with varying number of layers
L E {4, 8, 16, 32) on the Burgers' equation,
Table Test performance measured by nMSE ID simplified FNO on ID PDEs
Advection
Burgers'
aMSE
Layers
Tanh
ReLU
Tanh
ReLU
0.013
0.015
0.0055
0.00088
0.013
0.015
0.0069
0.0012
0² = 1.0
0.013
0.015
0.0068
0.0016
32
0.016
0.018
0.0071
0.0041
0.013
0.017
0.0036
0.00098
0.012
0.018
0.0034
0.0011
0 = 2.0
16
0.012
0.020
0.0050
0.0016
0.014
0.024
0.0062
0.0027
0.013
0.019
0.0044
0.00096
0.014
0.022
0.0042
0.0012
0² =
3.0
16
0.020
0.032
0.0060
0.0016
32
0.053
0.059
0.0093
0.0045
attainment parameters that achieve global minima. This issue could be mitigated by introducing more suitable regularization, appropriate optimizer, skip connection (Tran et al. 2022).
Our theory and experiments suggest that the training of deep FNOs has suffered from problems including gradient vanishing and exploding due to improper initialization, over-fitting caused by insufficient training data, and mnder-fitting caused by loss landscapes with strong non-convexity. While our edge of chaos initialization prevents the gradient vanishing and exploding, techniques to solve over-fitting and under-fitting problems are still needed in practice.
Visualization of Forward Propagation
We visualized the behavior parameters o² E {0.1, 1.0, 2.0, 4.0} and activation functions. The FNO, with width of D = 1024, was used and the input was sampled from the standard normal distribution with a spatial size of N 32. The results of the FNO with Tanh activation. both with and without mode truncation, shown in Figs. 7 and [ and Figs. and 10 respectiivyl Similarly, the results of the FNO with ReLU activation, both with and without mode truncation, and [eepppeciv In the ordered phase, all figures illustrate convergence to the fixed point Σ' where c" = I, with the rate of convergence increasing the parameter a² decreases. In the chaotic phase, the activation function dictates the covariance behavior. Without mode truncation, the covariance behavior of the FNO identical those of the DCN; otherwise non-uniform, FNO-specific periodic covariance is exhibited.
simplified FNO's covariance matrix [ with varying initialization
displayed in Figs. IT and 12 and Figs. B
29

---

Table 4: Test performance measured by nMSE of 2D original FNO with ReLU activation on Darcy Flow and NS equation.
Navier-Stokes le-4
aMSE
Layers
Darcy Flow
== 1e-3
V/
V 1e-5
0.025
0.0063
0.18
0.10
8
0.028
0.0047
0.14
0.094
1.0
16
0.035
0.0048
0.12
0.12
32
0.56
0.0057
0.13
0.16
4
0.029
0.0075
0.18
0.10
0.036
0.0057
0.14
0.11
0² 2.0
16
0.041
0.0057
0.12
0.11
32
0.042
0.0072
0.13
0.18
0.033
0.0089
0.18
0.10
8
0.040
0.0080
0.14
0.11
3.0
16
0.052
0.0098
0.13
0.13
0.16
0.028
0.14
1.00
e1
Horld
11
--20
-0.35
mabil poitee Ewo
4050 24800
rerld peic
-0.30
Twrl
11
-1175
gudd p4
uatal podiee
cotid
-100
Figure 7: Visualization of the correlation (1)
B.8r/ VV
(2) ( 4
3.. for the simplified FNO with Tanh
activation and no mode truncation.
30

---

3.0
175
sodid palle
qpd.dppaii 2-
Sadal pedce 11
1p4sa publoc
1.50
[-1
e1-13
125
CAM pmdi
volid pobloe
qua.dppeii 2-
1p4ad pedie
speia poblic
E
0.75
splld polloe
qutilpetti
speld publin
0.50
2=
0.25
queld
3
pubko
0.00
Figure 8: Visualization of the covariance ( for the simplified FNO with Tanh activation and no mode truncation.
quiid
e-10
vadial mulu
vodal podce
colld petice
5.00
Ewo
<- 11
t-
0.25
vdaal podiow
Vpile prlas
-1.53
Feo
Eoc
Zm 16
30
-075
(0)
2,,0 / /> (1) (() 3,8 33 for the simplified FNO with Tanh
Figure 9: Visualization of the correlation activation and the Fourier mode
3
31

---

3.0
175
sodid palle
qpd.dppaii 2-
Sadal pedce 11
1p4sa publoc
1.50
[-1
e1-13
125
Sodid poblo
qua-dppeiii 2-.8
1p4ad pedie
speia pillic
E
0.75
speld 20
qutilpetti
spel purkin
0.50
2=c
0.25
surld
20
3
pubkos
0.00
Figure 10: Visualization of the covariance ( for the simplified FNO with Tanh activation and the Fourier mode K = 5.
quiid
e-10
vadial mulu
vodal podco
5.00
Ewg
11
t-
-025
vddal podiew
Vpile prilas
-1.53
Fel
Eoc
Zu 15
-075
49/00
Figure 11: Visualization of the correlation activation and no mode truncation.
3,31 / /
23 3,00 3 for the simplified ENO with ReLU
(1)
32

---

X-4
C-11
C-41
2-4
C-36
C--22
vodad pmiss
(a) Initial parameters d² {0.1, 1.0, 2.0}.
4-5
(b) Initial parameter a² = 4.0
Figure 12: Visualization of the covariance (4) for the simplified FNO with ReLU activation and no mode truncation.
vaa podie
spale polass
0.50
Fe.O
6-1
Eood
6-6
Z= 16
3.25
2.00
-100
Figure 13: Visualization of the correlation 2. activation and the Fourier mode K 5.
(0) g
for the simplified ENO with ReLU
33

---

X-4
paido
poc
soeld purko
o1=30
(a) Initial parameters o² 1.0,
20.4
(b) Initial parameter 0² 4.0
Figure 14: Visualization of the covariance ( for the simplified FNO with ReLU activation and the Fourier mode K 5.
34

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The contribution and scope of the paper are described in the abstract and introduction.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We mention the limitation of this paper in Section
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
35

---

Justification:
The full set of assumptions and a complete proof are provided in Section
and Appendices Andd
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: All the information needed to reproduce the main experimental results are provided in Section 4 and Appendix D
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
36

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes)
Justification:
We use existing public datasets and cite them appropriately.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cr
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ praips.cc/publiccguu..ee CoodeSumeeieeeeee ssi iessonPolll molyicyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All training and test details are provided in Section 4 Appendix D Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
Justification: In the experimental results depicted in Figs. and determined that the variation in the reported values has minimal impact. In fact, the inclusion of error bars seems obstruct intuitive understanding through visualization.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
37

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
show in tables or
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: All sufficient information on the computer resources are provided Section and Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics htttp:://neun: htt.ee/iihiiee//rroooooooooooooa
Answer: [Yes]
Justification:
We follow the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Our paper is foundational research and not tied to particular applications.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
38

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We have stated the license and copyright of the data used.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in package should be provided. For popular datasets, paperswithcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
39

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
The content of our public code is well described.
Guidelines:
The answer NA means that the paper does not release new assets. Researchers should communicate the details the dataset/code/moodl as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
40


Reviews:
[{"review_id": "weilmP9CMv", "reviewer": "Reviewer_Jbgf", "summary": "This work explores the expressivity and trainability of Fourier Neural Operators (FNOs) by establishing a mean-field theory and analyzing their behavior from an edge-of-chaos perspective.", "questions": "1. Can the derivation in this paper for FNOs be generalized to more general neural operators?\n2. Is it possible to identify (approximately) how the initialization depends on different datasets?", "limitations": "Yes.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The study examines the ordered-chaos phase transition of the network based on weight distribution, highlighting unique characteristics induced by mode truncation and similarities to densely connected networks. It also identifies a connection between expressivity and trainability, noting that ordered and chaotic phases correspond to regions of vanishing and exploding gradients, respectively. This insight provides a practical prerequisite for the stable training of FNOs. The experimental results support the theoretical findings.", "weaknesses": "1. The summary of related work can be further improved, particularly by including recent neural operator architectures from recent conferences or journals.\n2. While Figure 1.(c) confirms that the original setup for initialization is reasonable, our experience suggests that the performance of FNOs in both training and test accuracy is highly dependent on the initialization variances for different hyperparameters and tasks. The derivation and results in this paper do not provide sufficient insights into this issue.\n3. The code is not available.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> Q1. The summary of related work can be improved, particularly by including recent neural operator architectures\n\nWe have revised Section 2.1 as follows. (Only newly added references are cited in numbered form.)\n\n**line 48:** The FNO (Li et al., 2021) is one of the well-established methods for solving PDEs across many scientific problems (Hwang et al., 2022; Jiang et al., 2021; Pathak et al., 2022; [1, 2]).\n\n**lines 58 - 71:** Several FNO variants have been developed to address specific challenges, such as geo-FNO (Li et al., 2022a) for irregular regions and group equivariant FNO (G-FNO) (Helwig et al., 2023), which maintains equivariance to rotation, reflection, and translation. U-NO[3] and U-FNO[4] integrate FNO with U-Net for multiscale modeling. Additionally, WNO[5] utilizes wavelet bases, while CFNO[6] enhances the use of geometric relations between different fields and field components through Clifford algebras. Adaptive FNO[7, 8] and F-FNO (Tran et al., 2022) have improved computational and memory efficiency through incremental learning and architectural modifications. Other approaches for improving performance include methods with increasing physical inductive bias[9], data augmentation[10], and a variance-preserving weight initialization scheme (Poli et al., 2022). While numerous new models and learning methods have been proposed, relatively little research has been conducted to understand the intrinsic nature of these methods. Issues such as spectral bias[7] and training instability (Tran et al., 2022) have been reported. \n\n[1] Yang, Y., et al. Seismic wave propagation and inversion with neural operators.\n\n[2] Wen, G., et al. Accelerating carbon capture and storage modeling using fourier neural operators.\n\n[3] Rahman, Md Ashiqur, et al. U-no: U-shaped neural operators\n\n[4] Wen, Gege, et al. U-FNO—An enhanced Fourier neural operator-based deep-learning model for multiphase flow\n\n[5] Tripura, Tapas, et al. Wavelet neural operator: a neural operator for parametric partial differential equations\n\n[6] Brandstetter, Johannes, et al., Clifford Neural Layers for PDE Modeling\n\n[7] Zhao, Jiawei, et al. Incremental spectral learning in Fourier neural operator.\n\n[8] Guibas, John, et al. Adaptive fourier neural operators: Efficient token mixers for transformers.\n\n[9] Li, Zongyi, et al. Physics-informed neural operator for learning partial differential equations.\n\n[10] Brandstetter, Johannes, et al. Lie point symmetry data augmentation for neural PDE solvers.\n\n---\n\n> Q2. The code is not available. \n\nThe paper does not propose any new models or learning methods and most of the code used in the experiments is based on the code of PDEBench. The only modifications to the model are to multiply the outputs (variable ‘out_ft’ in code of class FNO1d and FNO2d) corresponding to mode k=0, N/2 by $\\sqrt{2}$ and to initialize the weights by Gaussian distribution, as described in Section 3.1. \nWe added the same information as above to the detailed experimental setup description in the paper.\n\n---\n\n> Q3. Can the derivation in this paper for FNOs be generalized to more general neural operators? \n\nOur theory may be applicable to several variants of the FNO. Examples include the following:\n* geo-FNO: For geo-FNO, the influence of deformation is limited to the beginning and the end of the architecture, with the majority of the architecture consisting of the original FNO. Thus, our theory can be directly applied.\n* G-FNO: For G-FNO, similar results can be obtained if the initial distribution of the transformed weight $\\hat{L}_{s_g} R$ by elements $s_g$ from the stabilizer $S_g$ satisfies the assumptions of our theory.\n* WNO: Since the wavelet transform is a linear operator, the same analysis can potentially be applied to the WNO.\n\nConversely, applying our theory directly to the architectures of U-FNO and F-FNO is challenging due to the presence of skip-connections. In architectures with skip-connections, the degeneration of input information in the ordered phase does not occur, and chaotic behavior is expected to be prominent. Therefore, it is necessary to analyze the behavior of the input in the chaotic phase in detail, rather than focusing on the discussion of the fixed point. It may be possible to extend our theory by incorporating mean-field theory for DCNs with skip-connections (Yang & Schoenholz, 2017) alongside our theoretical framework. Once the discussion on skip-connections is resolved, similar arguments can be made as follows: \n* U-FNO: For U-FNO, it is considered possible to discuss by integrating the results related to CNN, DCN, and simplified FNO, as discussed in Section 3.3 on the original FNO. Note that skip-connections are included in U-Net.\n* F-FNO: The factorized Fourier convolution module (Eq. 8) can be discussed similarly to the conventional FNO. However, since the output of the inverse FFT is summed in the dimensional direction, it seems appropriate to consider the weight scale as 1/D.\n\n---\n\n> Q4. Is it possible to identify (approximately) how the initialization depends on different datasets?\n\nOur theory and experiments provide insight that, regardless of the dataset, loss function, optimizer, or scheduler, initialization around $\\sigma^2 = 2.0$ is one of the necessary conditions (not sufficient condition) for stable training, especially for deep FNOs. \n\nAt the same time, we also agree with your experience that the performance is dependent on hyperparameters and datasets (tasks). As discussed in Appendix D, underfitting caused by a loss landscape with strong non-convexity and overfitting caused by an inappropriate setting of model capacity relative to the amount of data are both hyperparameter- and data-dependent. Our findings suggest that the acceptable initial variance of the weights is already limited by the number of layers, so these problems need to be addressed by optimisers, schedulers, devising loss functions, introducing regularization, etc. Thank you for your valuable comments.", "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": "TkTzyzte8O", "reviewer": "Reviewer_xhAf", "summary": "This paper presents a comprehensive theoretical analysis of the Fourier Neural Operator (FNO) using mean-field theory. The main contributions to me at least are:\n\n1. Establishing a mean-field theory for FNO, analyzing its behavior from an edge of chaos perspective.\n2. Examining the ordered-chaos phase transition of random FNO based on weight distribution.\n3. Identifying a connection between expressivity and trainability, showing that ordered and chaotic phases correspond to regions of vanishing and exploding gradients respectively.\n4. Providing practical initialization requirements for stable FNO training.", "questions": "1. I know this paper tackled FNO's, but I am curious if this could easily be adapted to other variants of FNOs like U-FNO, F-FNO etc\n2. The paper doesn't provide a detailed analysis of how different PDE characteristics (e.g., stiffness, non-linearity, multi-scale behavior) might affect the expressivity and trainability of FNOs. Could you explain more on this point?", "limitations": "The authors have described their limitations!", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "1. Theoretical Depth: The paper provides a rigorous mathematical analysis of FNOs, which has been lacking in the literature. The authors adapt and extend existing mean-field theory techniques to the specific structure of FNOs.\n2. Novel Insights: The paper reveals FNO-specific characteristics in the phase transition, induced by mode truncation, while also showing similarities with densely connected networks. This bridges the gap between FNOs and more traditional neural network architectures. This was very very insighting to me!\n3. Practical Implications: The theoretical analysis leads to practical guidelines for initializing FNOs, addressing a known issue of training instability in deep FNOs. This demonstrates the real-world value of the theoretical work.\n4. Comprehensive Experiments: The authors validate their theoretical findings through extensive experiments on various PDE-solving tasks, showing the practical relevance of their analysis.\n5. Clarity: Despite the technical nature of the content, the paper is well-structured and clearly written, making it accessible to a broader audience in the machine learning community.", "weaknesses": "1. There are no specific weaknesses that I see however, while the experiments are comprehensive, the authors don't fully explore very deep FNOs or a wide range of hyperparameters (in particular sweeping over the number of modes). Additional experiments on larger-scale problems or more complex PDEs could further strengthen the empirical validation.\n2. Insufficient analysis of the impact of different Fourier modes: The paper doesn't provide a detailed analysis of how the number and distribution of Fourier modes affect the expressivity and trainability. This is crucial as the choice of modes is a key hyperparameter in FNOs.\n3. Insufficient exploration of error propagation: For time-dependent PDEs, there's limited analysis of how errors propagate through time steps and how this relates to the expressivity and trainability of FNOs.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> Q1. The authors don't fully explore very deep FNOs or a wide range of hyperparameters (in particular sweeping over the number of modes). Additional experiments on larger-scale problems or more complex PDEs could further strengthen the empirical validation.\n\nThank you for your valuable suggestions. We are currently conducting additional experiments. If the results are ready in time for the discussion period, we will report them in the thread.\n\n---\n\n> Q2. I know this paper tackled FNO's, but I am curious if this could easily be adapted to other variants of FNOs like U-FNO, F-FNO etc\n\nOur theory may be applicable to several variants of the FNO. Examples include the following:\n* geo-FNO: For geo-FNO, the influence of deformation is limited to the beginning and the end of the architecture, with the majority of the architecture consisting of the original FNO. Thus, our theory can be directly applied.\n* G-FNO: For G-FNO, similar results can be obtained if the initial distribution of the transformed weight $\\hat{L}_{s_g} R$ by elements $s_g$ from the stabilizer $S_g$ satisfies the assumptions of our theory.\n* WNO: Since the wavelet transform is a linear operator, the same analysis can potentially be applied to the WNO.\n\nConversely, applying our theory directly to the architectures of U-FNO and F-FNO is challenging due to the presence of skip-connections. In architectures with skip-connections, the degeneration of input information in the ordered phase does not occur, and chaotic behavior is expected to be prominent. Therefore, it is necessary to analyze the behavior of the input in the chaotic phase in detail, rather than focusing on the discussion of the fixed point. It may be possible to extend our theory by incorporating mean-field theory for DCNs with skip-connections (Yang & Schoenholz, 2017) alongside our theoretical framework. Once the discussion on skip-connections is resolved, similar arguments can be made as follows: \n* U-FNO: For U-FNO, it is considered possible to discuss by integrating the results related to CNN, DCN, and simplified FNO, as discussed in Section 3.3 on the original FNO. Note that skip-connections are included in U-Net.\n* F-FNO: The factorized Fourier convolution module (Eq. 8) can be discussed similarly to the conventional FNO. However, since the output of the inverse FFT is summed in the dimensional direction, it seems appropriate to consider the weight scale as 1/D.\n\n(Yang & Schoenholz, 2017) Yang, Ge, and Samuel Schoenholz. \"Mean field residual networks: On the edge of chaos.\" Advances in neural information processing systems 30 (2017).\n\n---\n\n> Q3. The paper doesn't provide a detailed analysis of how different PDE characteristics (e.g., stiffness, non-linearity, multi-scale behavior) might affect the expressivity and trainability of FNOs. Could you explain more on this point?\n\nIn terms of expressivity, it depends on the architecture of the FNO and not on the characteristics of the dataset (PDE). Regarding trainability, our theory provides a universal insight, independent of the dataset, loss function, optimizer, and scheduler. While it does not offer sufficient conditions for stable training, it does provide a necessary condition for achieving stability in training. As discussed in Appendix D, the characteristics of the dataset may, for example, lead to strong non-convexity problems in the loss landscape, but this is not explained by our theory. Our findings suggest that the acceptable initial variance of the weights is already limited by the number of layers, so these problems need to be addressed by optimizers, schedulers, devising loss functions, introducing regularization, etc.", "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": "0trxv1BoD9", "reviewer": "Reviewer_ePA1", "summary": "This paper leverages mean-field theory, developed for neural networks, to study the expressivity and trainability of random Fourier neural operators (FNOs). It provides practical insights for the appropriate initialization of FNOs.", "questions": "- Does the fixed point of the iterated map always exist?\n- The representation of Fig.1(a) is unclear. Does it depict the gradient norm averaged over the training epochs? I assume it shows the values only at initialization, averaged over multiple seeds. Additionally, is the standard deviation relatively small? Clarification on these points would be helpful.", "limitations": "I understand that error bars are omitted from the figures for the sake of clarity. However, including standard deviations in the tables would enhance the presentation of the results. Additionally, it is concerning that the results are averaged over only three seeds, which is insufficient for robust statistical analysis.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The authors address a significant question in the neural operators literature: why does the performance of Fourier neural operators degrade as the number of layers increases? The authors take a first step toward answering this question by adopting an original theoretical perspective.", "weaknesses": "- The main body of the paper is not self-contained. Additionally, the technical section can be difficult to follow due to the need for back-and-forth referencing between multiple sections and references. Overall, the paper could benefit greatly from improved clarity. To enhance readability, I recommend reorganizing some parts of the text.\n  - For example, the covariance matrix $\\Sigma_{\\alpha,\\alpha'}^{(l)}$ should be introduced in section 3.2 where it is first examined.\n  - At first glance, one might question if the iterated map in equation (6) is a definition. However, its origin becomes clear upon reviewing the supplementary material.\n  - Lemma A.1 should be included in the main text, as it is a central result.\n  - Some parts, such as the detailed description of the datasets, can be moved to the appendix to make space.\n- I recommend providing further motivation for all concepts and referencing previous works where appropriate. For instance, the importance of studying the fixed point of the iterated map should be emphasized. Most importantly, please refer to the supplementary material for proofs when available.\n- The proposed study focuses solely on the expressivity of random FNO at initialization, a rather limited scenario. This limitation should be more explicitly emphasized, particularly in the title. I am not convinced that the findings offer practical insights into the stable training of FNO, despite the preliminary numerical experiments provided. Additionally, the term _trainability_, as used by the authors, only pertains to the vanishing or exploding gradients of random FNO.\n- One might question the impact of the optimizer, the initial learning rate, and its scheduler on the results. An ablation study addressing these factors would have been beneficial.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Due to questions similar to those from other reviewers, we have addressed these in the global rebuttal. Please refer to our global rebuttal. We appreciate your attention and review.\n\n> Q1. The covariance matrix should be introduced in section 3.2. \n\nThe covariance $\\Sigma_{\\alpha, \\alpha’}^{(\\ell)}$ is defined in Section 3.1, lines 135-136. For further clarity, definitions have also been added to the statements in Lemma 3.1 (Section 3.2). \n\n---\n\n> Q2. Lemma A.1 should be included in the main text. Some parts, such as detailed description of datasets can be moved to the appendix. \n\nI appreciate your valuable suggestions. We moved Lemma A.1 to the line 155 in the main and some details about the dataset to the appendix. \n\n---\n\n> Q3. The proposed study focuses solely on the expressivity of random FNO at initialization, a rather limited scenario. This limitation should be more explicitly emphasized, particularly in the title. \n\nThank you for your helpful suggestions. We are considering the following titles. Which proposed title adequately reflects the scenario limitation?\n* Understanding the expressivity and trainability of Fourier Neural Operator at initialization\n* Understanding the impact of initialization on expressivity and trainability in Fourier Neural Operator\n* Understanding the expressivity and trainability of Mean Field Fourier Neural Operator\n\n---\n\n> Q4. I am not convinced that the findings offer practical insights into the stable training of FNO, despite the preliminary numerical experiments provided. \n\nOur theory provides a universal insight, independent of the dataset, loss function, optimizer, and scheduler. While it does not offer sufficient conditions for stable training, it does provide a necessary condition for achieving stability in training. We believe that our results are of practical importance as they prevent unnecessary exploration of initialization. We have revisited the following. \n\t\n**line 44:** This discovery provides the necessary conditions (not sufficient conditions) for stable training, which is practically significant as it eliminates the unnecessary initialization searches.\n\t\n**line 310:** While we do not provide sufficient conditions for stable training, we do offer one necessary condition for achieving stable training. \n\n---\n\n> Q5. The term trainability, as used by the authors, only pertains to the vanishing or exploding gradients of random FNO.\n\nThe term trainability is used in the following related references with the same meaning as ours:\n* Xiao, Lechao, et al. \"Dynamical isometry and a mean field theory of CNNs: How to train 10,000-layer vanilla convolutional neural networks.\" \n* Yang, Ge, and Samuel Schoenholz. \"Mean field residual networks: On the edge of chaos.\" \n* Yang, Greg, et al. \"A Mean Field Theory of Batch Normalization.\" \n* Chen, Minmin, et al. \"Dynamical isometry and a mean field theory of RNNs: Gating enables signal propagation in recurrent neural networks.\"\n\nIn these references, **trainability** refers to a necessary condition for successful training that is determined by the architecture, independent of the dataset, loss function, optimizer, or scheduler, primarily focusing on discussions related to gradient vanishing/exploding. However, we acknowledge your point that, in general, the term **trainability** could be interpreted to encompass the overall optimization process and the success of training. Therefore, we have revised the following sentences in the introduction for clarity. Thank you for your objective and valuable feedback.\n\n**Line 29:** We analyze the exponential expressivity (how far representations of different positions can be pulled apart) and trainability (how much gradient vanishing/explosion on average) of the random FNO from the perspective of whether the network is ordered or chaotic.\n\n---\n\n> Q6. One might question the impact of the optimizer, the initial learning rate, and its scheduler on the results. An ablation study addressing these factors would have been beneficial.\n\nThank you for your valuable suggestions. We are currently conducting additional experiments. If the results are ready in time for the discussion period, we will report them in the thread. However, we would like to reiterate that our theory does not tackle the analysis of the optimization process and our results represent necessary conditions for stable training, independent of these factors. It is important to emphasize that the failure of learning due to inappropriate choices of optimizer, learning rate, or scheduler is not inconsistent with our findings.\n\n---\n\n> Q7. Does the fixed point of the iterated map always exist?\n\nLemma A.1 guarantees the existence of a fixed point for the FNO when the DCN has a fixed point with $(q^*, c^* = 1)$. The existence of the fixed point for the DCN depends on the activation function, so we need to check it on a case-by-case basis. For example, despite both being bounded activation functions, Tanh has a fixed point, whereas the sigmoid function does not. Among unbounded activation functions, ReLU has a fixed point, but softplus (a smoothed version of ReLU) does not. The fixed point analysis for bounded activation functions is conducted in [1], and the general case is discussed in [2].\n\n[1] Poole, Ben, et al. \"Exponential expressivity in deep neural networks through transient chaos.\"\n\n[2] Roberts, Daniel A., Sho Yaida, and Boris Hanin. \"The principles of deep learning theory. \"\n\n---\n\n> Q8. The representation of Fig.1(a) is unclear. ... \n\nFig. 1(a) experimentally verified Theorem 3.4. Instead of averaging over training epochs, It averages over multiple seeds at initialization (not over training epochs), with a relatively small standard deviation. We revised the first sentence of the title of Fig. 1 as follows.\n\n**fig1- title:** Average gradient norm $\\operatorname{Tr}(\\tilde{\\Sigma}^{(\\ell)})$ at initialization over multiple seeds during the backpropagation of several FNOs plotted as a function of layer $\\ell$", "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": "bUjcz0BomQ", "reviewer": "Reviewer_gSyo", "summary": "In this work, the authors investigate the expressive power and trainability of the Fourier Neural Operator (FNO). A mean-field theory is adapted for the FNO to examine the behavior of randomly initialized FNOs from the perspective of the 'edge of chaos'. The investigation focuses on understanding the expressive behavior of randomly initialized FNOs by examining the transition between ordered and chaotic phases. This phase transition demonstrates specific characteristics unique to the FNO, influenced by mode truncation, and shares similarities observed in densely connected networks and CNNs. Additionally, the authors establish a correlation between expressive power and trainability: the ordered and chaotic phases correspond to regions where gradients vanish or explode, respectively. The authors identify the necessity of initializing FNO near the edge of chaos for stable training theoretically and experimentally.", "questions": "**Major Questions:**\n\n1. If we treat FNO simply as a CNN, where global convolution occurs in the frequency domain, with Parseval theorems, it also implies, though less rigorously, that we can initialize the (complex) weights based on the He Distribution. This also hints on the similarities between FNOs and CNNs. Could the authors comment on this perspective?\n\n2. In line 122, it states:\n   > \"two learnable weights $\\Theta^{(\\ell, k)} \\in \\mathbb{R}^{D \\times D}$ and $\\Xi^{(\\ell, k)} \\in \\mathbb{R}^{D \\times D}$\". If I understand correctly, these are complex weights for spectral convolution.\n\n   Firstly, why is $\\Xi^{(\\ell, k)} \\in \\mathbb{R}^{D \\times D}$? Is this a typo where a $K$ is missing? Secondly, why are there two sets of learnable weights? Isn't there just one for the Fourier transform of the convolution kernel? Do you treat the real and complex parts separately, resulting in one set of parameters for the real and another for the complex parts?\n\n3. The weights in spectral convolution are complex. By initializing the weights according to $\\mathcal{N}\\left(0, \\frac{\\sigma^2}{4 D}\\right)$, do you treat the real and complex parts separately?\n\n4. The FNO is implemented with $L = 64$ layers, which is particularly deep for FNOs involving $64$ forward and backward FFTs. Is such a large number of layers practical in real-world applications? If not, how are your findings useful for shallower FNOs?\n\n**Minor/Optional Questions and Suggestions:**\n\n1. Some abbreviations are introduced before their full names, such as DCN in line 25, whereas its full name is introduced in line 37 (densely connected network, DCN).\n\n2. Given your insights into FNO expressivity and trainability, do you have any thoughts on potential architectural improvements or variations that could enhance FNO performance?\n\n3. It would be beneficial to provide a clearer explanation and discussion of the theoretical results and their implications for readers without a background in mean-field theory.", "limitations": "Yes, the authors adequately addressed the limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. **Theoretical novelty:** The paper adapts mean-field theory specifically for Fourier Neural Operators (FNOs), extending existing theories from other neural network architectures such as densely connected networks and CNNs.\n2. **Practical implications:** The theoretical findings lead to concrete recommendations for initializing FNOs to ensure stable training, particularly for deep FNOs.\n3. **Comprehensive experiments:** The authors validate their theoretical results across multiple PDE-solving tasks, including advection, Burgers, Darcy Flow, and the Navier-Stokes equations. All these datasets are publicly available.", "weaknesses": "1. The adaptation of mean-field theory to FNOs, while a notable endeavor, may not sufficiently break new ground in the field.\n2. The significance of the theoretical aspects that FNOs share similarities with CNNs should be further discussed.\nOverall, this work addresses how the network should be initialized for stable training of FNOs; however, its contribution might be somewhat limited. A stability analysis beyond the initialization phase, such as on the architecture, would be interesting to see.\n4. There are barely any experimental results shown in the main text, and the discussion of how their theoretical analysis is useful in practical applications is lacking. The discovery in their work, the He distribution for FNOs, has already been applied in prior works such as [1] (see their GitHub implementation, `reset_parameters()` method under models); hwoever, it is always good to have a rigorous justification, which is shown in this work.\n\n[1] Helwig, J., Zhang, X., Fu, C., Kurtin, J., Wojtowytsch, S., and Ji, S. Group equivariant Fourier neural operators for partial differential equations.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Due to the character limit, we have addressed some questions in the global rebuttal. Please refer to our global rebuttal. We appreciate your attention and review.\n\n> Q1. The discovery in their work, ...\n\nThe following part of the code is different from our He initialization, as it only initializes using the ‘nn.init.kaiming_uniform_’ function. The weights are initialized by a uniform distribution with different scales and the bias is not fixed to zero. \n\n```python\ndef reset_parameters(self):\n        if self.Hermitian:\n            for v in self.W.values():\n                nn.init.kaiming_uniform_(v, a=math.sqrt(5))\n        else:\n            nn.init.kaiming_uniform_(self.W, a=math.sqrt(5))\n        if self.B is not None:\n            nn.init.kaiming_uniform_(self.B, a=math.sqrt(5))\n```\n---\n\n> Q2. If we treat FNO simply as a CNN ...\n\nThank you very much for your insightful and valuable comments. We added this discussion of the similarities between CNN and FNO to the appendix. The following is our answer to this question and also the part of appendix.\n\nThe iterated map for CNNs (Xiao et al., 2018) is obtained as follows:\n$$\n\\Sigma^{(\\ell)}\\_{\\alpha, \\alpha'} = \\frac{\\sigma^2}{2r+1} \\sum\\_{\\beta \\in \\mathrm{ker}} \\mathbb{E}\\_{\\boldsymbol{H}\\_{:. d} \\sim \\mathcal{N}(\\boldsymbol{0}, \\boldsymbol{\\Sigma}^{(\\ell-1)})}\\left[ \\phi(H\\_{\\alpha+\\beta, d})\\phi(H\\_{\\alpha'+\\beta, d}) \\right] + \\sigma\\_b^2 \n$$\nwhere $2r + 1$ is the number of filter width and $\\operatorname{ker} = \\{ \\beta \\in \\mathbb{Z} \\mid |\\beta| \\leq r\\}$ is the set of indices referring to the elements of the filter.\nWhen we consider the FNO without mode truncation ($K=N/2 + 1$), the propagation of the diagonal components $\\Sigma\\_{\\alpha, \\alpha}$ for any $\\alpha \\in [N]$ is equivalent to the propagation in a CNN with filter size $N$ performing global convolution.\n\n$$\n\\begin{aligned}\n \\Sigma^{(\\ell)}\\_{\\alpha, \\alpha} &= \\sigma^2 \\sum\\_{k=0}^{\\frac{N}{2}} c\\_k  \\mathbb{E}\\_{\\boldsymbol{H}\\_{:. d} \\sim \\mathcal{N}(\\boldsymbol{0}, \\boldsymbol{\\Sigma}^{(\\ell-1)})}\\left[\\left|\\left[\\mathcal{F}\\phi\\left(\\boldsymbol{H}\\_{:, d}\\right)\\right]\\_{k}\\right|^2 \\right] \\underbrace{\\cos\\left( \\theta^{(k)}\\_{\\alpha, \\alpha}\\right)}\\_{=1}  + \\sigma^2\\_b \\\\\\\\\n    &= \\frac{\\sigma^2}{N} \\sum\\_{\\beta=0}^{N-1} \\mathbb{E}\\_{\\boldsymbol{H}\\_{:. d} \\sim \\mathcal{N}(\\boldsymbol{0}, \\boldsymbol{\\Sigma}^{(\\ell-1)})}\\left[ \\left|\\phi(H\\_{\\beta, d})\\right|^2\\right] + \\sigma^2\\_b \\\\\\\\\n    &= \\frac{\\sigma^2}{N} \\sum\\_{\\beta=0}^{N-1} \\mathbb{E}\\_{\\boldsymbol{H}\\_{:. d} \\sim \\mathcal{N}(\\boldsymbol{0}, \\boldsymbol{\\Sigma}^{(\\ell-1)})}\\left[ \\phi(H\\_{\\alpha+\\beta, d})\\phi(H\\_{\\alpha+\\beta, d})\\right] + \\sigma^2\\_b.\n\\end{aligned}\n$$\nThe first equality follows from Perseval's equality and the second from the periodic boundary condition. In contrast, the propagation of the off-diagonal components differs between CNN and FNO, and the iterated map is different even in the presence of mode truncation. Similarly, the iterated map during the backpropagation for FNO can be recognized as the global CNN (see Eq. (2.16) in (Xiao et al., 2018)) with $v\\_{\\beta} = \\frac{1}{N}$.\n$$\n\\tilde{\\Sigma}\\_{\\alpha, \\alpha}^{(\\ell)} \\approx \\frac{1}{N^2} \\sum\\_{\\beta, \\beta'=0}^{N-1} \\tilde{\\Sigma}\\_{\\beta, \\beta'}^{(\\ell+1)}  \\chi\\_{c^*} \\underbrace{\\left(\\sum\\_{k=0}^{N/2} c\\_k \\cos\\left( \\theta^{(k)}\\_{\\beta-\\beta', \\alpha-\\alpha} \\right) \\right)}_{=N \\delta\\_{\\beta, \\beta'}} = \\chi\\_{c^*} \\sum\\_{\\beta=0}^{N-1} \\frac{1}{N} \\tilde{\\Sigma}\\_{\\beta, \\beta}^{(\\ell+1)}.\n$$\nThis equivalence suggests that the edge of chaos initialization (e.g., He initialization) is also valid for the FNO since the problems of gradient vanishing and explosion are determined by the diagonal components of $\\tilde{\\boldsymbol{\\Sigma}}$.\n\n(Xiao et al., 2018) Xiao, Lechao, et al.\"Dynamical isometry and a mean field theory of cnns: How to train 10,000-layer vanilla convolutional neural networks.\"\n\n---\n\n> Q3.  In line 122, it states: ... & The weights in spectral convolution ...\n\nFor ease of handling in proofs, we define the complex weights $R^{(\\ell)} \\in \\mathbb{C}^{D \\times D \\times K}$ of the $\\ell$-th layer of the FNO by splitting them into their real and complex parts. We further divide these weights along the mode dimension using the index k, as $R^{(\\ell)} = \\bigg|\\bigg|\\_{k=1}^K \\Theta^{(\\ell, k)} + \\sqrt{-1} \\Xi^{(\\ell, k)}$ where $\\bigg|\\bigg|$ denotes the concatenation operation along the mode axis k.\nFor all $\\ell \\in [L]$ and $k \\in [K]$, every element of the weights $\\Theta^{(\\ell, k)}$ and $\\Xi^{(\\ell, k)}$ is independently initialized from \\$\\mathcal{N}(0, \\sigma^2/2D)$.\n\n---\n\n> Q4: The FNO is implemented with 𝐿=64 layers,  ... Is such a large number of layers practical in real-world applications?\n\nNetworks with a large number of layers are practical for real-world applications. For instance, as demonstrated in Appendix D.2, deeper FNOs with 16 or 32 layers achieved better performance than the conventional 4-layer FNO in learning complex operators. The study by Tran et al. (2023) on F-FNO also explores the training of relatively deep FNOs. In the field of PDE solving, it is possible to generate a large amount of data through numerical computation. Therefore, we believe that in the future, architectures with larger model capacities will be employed for more complex operator approximations and the construction of large pre-trained models as in the fields of CV or NLP. Our theory provides valuable insights for such multi-layer FNOs and architectures that partially use FNO layers.\n\n---\n\n> Q5. Some abbreviations ...\n\nThank you for pointing this out. We have revised this point.\n\n---\n\n> Q6. Given your insights into FNO expressivity and trainability, do you have any thoughts on potential architectural improvements or variations that could enhance FNO performance?\n\nIt seems difficult to obtain suggestions for new architectures. We may be able to give theoretical guarantees for newly proposed architectures.", "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:
> Q1. **(reviewer ePA1)** I recommend providing further motivation for all concepts and referencing previous works where appropriate.  
**(reviewer gSyo)** It would be beneficial to provide a clearer explanation and discussion of the theoretical results and their implications for readers without a background in mean-field theory. 

We have included additional texts on the motivation for fixed point analysis and references to proofs located in their respective appendices. Is this explanation sufficiently clear for our readers? If you have any further suggestions or areas where you believe improvements could be made, please do not hesitate to inform us. Your feedback is greatly appreciated.

**line 152:** By linearizing the dynamics of signal propagation around this fixed point and analyzing the stability and rate of convergence to the fixed point, we can determine the depth to which each component of the input can propagate. 

**line 155:** To analyze the stability and convergence rate, we linearly approximate the C-map around the fixed point $\Sigma^*$, i.e., $\mathcal{C}(\Sigma) \approx \Sigma^* + J\_{\Sigma^*} (\Sigma - \Sigma^*)$, where $J\_{\Sigma^*}$ is the Jacobian linear map of the iterated map defined in Eq. (15).

**line 161:** From Lemma A.4, $K-1$ matrices in $\{ \psi^{(k)} \}\_{k \in [K] \backslash \{0\}}$ are eigenbases with the eigenvalue $\chi_{c^*}$ of the Jacobian linear map. From Lemma A.5, the matrix $\psi$ is the eigenbases with the eigenvalue $\chi$ of the Jacobian linear map. 
Since the rank of the Jacobian linear map is at most K (Lemma A.3), the deviation from the fixed point $\Sigma^{(\ell)} - \Sigma^*$ is spanned by K-dimensional eigenspace $\operatorname{span}\left(\\{\psi^{(k)} \\}_{k \in [K] \backslash \{0\}} \cup \\{ \psi \\} \right)$. 

We have also added a discussion on the similarity with CNN to improve the interpretability of the results of Theorem 3.3 and 3.4. 

---

> Q2. **(reviewer gSyo)** The significance of the theoretical aspects that FNOs share similarities with CNNs should be further discussed. 

We have revised the text to emphasize that FNO has similarities not only with DCN but also with CNN. In particular, we have revised the text in Section 3.2 as follows. 

**line 154:** (Xiao et al., 2018) showed that any fixed point for the iterated map of the DCN is also a fixed point for that of CNN. 

**line 155:** Lemma A.1 means that the fixed point for the iterated map $\Sigma^*$ of the DCN serves as a fixed point for the iterated map of the simplified FNO (as well as CNN). 

**line 178:** The convergence rates $\chi_{q^*}$ and $\chi_{c^*}$ are the same as the convergence rates of the variance and correlation to the fixed point for DCN (Schoenholz et al., 2016) and CNN (Xiao et al., 2018).

**line 182:** CNN possess diagonal eigenspaces associated with eigenvalues $\chi_{q^*}$ and non-diagonal eigenspaces associated with eigenvalues $\chi_{c^*}$. In contrast, FNOs without mode truncation exhibit a similarity, possessing eigenspaces $\chi_{q^*}$ for zero-frequency and eigenspaces $\chi_{c^*}$ for k-frequencies with diagonal components removed.

**line 200:** Despite the architectural and iterative map differences between FNO, DCN, and CNN, Theorems 3.3 and 3.4 demonstrate the similarity in the random behavior of FNO, DCN, and CNN.

Meta-review: The article investigates the expressivity and trainability of Fourier Neural Operator. 

Reviewers find the article addresses a significant question and that it offers theoretical depth, novel insight, practical implications, comprehensive experiments, and clarity. In particular, that it offered rigorous analysis of FNO and insights into practical aspects for training. Criticism included that the discussion of related work and discussion of certain experimental setups could be improved, that parts of the work were not self contained, and that the code was not made available. Authors offered additions to the references, additional experiments, and also explained that the article does not propose any new models or learning methods and the code used is based on other existing code. The responses clarified many of the questions and motivated reviewers to maintain positive scores or increase them. 

At the end of the discussion period all reviewers recommend (borderline, weak, strong) accept. In view of the positive assessment at conclusion of the discussion, I recommend accept.
Decision: Accept (poster)