1 Introduction↩︎

Fourier neural operators (FNOs) have emerged as a powerful architecture for operator learning, where the goal is to learn maps between infinite-dimensional function spaces from data. Many important operators in scientific computing arise as solution maps of evolution equations: an input function \(u_0\), representing the initial state of the system, is mapped to the corresponding state \(u(T)\) at a later time \(T\). In this paper, we focus on time-\(T\) solution maps \(u_0\mapsto u(T)\) associated with dissipative evolution equations. Despite the empirical success of FNOs, a key theoretical question remains: under what structural conditions can FNOs approximate and learn such solution operators efficiently?

This paper addresses this question through a guiding principle rooted in numerical analysis: we show that if a solution operator admits a stable and accurate Fourier spectral discretization, then it can also be represented efficiently by an FNO and learned from data with polynomial sample complexity. Thus, rather than treating solution operators as arbitrary maps between function spaces, we exploit the structure inherited from numerical schemes for the underlying evolution equation. This perspective identifies stable spectral approximability as a broad structural condition that enables FNOs to learn nonlinear evolution operators efficiently.

To make this principle precise, we study solution operators generated by dissipative evolution equations on the \(d\)-dimensional torus \(\mathbb{T}^d\) of the form \[\frac{du}{dt}+\mathcal{A}u+\mathcal{D}_s\bigl(\mathcal{G}(u)\bigr)=f,\qquad t>0,\] with initial condition \(u(0) = u_0\). Here \(\mathcal{A}\) is a dissipative linear operator diagonalizable in the Fourier basis, \(\mathcal{D}_s\) is a differential operator, \(\mathcal{G}\) is a nonlinear reaction or transport term, and \(f\) is a forcing term. For fixed \(T>0\), the associated time-\(T\) solution operator is the map \(\mathcal{S}: u_0 \mapsto u(T)\). Our goal is to prove quantitative approximation and learning guarantees for such operators using FNOs.

The analysis proceeds in four steps. First, we introduce classes of semi-implicit Fourier spectral schemes that discretize these evolution equations in space and time. Second, we define classes of solution operators that can be approximated by iterating stable and accurate spectral steps. Third, we construct FNOs that emulate these spectral steps and compose them to approximate the full time-\(T\) solution map. This yields quantitative FNO approximation bounds. Finally, we combine these approximation estimates with metric-entropy bounds to establish learning guarantees for empirical risk minimization over suitable FNO classes.

A key feature of our results is that they hold uniformly over operator classes, rather than for a single fixed partial differential equation (PDE). For equations with polynomial nonlinearities, the learning rates depend primarily on the Sobolev smoothness of the input space and the dimension of the physical domain. We verify that the framework applies to the time-\(T\) solution maps of the Navier–Stokes, Allen–Cahn, and Cahn–Hilliard equations. We also extend the theory to nonlinearities satisfying general smoothness assumptions. In this more general setting, polynomial sample complexity still holds, and the rates additionally depend on the smoothness of the nonlinearity and the strength of dissipation.

1.1 Related Works↩︎

For overviews of the operator learning literature, we refer to [1]–[4] and the references therein. Here, we review relevant literature on FNOs and statistical guarantees for operator learning.

1.1.0.1 Fourier Neural Operators

The FNO architecture was introduced in [5] for learning PDE solution operators over the \(d\)-dimensional torus. Further variants have been proposed to handle more complicated geometries [6]–[8]. FNOs have been used effectively in applications such as numerical weather prediction [9], biomedical design [10], plasma modeling [11], and carbon capture and storage [12]. On the theory side, FNOs were shown to have the universal approximation property and to efficiently approximate the Darcy flow and Navier–Stokes solution operators in [13]. Approximation error bounds for general operator learning architectures, including FNOs, were shown for nonlinear parabolic PDEs in [14]. The work [15] derives quantitative approximation guarantees with neural operators for a class of nonlinear evolution equations defined via Picard iteration. The discretization error of FNOs was characterized in [16], and bounds on the Rademacher complexity of FNOs were established in [17]. In the special case of linear FNOs, [18] decomposes the error into statistical, discretization, and truncation components and bounds each component for a particular least-squares estimator.

1.1.0.2 Statistical Guarantees for Operator Learning

The sample complexity of operator learning tasks has received significant attention in a variety of settings. For holomorphic operators, it has been shown that the number of samples required to approximate to \(\varepsilon\) accuracy scales algebraically in \(\varepsilon^{-1}\) [19], [20]. In contrast, learning \(k\)-times Fréchet differentiable or Lipschitz operators requires exponentially many samples in the worst case. This negative result holds for input measures with support on a compact set [21] and with Gaussian input measures [22]. Consequently, understanding the success of operator learning requires developing assumptions beyond those of classical smoothness. Operator Barron spaces [23] may be one such assumption; however, sample complexity guarantees for such operators have yet to be developed.

Several positive results are known for more structured operator classes. Optimal rates for learning linear operators between Sobolev spaces are established in [24]. Restricting to the case of diagonalizable linear operators, [25] characterizes the sample complexity and provides generalization bounds under distribution shifts. Random feature methods have been shown to achieve Monte Carlo estimation rates if the target operator lies in a reproducing kernel Hilbert space [26]. In [27], the authors again consider learning Lipschitz operators and demonstrate that classes of encoder-decoder networks achieve the optimal, but slow, logarithmic convergence rate. They further show that, if the input function space is restricted to be finite-dimensional, then the sample complexity scales polynomially in the dimension of the input space.

Other work focuses on guarantees specifically for learning PDEs. For elliptic PDEs, [28] proves that the required sample complexity is logarithmic, meaning that the error decreases exponentially fast with the sample size. It is important to note, however, that [28] does not consider noisy observational data. Instead, the authors assume access to PDE queries with deliberately chosen randomized inputs and corresponding exact solution outputs. For linear parabolic PDEs and in a similar setting, [29] demonstrates an algorithm achieving a polynomial sample complexity that is increasing in the dimension of the physical domain. In the more challenging setting of linear hyperbolic PDEs, [30] proves that polynomially many samples can still suffice. Convergence rates for learning pseudo-differential operators from noisy input-output data are derived in [31]. There has been less work on nonlinear PDEs. In [32], the authors demonstrate that the solution operator for viscous Burgers’ equation, along with several linear PDEs, can be learned with polynomial sample complexity. The work most similar to our own is [33], which considers learning time-stepping schemes for evolution equations with an encoder-decoder neural network. In contrast to our work, [33] assumes that the training data is generated from the numerical methods, not from the true PDE.

1.1.0.3 Spectral Methods

We refer to [34] for a classical reference on spectral methods. More recently, and of particular relevance to our results, there has been analysis of so-called “large time-stepping” implicit-explicit (IMEX) methods for dissipative PDEs that do not require restrictive CFL conditions. For the 2D Navier–Stokes equations, the IMEX Euler scheme was shown to be convergent with a step-size requirement independent of the spatial discretization for initial data in \(H^2\) [35], and with mild dependence on the spatial discretization for initial data in \(L^2\) or \(H^1\) [36]. Similar convergence guarantees were derived for a stabilized IMEX Euler scheme for the Cahn–Hilliard equation in [37]. However, these bounds are conditional on the stabilization term being large enough relative to the \(L^{\infty}\) norm of the numerical solution, which is not known a priori.

In [38], the authors demonstrate that an IMEX Euler scheme for the Cahn–Hilliard equation with logarithmic potentials is unconditionally energy stable. They develop a bootstrap technique to obtain a priori control over the \(L^{\infty}\) norm of the numerical solution, which in turn yields unconditional large time-stepping guarantees. Their analysis leverages recent tools from harmonic analysis [39]. These techniques have subsequently been used to obtain unconditional large time-stepping guarantees for the standard Cahn–Hilliard equation [40] and the Allen–Cahn equation [41].

1.2 Main Contributions and Outline↩︎

We now summarize the main contributions of the paper and indicate where they are established. All required background on FNOs will be provided in Section 2.

• A spectral-method-based operator class. In Section 3, we introduce a class of time-\(T\) solution operators defined through stable and accurate semi-implicit Fourier spectral approximations. This class is designed to capture the numerical structure shared by many dissipative evolution equations, while remaining broad enough to include several standard nonlinear PDEs.

• Uniform FNO approximation bounds for polynomial nonlinearities. For equations whose nonlinearities are polynomial, we prove that every operator in the resulting class can be uniformly approximated by an FNO. The construction emulates one step of the spectral method and then composes the resulting FNOs to approximate the time-\(T\) solution map. The resulting architecture sizes are quantified in terms of the desired accuracy and the parameters defining the operator class. This is proved in Theorem 2.

• Polynomial sample complexity for FNO empirical risk minimization. We combine the constructive approximation bounds with metric-entropy estimates to prove learning guarantees for empirical risk minimization over FNO classes. In the polynomial-nonlinearity setting, the resulting rates depend on the input space smoothness and the physical dimension. This yields polynomial sample complexity for learning all operators in the class; see Theorem 3.

• Applications to classical dissipative PDEs. In Section 3.4, we verify that the time-\(T\) solution maps of the Navier–Stokes, Allen–Cahn, and Cahn–Hilliard equations fit into the abstract framework and apply Theorem 3. This gives concrete learning guarantees for these PDEs and illustrates how the assumptions arise from standard spectral stability and convergence estimates.

• Extension to general smooth nonlinearities. In Section 4, we extend the framework beyond polynomial nonlinearities. Under non-parametric smoothness assumptions on \(\mathcal{G}\), we again obtain FNO approximation bounds and polynomial sample complexity. The resulting rates are generally slower and depend additionally on the smoothness of the nonlinearity and the strength of the dissipation; see Theorems 5 and 6. In Section 4.4, we apply this theory to learn the time-\(T\) solution map of the Cahn-Hilliard equation with a logarithmic potential. This equation presents additional difficulties for approximation and learning due to the possibility of singularities in the potential.

2 Fourier Neural Operators↩︎

This section introduces the class of FNOs considered in this paper. Intuitively, an FNO \(\Psi\) with \(L\) hidden layers is a nonlinear map \(\Psi: L_N^2 (\mathbb{T}^d, \mathbb{R}^{d_{\rm in}}) \to L_N^2(\mathbb{T}^d, \mathbb{R}^{d_{\rm out}})\) of the form \[\begin{align} \Psi = \mathcal{Q} \circ \mathcal{L}_L\circ \cdots \circ \mathcal{L}_1 \circ \mathcal{R}, \end{align}\] where the operator \(\mathcal{R}\) lifts the space \(L_N^2 (\mathbb{T}^d, \mathbb{R}^{d_{\rm in}})\) of \(\mathbb{R}^{d_{\rm in}}\)-valued degree-\(N\) trigonometric polynomials on the torus \(\mathbb{T}^d = [0,2\pi]^d\) to a latent space \(L_N^2 (\mathbb{T}^d, \mathbb{R}^{d_{\rm lat}}),\) the hidden layers \(\mathcal{L}_\ell\), \(1 \le \ell \le L\), map the latent space to itself, and \(\mathcal{Q}\) maps the latent space to \(L_N^2(\mathbb{T}^d, \mathbb{R}^{d_{\rm out}}).\) In Section 2.1, we introduce the function spaces and operators involved in the lifting and projection layers \(\mathcal{R}\) and \(\mathcal{Q}\), as well as the Fourier hidden layers \(\mathcal{L}_\ell\). We then formally introduce the class of FNOs in Section 2.2.

2.1 Background: Function Spaces and Operators↩︎

FNOs perform finite-dimensional function approximation through finite Fourier representations. The basic idea is to represent functions by finitely many Fourier modes, perform operations on these finite representations, and then move between Fourier coefficients and grid values using the discrete Fourier transform. Hence, we begin by defining finite-dimensional function spaces consisting of finite Fourier series. Throughout, we use \(N \in \mathbb{N}\) to denote the Fourier cutoff for these spaces. For \(N\in \mathbb{N}\), let \(\mathsf{K}_N :=\{k\in \mathbb{Z}^d:\|k\|_{\infty}\leq N\}\) and let \[\begin{align} \label{eq:trigpoly} L^2_N(\mathbb{T}^d) : = \biggl\{v_N:\mathbb{T}^d\to \mathbb{R}\, \Big\vert \, v_N(x)=\sum_{k \in \mathsf{K}_N}c_ke^{i\langle k,x\rangle} \, \forall x \in \mathbb{T}^d,\, c_{-k}=c_k^* \, \, \forall k \in \mathsf{K}_N \biggr\} \end{align}\tag{1}\] be the space of real-valued trigonometric polynomials of degree \(N\) equipped with the usual \(L^2\) inner product. The condition \(c_{-k}=c_k^*\) in 1 guarantees that \(v_N(x) \in \mathbb{R}\) for all \(x\in \mathbb{T}^d.\) The space \(L_N^2(\mathbb{T}^d, \mathbb{R}^{{d_*}})\) of \(\mathbb{R}^{{d_*}}\)-valued trigonometric polynomials is defined as the space of functions from \(\mathbb{T}^d\) to \(\mathbb{R}^{d_*}\) with components in \(L^2_N(\mathbb{T}^d).\)

Consider the set of indices \(\mathsf{J}_N=\{0,1,\ldots,2N\}^d\) and the corresponding uniform grid on the torus \(\{x_j\}_{j\in \mathsf{J}_N}\), where \(x_j=\frac{2\pi j}{2N+1}\). Note that \(|\mathsf{J}_N|=|\mathsf{K}_N|=(2N+1)^d.\) The discrete Fourier transform (DFT) is the map \(\mathcal{F}_N:L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\to \mathbb{C}^{(2N+1)^d\times {d_*}},\) which assigns to \(v_N \in L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\) the matrix \(\bigl(\mathcal{F}_N v_N(k) \bigr)_{k \in \mathsf{K}_N} \in \mathbb{C}^{(2N+1)^d\times {d_*}}\) with rows \[\begin{align} \label{eq:DFT} (\mathcal{F}_N v_N)(k):=\frac{1}{(2N+1)^d}\sum_{j\in \mathsf{J}_N}v_N(x_j)e^{-i \langle x_j , k\rangle} =: \hat{v}_N(k) \in \mathbb{C}^{d_*}, \quad \forall\: k \in \mathsf{K}_N. \end{align}\tag{2}\] Denoting \(\vec{v}_N :=\{v_N(x_j)\}_{j\in \mathsf{J}_N},\) the DFT can be written as a map (matrix) \({\boldsymbol{F}} :\mathbb{R}^{(2N+1)^d\times {d_*}}\to \mathbb{C}^{(2N+1)^d\times {d_*}}\) from function values at the grid points to Fourier coefficients, given row-wise by \[\begin{align} ({\boldsymbol{F}} \vec{v}_N)(k) = \frac{1}{(2N+1)^d}\sum_{j\in \mathsf{J}_N}v_N(x_j)e^{-i \langle x_j , k\rangle} =:\hat{v}_N(k) \in \mathbb{C}^{d_*}, \quad \forall k \in \mathsf{K}_N. \end{align}\] The inverse DFT \(\mathcal{F}_N^{-1}:\mathbb{C}^{(2N+1)^d\times {d_*}}\to L^2_N(\mathbb{T}^d,\mathbb{C}^{{d_*}})\) is given by \[\begin{align} (\mathcal{F}_N^{-1} \hat{v}_N)(x)=\sum_{k\in \mathsf{K}_N}\hat{v}_N(k)e^{i \langle k,x\rangle}, \quad \forall x\in \mathbb{T}^d. \end{align}\] As a map from conjugate-symmetric Fourier coefficients to function values on the uniform grid \(\{x_j\}_{j\in \mathsf{J}_N}\), we write the inverse DFT as \({\boldsymbol{F}} ^{-1}:\mathbb{C}^{(2N+1)^d\times {d_*}}\to \mathbb{R}^{(2N+1)^d\times {d_*}}.\) We have that, row-wise, \[\begin{align} ({\boldsymbol{F}} ^{-1}\hat{v}_N)(j) :=\sum_{k\in \mathsf{K}_N}\hat{v}_N(k)e^{i \langle k,x_j \rangle} \in \mathbb{R}^{d_*}, \quad \forall j\in \mathsf{J}_N. \end{align}\]

We next introduce an interpolation operator to formalize the usual pseudo-spectral identification between a function and its values on the uniform grid. Specifically, we let \(\mathcal{I}_N:C(\mathbb{T}^d)\to L_N^2(\mathbb{T}^d)\) be the pseudo-spectral interpolation operator mapping \(v \in C(\mathbb{T}^d)\) to the unique trigonometric polynomial \(\mathcal{I}_N v\in L^2_N(\mathbb{T}^d)\) that satisfies \[\begin{align} \label{eq:interpolationoperator} (\mathcal{I}_N v)(x_j)=v(x_j), \quad \forall j\in \mathsf{J}_N. \end{align}\tag{3}\] Extending the DFT to \(C(\mathbb{T}^d)\) by the same formula as in 2 , we can write, for \(v \in C(\mathbb{T}^d),\) \[\mathcal{I}_Nv=\mathcal{F}_N^{-1}\mathcal{F}_Nv.\] For vector-valued functions, the operator \(\mathcal{I}_N\) is applied component-wise. Note that for \(v_N\in L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\), \(\mathcal{I}_Nv_N=v_N.\) Hence, we can identify functions \(v_N\in L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\) with their values at the grid points \(\{x_j\}_{j \in \mathsf{J}_N}.\) See [34] for a detailed discussion of Fourier pseudo-spectral approximation.

2.2 Fourier Neural Operators↩︎

Let \(d_{\rm in},d_{\rm out},d_{\rm lat}\in \mathbb{N}\) be given. A Fourier neural operator (FNO) with \(L\in \mathbb{N}\) hidden layers is a nonlinear map \(\Psi:L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm in}})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm out}})\) of the form \[\begin{align} \Psi = \mathcal{Q} \circ \mathcal{L}_L\circ \cdots \circ \mathcal{L}_1 \circ \mathcal{R}, \end{align}\] where:

1. \(\mathcal{R}: L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm in}})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) and \(\mathcal{Q}:L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\to L^2_N(\mathbb{T}^{d}, \mathbb{R}^{d_{\rm out}})\) are linear “lifting” and “projection” layers mapping \(u \in L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm in}})\) and \(v \in L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) as follows: \[\begin{align} (\mathcal{R}u)(x)= {\boldsymbol{R}} u(x), \quad \forall x\in \mathbb{T}^d, \\ (\mathcal{Q}v)(x)={\boldsymbol{Q}} v(x), \quad \forall x\in \mathbb{T}^d, \end{align}\] where \({\boldsymbol{R}}\in \mathbb{R}^{d_{\rm lat}\times d_{\rm in}}\) and \({\boldsymbol{Q}}\in \mathbb{R}^{d_{\rm out}\times d_{\rm lat}}\) are learnable matrices.

2. \(\mathcal{L}_{\ell}:L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\), for \(1 \le \ell \le L,\) are nonlinear hidden layers mapping \(v \in L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) as follows: \[\mathcal{L}_{\ell}(v)(x)=\mathcal{I}_{N}\sigma\bigl( \mathcal{W}_{\ell}v+ \mathcal{K}_{\ell}(v)+b_{\ell}\bigr)(x), \quad \forall x \in \mathbb{T}^d.\] Here:

   1. \(\mathcal{W}_{\ell}:L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) is linear and maps \(v \in L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) as follows: \[\begin{align} (\mathcal{W}_{\ell} v)(x)={\boldsymbol{W}}_{\ell}v(x), \quad \forall x\in \mathbb{T}^d, \end{align}\] where \({\boldsymbol{W}}_{\ell}\in \mathbb{R}^{d_{\rm lat}\times d_{\rm lat}}\) is a learnable matrix.

   2. \(\mathcal{K}_{\ell}:L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) is a nonlocal linear operator mapping \(v \in L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) as follows: \[(\mathcal{K}_{\ell}v)(x)=\Bigl(\mathcal{F}^{-1}_N\bigl(\hat{P}_{\ell}(\mathcal{F}_N v)\bigr)\Bigr)(x)=\sum_{k\in \mathsf{K}_N}\hat{P}_{\ell}(k)\hat{v}(k)e^{i \langle k,x\rangle},\] where \(\mathcal{F}_N\) is the DFT and \(\mathcal{F}_N^{-1}\) its inverse, and the \(\hat{P}_\ell(k)\in \mathbb{C}^{d_{\rm lat}\times d_{\rm lat}}\) for \(k\in \mathsf{K}_N\) are a learnable sequence of matrices required to satisfy \(\hat{P}_{\ell}(-k)=\mkern 1.5mu\overline{\mkern-1.5mu\hat{P}_{\ell}(k)\mkern-1.5mu}\mkern 1.5mu\) for all \(k\in \mathsf{K}_N.\)

   3. \(b_{\ell}\in L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) is a learnable function.

   4. \(\sigma:L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\to C(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}})\) is a ReLU activation function, which is applied component-wise and point-wise; that is, for \(v \in L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm lat}}),\) \[\begin{align} \sigma(v)(x)=\begin{bmatrix} \max\{0,v_1(x)\}\\ \vdots \\ \max\{0,v_{d_{\rm lat}}(x)\} \end{bmatrix} \,, \qquad \forall x\in \mathbb{T}^d. \end{align}\]

   5. \(\mathcal{I}_N:C(\mathbb{T}^d,\mathbb{R}^d_{\rm lat})\to L_N^2(\mathbb{T}^d,\mathbb{R}^d_{\rm lat})\) is the pseudo-spectral projection defined in Section 2.1.

We have defined FNOs as acting on functions in \(L^2_N(\mathbb{T}^d,\mathbb{R}^{d_{\rm in}}).\) For a general continuous function, an FNO acts on its trigonometric interpolant: \(\Psi(u)=\Psi(\mathcal{I}_N u)\) for \(u\in C(\mathbb{T}^d,\mathbb{R}^{d_{\rm in}})\). We characterize the size of an FNO \(\Psi\) through the following quantities:

1. The depth \(\mathscr{D}(\Psi):=L;\)

2. The channel width \(\mathscr{W}(\Psi):=\max \bigl\{d_{\rm in},d_{\rm out},d_{\rm lat}\bigr\};\)

3. The weight magnitude \(\mathscr{S}(\Psi):=\max_{k\in \mathsf{K}_N,1\leq \ell \leq L}\{\|{\boldsymbol{R}}\|_{\infty},\|{\boldsymbol{Q}}\|_{\infty}, \|{\boldsymbol{W}}_{\ell}\|_{\infty},\|\hat{P}_k^{\ell}\|_{\infty},\,\|\hat{b}^{\ell}_k\|_{\infty}\}\);

4. The Fourier cutoff, \(\mathscr{C}(\Psi) = N\).

We define \(\mathsf{FNO}(\mathscr{D},\mathscr{W},\mathscr{S},\mathscr{C})\) to be the set of FNOs with depth \(\mathscr{D}\), width \(\mathscr{W}\), and Fourier cutoff \(\mathscr{C}\) with weights constrained to be bounded by \(\mathscr{S}\). We will refer to the vector of all trainable parameters in an FNO as \(\theta=\text{vec}({\boldsymbol{R}}, {\boldsymbol{W}}_1,\hat{P}_1^1,\hat{b}_1^1,\ldots,\hat{P}_{ (2N+1)^d }^L,\hat{b}_{(2N+1)^d}^L,{\boldsymbol{Q}})\in \mathbb{R}^{d_{\theta}}\), where \(d_{\theta}\leq 5(2N+1)^dLd_{\rm lat}^2\).

3 Dissipative Models with Polynomial Nonlinearity↩︎

This section investigates how well FNOs can approximate and learn the time-\(T\) solution operator for evolution equations on the \(d\)-dimensional torus \(\mathbb{T}^d\) of the form \[\begin{align} \label{eq:evolution32PDE32again} \begin{aligned} \frac{du}{dt}+\mathcal{A}u+\mathcal{D}_s\bigl(\mathcal{G}(u)\bigr) &=f, \qquad t>0. \end{aligned} \end{align}\tag{4}\] We assume \(\mathcal{A}\) to be a dissipative linear operator diagonalizable in the Fourier basis, \[\begin{align} \mathcal{A}e^{i\langle k,x\rangle}=\lambda_k(\mathcal{A})e^{i\langle k, x\rangle}, \quad \forall k\in \mathbb{Z}^{d}, \end{align}\] with eigenvalues scaling as \[\begin{align} \label{eq:dissipation32scaling} \lambda_k(\mathcal{A})\asymp \|k\|_{2}^{\beta}, \quad \forall k\in \mathbb{Z}^{d}, \end{align}\tag{5}\] for some \(\beta>0.\) For instance, \(\mathcal{A}\) may be a negative Laplacian \((\beta= 2)\), a bilaplacian \((\beta= 4)\), or a fractional Laplacian. We assume that the nonlinearity \(\mathcal{G}\) is a degree-\(p\) polynomial acting on \(u\), and that \(\mathcal{D}_s\) is a differential operator with each coordinate derivative order bounded by \(s.\) General smooth nonlinearities, not necessarily of polynomial type, will be considered in Section 4.

A key observation is that equations of the form 4 can be efficiently solved by spectral methods. This motivates us to introduce in Section 3.1 a class of operators that admit accurate spectral approximations. For this operator class, which includes solution maps of many equations of the form 4 , we establish FNO approximation bounds in Section 3.2 and learning bounds in Section 3.3. Applying these general results, we obtain in Section 3.4 learning bounds for the solution operators of three important evolution PDEs: the Navier–Stokes, Allen–Cahn, and Cahn–Hilliard equations.

3.1 Spectral Methods and Class of Operators↩︎

In this section, we introduce the general class of operators for which we will obtain FNO approximation and learning bounds. This class, defined in Subsection 3.1.3, consists of operators that can be well approximated by the family of spectral methods introduced in Subsection 3.1.2. We begin with background on spectral methods in Subsection 3.1.1.

3.1.1 Background on Spectral Methods↩︎

To solve an equation of the form 4 numerically, one must discretize the problem in both space and time. Semi-implicit Fourier spectral methods are widely used and broadly applicable to dissipative equations of the form 4 . A first-order-in-time semi-implicit method is given by \[\begin{align} \label{eq:generic32spectral32method32unstabilized} \frac{1}{\tau}\bigl(u_N^{j+1}-u_N^j\bigr)+\mathcal{A}u_N^{j+1}+ \mathcal{P}_N\mathcal{D}_s\bigl(\mathcal{G}(u_N^j)\bigr)-f_N=0, \end{align}\tag{6}\] for a spatial resolution \(N\in \mathbb{N}\) and timestep \(\tau>0.\) Here \(\mathcal{P}_N\) denotes an \(L^2\)-orthogonal projection onto \(L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\) or onto a shift-invariant subspace of \(L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\). The subspaces used in our examples in Section 3.4 include spaces of zero-mean functions and divergence-free functions. Our learning-theoretic guarantees in Section 3.3 do not require a priori knowledge of the subspace. Finally, we define \(f_N=\mathcal{P}_Nf\) to be the \(L^2\)-orthogonal projection of the forcing function \(f\) onto the same subspace.

Notice that in 6 , the linear dissipation term is treated implicitly, while the nonlinearity is treated explicitly. Semi-implicit methods typically have better stability than explicit methods, while avoiding the fixed-point iterations required by fully implicit methods [42]. For both theoretical and practical reasons, semi-implicit methods are often implemented with an additional stabilization term controlled by parameters \(\gamma\geq 0\) and \(\eta\in [0,1),\) leading to the update \[\begin{align} \label{eq:generic32spectral32method} \frac{1}{\tau}\bigl(u_N^{j+1}-u_N^j\bigr)+\gamma\mathcal{A}^{\eta}\bigl(u_N^{j+1}-u_N^j\bigr)+\mathcal{A}u_N^{j+1}+ \mathcal{P}_N\mathcal{D}_s\bigl(\mathcal{G}(u_N^j)\bigr)-f_N=0. \end{align}\tag{7}\] Equation 7 defines a mapping \(u_N^j\mapsto u_N^{j+1}=:\mathcal{S}_{\tau,N}(u_N^j).\) Thus, for any \(u_N \in L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}}),\) we define \[\begin{align} \label{eq:spec} \begin{aligned} \mathcal{S}_{\tau,N}(u_N) :&=\Bigl(\frac{1}{\tau}I+\gamma\mathcal{A}^{\eta}+\mathcal{A}\Bigr)^{-1}\Bigl((\frac{1}{\tau}I+\gamma\mathcal{A}^{\eta})u_N-\mathcal{P}_N\mathcal{D}_s\bigl(\mathcal{G}(u_N)\bigr)+f_N\Bigr)\\ & = \mathcal{B}\Bigl(\bigl(\frac{1}{\tau}I+\gamma \mathcal{A}^{\eta}\bigr)u_N- \mathcal{D}_s\bigl(\mathcal{G}(u_N)\bigr)+f_N\Bigr), \end{aligned} \end{align}\tag{8}\] where \(\mathcal{B}:= \Bigl(\frac{1}{\tau}I +\gamma\mathcal{A}^{\eta}+ \mathcal{A}\Bigr)^{-1} \mathcal{P}_N.\) We refer to \(\mathcal{S}_{\tau,N}\) as a single step of the spectral method. To approximate the time-\(T\) solution operator of the PDE, one can choose \(J \in \mathbb{N}\) and \(\tau>0\) such that \(J=T/\tau,\) and then iterate this single-step map: \[\begin{align} u(T) = \mathcal{S}(u) \approx \mathcal{S}_{\tau,N}^J(\mathcal{I}_Nu(0))=\underbrace{\mathcal{S}_{\tau,N}\circ \cdots \circ \mathcal{S}_{\tau,N}}_{J \text{ times}}(u_N^0). \end{align}\] In the next subsection, we introduce the class of single-step spectral methods considered in this work. Our operator class, introduced in Subsection 3.1.3, consists of operators that can be well approximated by iterating these single-step spectral methods.

3.1.2 Family of Single-Step Spectral Methods↩︎

Recall the generic single-step spectral method \(\mathcal{S}_{\tau,N}\) introduced in 8 . In this subsection, we define a family of single-step spectral methods by imposing conditions on the linear operator \(\mathcal{B},\) the lower-order differential operator \(\mathcal{D}_s,\) the nonlinearity \(\mathcal{G},\) and the forcing \(f_N\).

For timestep \(\tau>0\) and degree \(N \in \mathbb{N},\) we introduce parameter \(\vartheta:=( \beta,\gamma,\eta,c, s,C_{\mathcal{D}},p, C_{\mathcal{G}}, F)\) and define the family of single-step spectral methods \[\begin{align} \mathsf{S}_{\tau,N}(\vartheta) :=\Biggl\{\mathcal{S}_{\tau,N}:L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}}) \, \bigg| & \, \mathcal{S}_{\tau,N}(u_N)=\mathcal{B}\Bigl(\bigl(\frac{1}{\tau}I+\gamma\mathcal{A}^{\eta} \bigr)u_N-\mathcal{D}_s\bigl(\mathcal{G}(u_N)\bigr)+f_N\Bigr), \\ &\mathcal{B}\in \mathsf{B}(\beta,\gamma,\eta,c), \, \mathcal{D}_s\in \mathsf{D}(s,C_{\mathcal{D}}), \, \mathcal{G}\in \mathsf{G}(p,C_{\mathcal{G}}), \, \, f_N\in \mathsf{F}(F) \Biggr\}. \end{align}\]

In the above, \(\mathsf{B}(\beta,\gamma,\eta,c)\) is a class of linear operators, \(\mathsf{D}(s,C_{\mathcal{D}})\) is a class of differential operators, \(\mathsf{G}(p,C_{\mathcal{G}})\) is a class of polynomial nonlinearities, and \(\mathsf{F}(F)\) is a class of forcings. Before formally defining these classes, we provide an intuitive interpretation of their parameters. The parameters \((\beta,\gamma,\eta,c)\) in \(\mathsf{B}\) control respectively the rate of dissipation, the strength of the stabilization term, the order of the stabilization term, and an upper bound on the inverse term; the parameters \((s,C_{\mathcal{D}})\) in \(\mathsf{D}\) represent the order of derivatives and a bound on their coefficients; the parameters \((p,C_{\mathcal{G}})\) in \(\mathsf{G}\) represent the degree of the polynomial nonlinearity and a bound on its coefficients; and the parameter \(F\) in \(\mathsf{F}\) represents a pointwise bound on the forcing term. The following table summarizes the values of these parameters for our PDE examples in Section 3.4.

We now define the classes \(\mathsf{B}, \mathsf{G},\) \(\mathsf{D},\) and \(\mathsf{F}\):

1. \(\mathsf{B}(\beta,\gamma,\eta,c)\) is a class of linear operators defined for \(\beta>0,\) \(\gamma \ge 0,\) \(\eta\in [0,1)\) and \(c>0.\) It contains operators \(\mathcal{B}:L^2(\mathbb{T}^d,\mathbb{R}^{{d_*}})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\) of the form \[\mathcal{B}= \Bigl(\frac{1}{\tau}I+\gamma \mathcal{A}^{\eta} + \mathcal{A}\Bigr)^{-1} \mathcal{P}_N,\] where:

   1. \(\mathcal{P}_N : L^2(\mathbb{T}^d,\mathbb{R}^{d_*}) \to L^2_{N}(\mathbb{T}^d,\mathbb{R}^{d_*})\) is an \(L^2\)-orthogonal projection onto any shift invariant subspace of \(L^2_N(\mathbb{T}^{d},\mathbb{R}^{{d_*}})\), so that for any \(u\in L^2(\mathbb{T}^d,\mathbb{R}^{{d_*}}),\) \[\begin{align} (\mathcal{P}_Nu)(x)=\sum_{k\in \mathsf{K}_N}\hat{P}_{\mathcal{P}_N}(k)\hat{u}(k)e^{i\langle k,x\rangle}, \end{align}\] where \(\hat{P}_{\mathcal{P}_N}(k)\in \mathbb{C}^{{d_*}\times {d_*}}\) are orthogonal projections \(\hat{P}_{\mathcal{P}_N}(k)\in \mathbb{C}^{{d_*}\times {d_*}}\) satisfying \(\hat{P}_{\mathcal{P}_N}(k)=\mkern 1.5mu\overline{\mkern-1.5mu\hat{P}_{\mathcal{P}_N}(-k)\mkern-1.5mu}\mkern 1.5mu\) for all \(k\in \mathsf{K}_N.\)

   2. \(\mathcal{A}:\text{Range}(\mathcal{P}_N)\subseteq L^2_{N}(\mathbb{T}^d,\mathbb{R}^{{d_*}})\to L^2_{N}(\mathbb{T}^d,\mathbb{R}^{{d_*}})\) is positive semi-definite, diagonalizable in Fourier space, and its eigenvalues satisfy 5 with parameter \(\beta,\) so that \[\begin{align} (\mathcal{A}u)(x)=\sum_{k\in \mathsf{K}_N}\hat{P}_{\mathcal{A}}(k)\hat{u}(k)e^{i\langle k,x\rangle}, \end{align}\] for positive definite matrices \(\hat{P}_{\mathcal{A}}(k)\in \mathbb{C}^{{d_*}\times {d_*}}\) satisfying \(\hat{P}_{\mathcal{A}}(k)=\mkern 1.5mu\overline{\mkern-1.5mu\hat{P}_{\mathcal{A}}(-k)\mkern-1.5mu}\mkern 1.5mu\), \(\hat{P}_{\mathcal{A}}(k)\hat{P}_{\mathcal{P}_N}(k)=\hat{P}_{\mathcal{P}_N}(k)\hat{P}_{\mathcal{A}}(k)\), and \[\begin{align} \label{eq:A32scaling32assumption} 0\prec \|k\|_2^{\beta}c\preceq \hat{P}_{\mathcal{A}}(k)\preceq \|k\|_2^{\beta}I, \end{align}\tag{9}\] for some \(c\in (0,1]\) for all \(k\neq 0\) in \(\mathsf{K}_N\) and \(\hat{P}_{\mathcal{A}}(0)=0.\)

2. \(\mathsf{D}(s,C_{\mathcal{D}}),\) with \(\beta/2 \ge s \in \mathbb{Z}_+,\) denotes the class of differential operators with each coordinate derivative order bounded by \(s\) given by \[\begin{align} \mathsf{D}(s,C_{\mathcal{D}})& := \Biggl\{\mathcal{D}_s:L^2_{pN}(\mathbb{T}^d,\mathbb{R}^{\tilde{d}})\to L^2_{pN}(\mathbb{T}^d,\mathbb{R}^{{d_*}})\bigg|\\ & \mathcal{D}_s(v)(x)=\begin{bmatrix} \sum_{\ell=1}^{\tilde{d}} \sum_{j=1}^{d}c_{j\ell 1}\partial_{x_j}^{s_{j\ell 1}}v_{\ell}(x)\\ \vdots \\ \sum_{\ell=1}^{\tilde{d}}\sum_{j=1}^dc_{j\ell {d_*}}\partial_{x_j}^{s_{j\ell{d_*}}} v_{\ell}(x) \end{bmatrix}, \, \|\vec{c}\|_{\infty}\leq C_{\mathcal{D}}, \,\|\vec{s}\|_{\infty}\leq s \,\Biggr\} \, , \end{align}\] where \(\vec{c}=(c_{111},\ldots, c_{d\tilde{d}{d_*}})\in \mathbb{R}^{d\tilde{d}{d_*}}\) and \(\vec{s}=(s_{111},\ldots,s_{d\tilde{d}{d_*}})\in \mathbb{Z}_{+}^{d\tilde{d}{d_*}}.\) Here \(\|\vec{c}\|_{\infty}\) denotes the largest coefficient and \(\|\vec{s}\|_{\infty}\) denotes the highest order derivative in \(\mathcal{D}_s.\)

3. \(\mathsf{G}(p,C_{\mathcal{G}})\) is the space of degree-\(p\) polynomial nonlinearities (applied component-wise and point-wise) of \({d_*}\) variables with coefficients bounded by \(C_{\mathcal{G}}\); that is, \[\begin{align} \mathsf{G}(p,C_{\mathcal{G}})&:=\Biggl\{\mathcal{G}:L^2_{N}(\mathbb{T}^d,\mathbb{R}^{{d_*}})\to L^2_{pN}(\mathbb{T}^d,\mathbb{R}^{\tilde{d}}) \, \bigg\vert \, \\ & \mathcal{G}(u_N)(x)=\begin{bmatrix} g_{1}(u_N(x))\\ \vdots \\ g_{\tilde{d}}(u_N(x)) \end{bmatrix} \,\,\, \forall x\in \mathbb{T}^d, g_{j}\in \mathsf{POLY}(p,C_{\mathcal{G}}), 1 \le j \le \tilde{d}\Biggr\}, \end{align}\] where \[\begin{align} \mathsf{POLY}(p,C_{\mathcal{G}}):=\biggl\{g:\mathbb{R}^{{d_*}}\to \mathbb{R}\, \Big|\, g(z)= c_{0}+\sum_{i_1=1}^{{d_*}}c_{i_1}z_{i_1}+\cdots+ \sum_{i_1,\ldots, i_p=1}^{{d_*}}c_{i_1,\ldots, i_p}z_{i_1}\ldots z_{i_p},\|\vec{c}\|_{\infty} \leq C_{\mathcal{G}}\biggr\} \end{align}\] and \(\|\vec{c}\|_{\infty}\) denotes the magnitude of the largest coefficient.

4. \(\mathsf{F}(F)\) denotes the class of forcing terms in the range of \(\mathcal{P}_N\) that are bounded pointwise by \(F:\) \[\begin{align} \mathsf{F}(F)&:=\biggl\{f_N\in \text{range}(\mathcal{P}_N) \Big\vert \sup_{x\in \mathbb{T}^d}\|f_N(x)\|_2\leq F\biggr\}. \end{align}\]

The conditions \(0\preceq \|k\|_2^{\beta}cI \preceq \hat{P}_{\mathcal{A}}(k)\preceq \|k\|_2^{\beta}\), and \(\hat{P}_{\mathcal{A}}(k)\hat{P}_{\mathcal{P}_N}(k)=\hat{P}_{\mathcal{P}_N}(k)\hat{P}_{\mathcal{A}}(k)\) for all \(k\in \mathsf{K}_N\) imply that for any mean-zero \(u\in \text{Range}(\mathcal{P}_N)\subseteq L^2_{N}(\mathbb{T}^d,\mathbb{R}^{{d_*}})\), \(\mathcal{A}\) satisfies the Poincaré inequality \[\begin{align} c\|u\|_{L^2}\leq \|\mathcal{A}u\|_{L^2}. \end{align}\] It is straightforward to extend the class \(\mathsf{D}(s,C_{\mathcal{D}})\) to include mixed partial derivatives, and our main results in Theorems 2 and 3 would continue to hold with only minor modifications to the proofs. For ease of notation, we do not pursue this extension, since the examples in Section 3.4 do not involve mixed derivatives.

3.1.3 Class of Operators↩︎

Next, we introduce the class of operators that can be well approximated by spectral methods and then show that this class can also be well approximated by FNOs. We consider operators acting on a Sobolev ball of radius \(U,\) denoted by \[\mathsf{U}_a:=\bigl\{u\in L^2(\mathbb{T}^d,\mathbb{R}^{{d_*}}):\|u\|_{H^a}\leq U\bigr\},\] where the Sobolev norm is defined via Fourier series as \[\begin{align} \|u\|_{H^{a}}^2:=\sum_{k\in \mathbb{Z}^d}\bigl(1+\|k\|_2^2\bigr)^a\|\hat{u}(k)\|_{2}^2. \end{align}\] For fixed \(T>0\) and parameter \(\bar{\vartheta}:=(\vartheta,a,U,V),\) with \(a>d/2,\) we define \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) to be the class of operators \(\mathcal{S}\) for which the following condition holds:

The parameter \(a>d/2\) will play a key role in our approximation and learning bounds. One can interpret this parameter as the Sobolev smoothness of the solution to a (semi)-dynamical system of the form 4 . For dissipative evolution equations, the dynamics are typically asymptotically bounded, so Condition 1 [item:boundedness] is natural. We assume, for simplicity of notation, that the solution upper bound in Condition 1 [item:boundedness] is the same as the spectral method stability bound appearing in the first item of Condition 1 [item:spectral32method]; if this is not the case, one can take the maximum of the two upper bounds. The second item in Condition 1 [item:spectral32method] captures the rate at which the spectral method converges to the true solution operator, where the spatial rate is determined by the smoothness parameter \(a\).

We emphasize that semi-implicit spectral methods \(\mathcal{S}_{\tau,N}\) typically do not require a restrictive CFL condition, as is often the case for explicit time-stepping schemes. For this reason, the accuracy bound in Condition 1 [item:spectral32method] does not require \(\tau\) to scale with \(N\). We will show in Section 3.4 that the time-\(T\) solution maps of the Navier–Stokes, Allen–Cahn, and Cahn–Hilliard equations belong to the class \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) for appropriate choices of the parameter \(\bar{\vartheta}.\)

3.2 Approximation Bound↩︎

In this section, we state our first main result: operators in \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) can be uniformly approximated to accuracy \(\varepsilon\) by sufficiently large FNOs. In particular, we provide quantitative estimates on the required network size. We show that any \(\mathcal{S}\in \mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) can be approximated to accuracy \(\varepsilon\) within the following class of FNOs:

\[\begin{align} \label{eq:sigmaclass} \begin{aligned} \Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})=\biggl\{\Psi\in \mathsf{FNO}(\mathscr{D},\mathscr{W},\mathscr{S},\mathscr{C}): \mathscr{D}&=r\varepsilon^{-2}\log (\varepsilon^{-1}),\\ \mathscr{W}&=r\tilde{d}{p+{d_*}-1\choose {d_*}-1},\\ \mathscr{S}&=r\varepsilon^{-{\max\{1,\frac{\beta \eta}{a},\frac{s}{a}\}}},\\ \mathscr{C}&=r\varepsilon^{-\frac{1}{a}},\\ &\sup_{u\in \mathsf{U}_a}\|\Psi(u)\|_{L^2}\leq 2V\biggr\}. \end{aligned} \end{align}\tag{10}\] Here, \(r\) is a positive constant, increasing in \(\beta^{-1},\gamma,\eta,c^{-1}, s,p, C_{\mathcal{G}}, F,U\) and \(V\), which can be traced through the proof of Theorem 2 in Section 8.1. Note that we also constrain the output of the FNO to lie in an \(L^2\)-ball of radius \(2V.\) We are now ready to state the first main result.

Theorem 2 provides a uniform approximation guarantee over the class of operators \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta}),\) with the required size of the FNO architecture controlled quantitatively in terms of the desired accuracy \(\varepsilon\), as specified in the definition of \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\). While quantitative FNO approximation bounds have been shown for specific PDEs, such as stationary Darcy flow and the Navier–Stokes equations [13], the constructions have been specific to the equations at hand. To our knowledge, Theorem 2 is the first approximation result to hold uniformly over a broad class of operators, including those arising from a range of practically relevant nonlinear evolution equations. We discuss particular examples of PDEs covered by our theory in Section 3.4. This approximation result is essential to the learning bound in Section 3.3, where we establish polynomial sample complexity for learning operators in \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) with FNOs.

We postpone the proof of Theorem 2 to Section 8.1. The idea is to construct an FNO that efficiently approximates one step of the spectral method for an appropriately chosen timestep \(\tau\) and spatial resolution \(N\), and then to compose \(j=\lceil T/\tau\rceil\) such single-step approximations. The resulting composition can be implemented as a single FNO with \(j\) times as many layers. The construction is inspired by, but more general than, that in [13], which is specific to the Navier–Stokes equations. The error analysis differs more significantly. In particular, we measure the single-step error in the nonlinearity in an \(\mathcal{A}^{-1/2}\)-weighted norm, which is both well suited to analyzing the error of \(j\) compositions and allows for faster approximation rates by exploiting the dissipation induced by \(\mathcal{A}\).

The proof proceeds as follows. In Lemmas 14, 16, and 17, we construct FNOs to approximate \(\mathcal{B}\in \mathsf{B}(\beta,\gamma,\eta,c),\) \(\mathcal{D}_s \in \mathsf{D}(s,C_{\mathcal{D}})\), and \(\mathcal{G}\in \mathsf{G}(p,C_{\mathcal{G}})\), respectively. Combining these FNOs allows us to efficiently approximate any spectral method in \(\mathsf{S}_{\tau,N}(\vartheta),\) and in turn approximate operators in \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\). Our proof technique for controlling the error incurred by iterating the single-step FNO differs from that of [13] for the Navier–Stokes case, and more closely resembles stability arguments for spectral methods, such as those in [35].

3.3 Learning Bound↩︎

We assume that we are given training data \(\bigl\{u^{(m)},\mathcal{S}(u^{(m)})\bigr\}_{m=1}^M\) for some unknown \(\mathcal{S}\in \mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\), where \(u^{(1)},\ldots,u^{(M)}\overset{\mathrm{i.i.d.}}{\sim} \mu,\) and \(\mu\) is a measure with \(\operatorname{supp}(\mu)\subseteq \mathsf{U}_a\). We consider learning with an FNO trained by empirical risk minimization. In particular, we consider the squared-error empirical risk functional \[\begin{align} \label{eq:empirical32risk} \hat{\mathcal{L}}(\Psi,\mathcal{S}):=\frac{1}{M}\sum_{m=1}^{M}\|\Psi( u^{(m)} )-\mathcal{S}( u^{(m)} )\|_{L^2}^2, \end{align}\tag{11}\] and the corresponding population risk \[\begin{align} \label{eq:population32risk} \mathcal{L}(\Psi,\mathcal{S}):=\mathbb{E}_{u\sim \mu}\|\Psi(u)-\mathcal{S}(u)\|_{L^2}^2. \end{align}\tag{12}\] Our learned estimator is given by a minimizer of 11 over the class of FNOs \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\), defined in the previous section, for some choice of \(\varepsilon\): \[\begin{align} \label{eq:ERMfortheorem} \hat{\Psi}_{\mathcal{S}}\in \mathop{\mathrm{arg\,min}}_{\Psi\in \Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})}\hat{\mathcal{L}}(\Psi,\mathcal{S}). \end{align}\tag{13}\] We refer to any such \(\hat{\Psi}_{\mathcal{S}}\) as an empirical risk minimizer (ERM). In the following remark, we briefly explain why such empirical risk minimizers exist. We denote by \(C(\mathsf{U}_{a},L^2(\mathbb{T}^{d},\mathbb{R}^{{d_*}}))\) the set of continuous operators from \(\mathsf{U}_{a}\) to \(L^2(\mathbb{T}^{d},\mathbb{R}^{{d_*}})\) equipped with the supremum norm. By definition of \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\subset C(\mathsf{U}_{a},L^2(\mathbb{T}^{d},\mathbb{R}^{{d_*}}))\), the trainable parameters of the FNO satisfy \(\theta\in[-r\varepsilon^{-1},r\varepsilon^{-1}]^{d_{\theta}}\), which is compact. Moreover, by [21], the parameter-to-operator map \(\theta\mapsto \Psi_{\theta}\) is Lipschitz continuous. Thus, \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\) is compact in \(C(\mathsf{U}_{a},L^2(\mathbb{T}^{d},\mathbb{R}^{{d_*}})),\) and a minimizer of 11 is guaranteed to exist.

The main result of this section, Theorem 3, shows that FNOs obtained by minimizing the empirical risk 11 over \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\), for an appropriate choice of \(\varepsilon\), learn operators in \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) with polynomial sample complexity.

This theorem demonstrates that operators in \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) can be learned with polynomial sample complexity using FNOs. In contrast, [21] shows that general \(k\)-times Fréchet differentiable or Lipschitz operators cannot be learned from a polynomial number of samples, regardless of the learning procedure. In our setting, operators in \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) can be well approximated by spectral methods and consequently admit parameter-efficient representations that can be learned data-efficiently by FNOs. As we will discuss in detail in Section 3.4, the class \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) contains the solution maps of a range of nonlinear dissipative evolution equations.

We present the proof of Theorem 3 in Section 8.2. The argument follows the general structure of [21], using ideas similar to those in [21], [43]. At a high level, the idea is classical. We first show that the population risk of an ERM can be controlled by the best approximation error and the worst-case empirical risk deviation over \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\) and \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta}).\) Since we minimize over \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta}),\) the approximation error is controlled by \(\varepsilon\) by construction. The remaining task is therefore to show that, with high probability, the worst-case empirical risk deviations are uniformly small. This is achieved by controlling the metric entropies of \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\) and \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta}).\) Finally, \(\varepsilon\) is chosen to balance the approximation error against the metric entropy of \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\).

In the setting considered in this section, where the nonlinearity is given by a degree-\(p\) polynomial, the statistical rate achieved by the FNO ERM is determined by the smoothness parameter \(a\) and by the dimension \(d\) of the physical domain \(\mathbb{T}^d\). This is because the depth required for an FNO to approximate a degree-\(p\) polynomial to accuracy \(\varepsilon\) scales like \(r\log(\varepsilon^{-1})\), where only the constant factor \(r\) increases with \(p\); see Lemma 6. The contribution of this part of the approximation to the metric entropy of \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\) therefore only introduces polynomial factors in \(r\log(\varepsilon^{-1})\). Iterating this construction \(j=\lceil T/\tau\rceil\) times and taking the input domain to be sufficiently large contributes \(r\varepsilon^{-4}\text{poly}\log(\varepsilon^{-1})\) to the metric entropy, where the parameters in \(\bar{\vartheta}\) only enter through the constant \(r\). The width contributes only a constant factor, and the weight bound only a logarithmic factor in \(\varepsilon^{-1}\) and a constant factor in \(\bar{\vartheta}\) (see Lemma 22). In contrast, the error incurred by truncating an \(H^{a}\) function to its first \(N\) Fourier modes scales like \(N^{-a}\), so achieving accuracy \(\varepsilon\) requires \(N\asymp \varepsilon^{-1/a}\). The corresponding contribution to the metric entropy scales like \(\varepsilon^{-\frac{d}{a}}\), which is the primary parametric dependence in the bound. In Section 4, we consider nonlinearities \(\mathcal{G}\) that satisfy general smoothness assumptions. The resulting rates for these operators depend explicitly on the smoothness of the nonlinearity and the strength of the dissipation, and are generally slower than the rates for polynomial nonlinearities. See Remark [remark:comparison32remark] for further discussion.

3.4 Examples: Learning the Solution Map of Dissipative PDEs↩︎

In this section, we verify that several classical dissipative evolution equations fit into the abstract framework developed above. We focus on the Navier–Stokes, Allen–Cahn, and Cahn–Hilliard equations, which illustrate how the assumptions defining \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) arise naturally for physically relevant PDEs. Together, these examples exhibit distinct nonlinear structures, dissipation mechanisms, and stabilization regimes. They are particularly well suited to FNO approximation because their dynamics admit efficient Fourier spectral discretizations on periodic domains.

3.4.1 Navier–Stokes Equations↩︎

The incompressible Navier–Stokes equations describe the evolution of viscous incompressible fluid flows under the combined effects of advection and diffusion. On the two-dimensional torus \(\mathbb{T}^2\) with periodic boundary conditions, the equations take the form \[\begin{align} \label{eq:NSE} \frac{\partial u}{\partial t}- \nu \mathcal{P}^{ \rm LH} \Delta u+\mathcal{P}^{ \rm LH} \bigl((u\cdot \nabla)u \bigr)=f \, . \end{align}\tag{14}\] Here \(u(x,t)\in\mathbb{R}^2\) denotes the velocity field, \(\nu\in(0,1]\) is the viscosity, \(\Delta\) is the Laplacian, and \(f\in L^\infty(\mathbb{T}^2,\mathbb{R}^2)\) is a divergence-free, mean-zero forcing term. The operator \(\mathcal{P}^{ \rm LH}\) denotes the Leray–Helmholtz projection onto divergence-free vector fields, enforcing the incompressibility constraint. The Laplacian induces second-order dissipation, corresponding to \(\beta=2,\) while the quadratic transport term gives the degree-two polynomial nonlinearity appearing in our general framework.

Let \(\mathcal{S}^{\mathrm{NS}}: L^2(\mathbb{T}^2, \mathbb{R}^2)\to L^2(\mathbb{T}^2,\mathbb{R}^2)\) denote the time-\(T\) solution map of the Navier–Stokes equations, which exists and is unique [44]. We now verify Condition 1 to show that \(\mathcal{S}^{\mathrm{NS}}\in \mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) with parameters \[\vartheta=(\beta,\gamma,\eta,c, s,C_{\mathcal{D}},p, C_{\mathcal{G}}, F)=(2,0,0,\nu,1,1,2,1,F)\] and \[\bar{\vartheta}=(\vartheta,a,U,V)=(\vartheta,2,U,U+\frac{F}{\nu})\] for any desired \(U,F>0\) and \(\nu\in (0,1]\).

We begin with the boundedness of \(\mathcal{S}^{\mathrm{NS}}\). Given a pointwise upper bound \(F>0\) on the forcing function \(f\), [45] shows that \[\begin{align} \sup_{u:\|u\|_{H^2}\leq U}\|\mathcal{S}^{\mathrm{NS}}(u)\|_{L^2}\leq U\exp(-\nu 2\pi^2T)+\frac{F}{\nu}\leq U+\frac{F}{\nu}:=V. \end{align}\] We consider the semi-implicit Euler scheme \[\begin{align} \label{eq:NS32euler} \frac{1}{\tau}\bigl(u_N^{j+1}-u_N^{j}\bigr)- \nu \Delta u_N^{j+1}+\mathcal{P}_N^{\mathrm{NS}} \nabla \cdot \bigl( u_N^j \otimes u_N^j\bigr) =\mathcal{P}_N^{\mathrm{NS}} f, \end{align}\tag{15}\] where \(\mathcal{P}_N^{\mathrm{NS}}:L^2(\mathbb{T}^2,\mathbb{R}^2)\to L^2_N(\mathbb{T}^2,\mathbb{R}^2)\) denotes the \(L^2\)-orthogonal projection onto the space of degree-\(N\), divergence-free, mean-zero trigonometric polynomials. This scheme can be written in the form of 8 with

1. \(\mathcal{B}^{\mathrm{NS}}:L^2(\mathbb{T}^2,\mathbb{R}^2)\to L^2_N(\mathbb{T}^2,\mathbb{R}^2)\) given by \[\begin{align} \mathcal{B}^{\text{NS}}&=\Bigl(\frac{1}{\tau} I+ \mathcal{A}^{\text{NS}}\Bigr)^{-1}\mathcal{P}_N, \end{align}\] which corresponds to \(\gamma=\eta=0\), where

   1. \(\mathcal{P}_N^{\mathrm{NS}}:L^2(\mathbb{T}^2,\mathbb{R}^2)\to L^2_N(\mathbb{T}^2,\mathbb{R}^2)\) is given by \[\begin{align} (\mathcal{P}_N^{\mathrm{NS}}u)(x)=\sum_{ k\in \mathsf{K}_N\backslash\{0\}}\Bigl(I-\frac{k\otimes k}{\|k\|_2^2}\Bigr)\hat{u}(k)e^{i\langle k,x\rangle}, \end{align}\] for any \(u\in L^2(\mathbb{T}^2,\mathbb{R}^2).\)

   2. \(\mathcal{A}^{\mathrm{NS}}:\text{Range}(\mathcal{P}_N^{\mathrm{NS}})\subset L^2_N(\mathbb{T}^2,\mathbb{R}^2)\to L^2_N(\mathbb{T}^2,\mathbb{R}^2)\) is given by \[\begin{align} (\mathcal{A}^{\mathrm{NS}}u_N)(x)=(-\nu\Delta u_N)(x)=\nu\sum_{ k\in \mathsf{K}_N\backslash\{0\}}\|k\|_2^2\hat{u}_N(k)e^{i\langle k,x\rangle}, \end{align}\] for any \(u_N\in \text{Range}(\mathcal{P}_N^{\mathrm{NS}})\), which implies \(\beta=2\) and \(c=\nu.\)

   Consequently, \(\mathcal{B}^{\mathrm{NS}}\in \mathsf{B}(2,0,0,\nu)\).

2. \(\mathcal{D}_s^{\mathrm{NS}}:L^2_{2N}(\mathbb{T}^2,\mathbb{R}^3)\to L^2_{2N}(\mathbb{T}^2,\mathbb{R}^2)\) given by \[\begin{align} \mathcal{D}_s^{\mathrm{NS}}(v_{2N})=\begin{bmatrix} \partial_{x_1}v_{2N}^1+\partial_{x_2}v_{2N}^2\\ \partial_{x_1}v_{2N}^2+\partial_{x_2}v_{2N}^3 \end{bmatrix}, \end{align}\] for any \(v_{2N}\in L^2_{2N}(\mathbb{T}^2,\mathbb{R}^3),\) which implies \(s=1\) and \(C_{\mathcal{D}}=1\), so \(\mathcal{D}_s^{\mathrm{NS}}\in \mathsf{D}(1,1).\)

3. \(\mathcal{G}^{\mathrm{NS}}:L^2_{N}(\mathbb{T}^2,\mathbb{R}^2)\to L^2_{2N}(\mathbb{T}^2,\mathbb{R}^3)\) given by \[\begin{align} \mathcal{G}^{\mathrm{NS}}(u_N)=\begin{bmatrix} (u_N^1)^2\\ u_N^1 u_N^2\\ (u_N^2)^2 \end{bmatrix}, \end{align}\] for any \(u_N\in L^2_{N}(\mathbb{T}^2,\mathbb{R}^2),\) which implies \(p=2\) and \(C_{\mathcal{G}}=1\). Consequently, \(\mathcal{G}^{\mathrm{NS}}\in \mathsf{G}(2,1).\)

It remains to verify the stability and accuracy requirements in Condition 1. In [35], this scheme is shown to satisfy a uniform \(L^2\) stability bound, verifying the first part of Condition 1 (II). In particular, they show that for any \(u_0\in \mathsf{U}_2\), if \(\tau\) is chosen sufficiently small relative to \(\nu\), \(F\), and \(U,\) then, for all \(j\geq 1,\)¹ \[\begin{align} \|\bigl(\mathcal{S}_{\tau,N}^{\text{NS}}\bigr)^j(\mathcal{I}_N u_0)\|_{L^2}\leq U+\frac{F}{\nu}=V. \end{align}\] Further, [35] shows the method yields the error bound \[\begin{align} \|u(j\tau)-\bigl(\mathcal{S}_{\tau,N}^{\text{NS}}\bigr)^j(\mathcal{I}_N u_0)\|_{L^2}\leq C_1e^{C_2 j\tau}\Bigl(N^{-2}+\tau\Bigr), \end{align}\] for \(\|u_0\|_{H^2}\leq U\), if \(\tau\) is sufficiently small relative to a constant depending only on \(U\), \(F\) and \(\nu\), where \(C_1\) and \(C_2\) also only depend on \(\nu\), \(F\), and \(U\). This verifies the second part of Condition 1 [item:spectral32method] with \(a=2\). Consequently, the Navier–Stokes solution operator belongs to \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\), and therefore satisfies the approximation and learning guarantees established in Theorem 3. We thus obtain the following corollary:

The rate \(M^{-1/8}\) results from the spatial smoothness parameter and the dimension of the physical space, which are \(a=2\) and \(d=2\), respectively. The constant \(r\), which we do not specify here but can be traced through the proof of Theorem 3, is increasing in \(\nu^{-1},F\) and \(U,\) capturing the natural intuition that the Navier–Stokes dynamics become harder to learn as the viscosity decreases or the forcing strength increases.

3.4.2 Allen–Cahn Equation↩︎

The Allen–Cahn equation is a reaction-diffusion model describing phase separation dynamics driven by interfacial energy minimization. On the two-dimensional torus \(\mathbb{T}^2\) and with periodic boundary conditions, it is given by \[\begin{align} \label{eq:Allen32Cahn} \frac{\partial u}{\partial t}-\nu \Delta u+u^3-u=0, \end{align}\tag{16}\] where \(u(x,t)\) represents a mixture of two phases [46]. For initial data satisfying \(-1\leq u_0\leq 1\), the maximum principle ensures that the solution remains in this range. The Laplacian induces second-order dissipation corresponding to \(\beta=2,\) while the cubic reaction term drives phase separation and produces the degree-three polynomial nonlinearity appearing in our general framework.

We let \(\mathcal{S}^{\mathrm{AC}}:L^2(\mathbb{T}^2, \mathbb{R})\to L^2(\mathbb{T}^2,\mathbb{R})\) denote the time-\(T\) solution map of the Allen–Cahn equation, mapping \(u_0\mapsto u(T)=\mathcal{S}^{\mathrm{AC}}(u_0).\) We now verify that \(\mathcal{S}^{\mathrm{AC}}\in \mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) for parameters \[\vartheta=(\beta,\gamma,\eta,c, s,C_{\mathcal{D}},p, C_{\mathcal{G}}, F)=(2,r(U^2+\nu^{-1}|\log(\nu)|^2),0,\nu,0,1,3,1,0)\] and \[\bar{\vartheta}=(\vartheta,a,U,V)=(\vartheta,4,U,\sqrt{ 2\nu U^2+r^4U^4+9\pi^2})\] for any \(U>0\) and \(\nu\in (0,1].\) The Allen–Cahn equation is an \(L^2\)-gradient flow of the energy functional \[\begin{align} \label{eq:energy} E(u):=\int_{ \mathbb{T}^2 }\frac{1}{2}\nu \|\nabla u\|_2^2+\frac{1}{4}(u^2-1)^2 \, dx, \end{align}\tag{17}\] which implies that \(E(u(t))\leq E(u(0))\) for all \(t>0\) [41]. Since for all \(u\in \mathbb{R}\), \(u^2\leq(u^2-1)^2+\frac{5}{4},\) we have that \(\|u\|_{L^2}^2\leq \int_{\mathbb{T}^2}(u^2-1)^2 \, dx+5\pi^2\leq 4E(u)+5\pi^2\). Consequently, \[\begin{align} \sup_{u: \|u\|_{H^4}\leq U}\|\mathcal{S}^{\mathrm{AC}}(u)\|_{L^2}\leq \sqrt{4E(u)+5\pi^2}\leq \sqrt{ 2\nu U^2+r^4U^4+9\pi^2}:=V, \end{align}\] where the constant \(r>0\) is from the Sobolev embedding \(\|u\|_{L^{\infty}}\leq r\|u\|_{H^4},\) which holds since \(a=4> 1=\frac{d}{2}.\) We consider the stabilized semi-implicit Euler scheme \[\begin{align} \label{eq:AC32Euler} \frac{1}{\tau}(u_N^{j+1}-u_N^{j}) =\nu \Delta u_N^{j+1} -\gamma(u_N^{j+1}-u_N^j) - \mathcal{P}_N\Bigl((u_N^j)^3-u_N^j\Bigr) \end{align}\tag{18}\] for \(\gamma= r(U^2+\nu^{-1}|\log(\nu)|^2)\), where \(\mathcal{P}_N\) is the \(L^2\)-orthogonal projection onto \(L^2_N(\mathbb{T}^2 ,\mathbb{R})\). This scheme can be written in the form of 8 with

1. \(\mathcal{B}^{\mathrm{AC}}:L^2(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) is given by \[\begin{align} \mathcal{B}^{\mathrm{AC}}&=\Bigl(\frac{1}{\tau} I+\gamma I+ \mathcal{A}^{\mathrm{AC}}\Bigr)^{-1}\mathcal{P}_N, \end{align}\] which implies \(\eta=0\), where

   1. \(\mathcal{P}_N:L^2(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) is given by \[\begin{align} (\mathcal{P}_Nu)(x)=\sum_{k\in \mathsf{K}_N}\hat{u}(k)e^{i\langle k,x\rangle}, \end{align}\] for any \(u\in L^2(\mathbb{T}^2,\mathbb{R});\)

   2. \(\mathcal{A}^{\mathrm{AC}}:L^2_N(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) is given by \[\begin{align} (\mathcal{A}^{\mathrm{AC}}u_N)(x)=(-\nu\Delta u_N)(x)=\nu\sum_{k\in \mathsf{K}_N}\|k\|_2^2\hat{u}_N(k)e^{i\langle k,x\rangle}, \end{align}\] for any \(u_N\in L^2_N(\mathbb{T}^2,\mathbb{R})\), which implies \(\beta=2\) and \(c=\nu.\)

   Consequently, \(\mathcal{B}^{\mathrm{AC}}\in \mathsf{B}(2,r(U^2+\nu^{-1}|\log(\nu)|^2),0,\nu).\)

2. \(\mathcal{D}_s^{\mathrm{AC}}:L^2_{3N}(\mathbb{T}^2,\mathbb{R})\to L^2_{3N}(\mathbb{T}^2,\mathbb{R})\) by \[\begin{align} \mathcal{D}_s^{\mathrm{AC}}(v_{3N})= v_{3N}, \end{align}\] for any \(v_{3N}\in L^2_{3N}(\mathbb{T}^2,\mathbb{R}^3),\) which implies \(s=0\) and \(C_{\mathcal{D}}=1\), so \(\mathcal{D}_s^{\mathrm{AC}}\in \mathsf{D}(0,1).\)

3. \(\mathcal{G}^{\mathrm{AC}}:L^2_{N}(\mathbb{T}^2,\mathbb{R})\to L^2_{3N}(\mathbb{T}^2,\mathbb{R})\) given by \[\begin{align} \mathcal{G}^{\mathrm{AC}}(u_N)= u_N^3-u_N, \end{align}\] for any \(u_N\in L^2_{N}(\mathbb{T}^2,\mathbb{R}),\) which implies \(p=3\) and \(C_{\mathcal{G}}=1\). Consequently, \(\mathcal{G}^{\mathrm{AC}}\in \mathsf{G}(3,1)\).

It is shown in [41] that with this choice of \(\gamma\), this method is energy-stable, meaning that for any \(j\geq 1,\) \[\begin{align} E\Bigl(\bigl(\mathcal{S}_{\tau,N}^{\mathrm{AC}}\bigr)^j(\mathcal{I}_Nu_0)\Bigr)\leq E\Bigl(\bigl(\mathcal{S}_{\tau,N}^{\mathrm{AC}}\bigr)^{j-1}(\mathcal{I}_Nu_0)\Bigr), \end{align}\] for \(\|u_0\|_{H^{4}}\leq U\) if \(\tau\) is chosen small enough relative to a constant depending only on \(U\), \(\gamma\), and \(\nu\). This implies a uniform stability bound in \(L^2\), \[\begin{align} \|\bigl(\mathcal{S}_{\tau,N}^{\mathrm{AC}}\bigr)^j(\mathcal{I}_N u_0)\|_{L^2} \leq \sqrt{4E(\mathcal{I}_Nu_0)+5\pi^2}\leq \sqrt{ 2\nu U^2+r^4U^4+9\pi^2}=V, \end{align}\] verifying the first part of Condition 1 (II). Further, the iterates have error bounded by \[\begin{align} \|u(j\tau)-\bigl(\mathcal{S}_{\tau,N}^{\mathrm{AC}}\bigr)^j(\mathcal{I}_Nu_0)\|_{L^2}\leq C_1 e^{C_2j\tau}\Bigl(N^{-4}+\tau\Bigr), \end{align}\] for \(\|u_0\|_{H^{4}}\leq U\), where \(C_1\) and \(C_2\) also only depend on \(\nu\), \(\gamma\), and \(U\). This implies that \(a=4,\) and verifies the second part of Condition 1 (II). Consequently, the Allen–Cahn solution operator belongs to \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\), and therefore satisfies the approximation and learning guarantees established in Theorem 3. We thus obtain the following corollary:

The rate \(M^{-1/7}\) is determined by the spatial smoothness parameter and the dimension of the physical space, which are \(a=4\) and \(d=2,\) respectively. The constant \(r\) is increasing in \(\nu^{-1}\) and \(U.\)

3.4.3 Cahn–Hilliard Equation↩︎

The Cahn–Hilliard equation is a mass-conserving phase-field model describing phase separation and coarsening in binary mixtures through fourth-order diffusion. On the two-dimensional torus \(\mathbb{T}^2\) and with periodic boundary conditions, it is given by \[\begin{align} \frac{du}{dt}- \Delta \bigl((- \nu\Delta) u+(u^3-u)\bigr)=0, \end{align}\] where \(u(x,t)\) is a real-valued function. The biharmonic operator induces fourth-order dissipation corresponding to \(\beta=4,\) while the cubic potential produces the degree-three polynomial nonlinearity governing phase separation dynamics. In contrast to the Allen–Cahn equation, whose dynamics are governed by second-order diffusion, the Cahn–Hilliard equation involves fourth-order diffusion, making stabilization particularly important for large-timestep numerical schemes.

We let \(\mathcal{S}^{\mathrm{CH}}:L^2(\mathbb{T}^2,\mathbb{R})\to L^2(\mathbb{T}^2,\mathbb{R})\) denote the time-\(T\) solution map of the Cahn–Hilliard equation, mapping \(u_0\mapsto u(T)=\mathcal{S}^{\text{CH}}(u_0).\) We now verify that \(\mathcal{S}^{\mathrm{CH}}\in \mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) for \[\vartheta=(\beta,\gamma,\eta,c, s,C_{\mathcal{D}},p, C_{\mathcal{G}}, F)=(4,r(U^2+\nu^{-1}|\log(\nu)|^2),\frac{1}{2},\nu,2,1,3,1,0)\] and \[\bar{\vartheta}=(\vartheta,a,U,V)=(\vartheta,4,U,\sqrt{ 2\nu U^2+r^4U^4+9\pi^2})\] for any \(U>0\) and \(\nu\in (0,1].\) By the energy dissipation estimate in [37], if \(u_0\in H^4(\mathbb{T}^2)\), then the solution satisfies \[\|\mathcal{S}^{\mathrm{CH}}(u_0)\|_{L^2} \leq \sqrt{4E(u_0)+5\pi^2}.\] Consequently, \[\sup_{\|u_0\|_{H^4}\leq U} \|\mathcal{S}^{\mathrm{CH}}(u_0)\|_{L^2} \leq \sqrt{2\nu U^2+r^4U^4+9\pi^2} =:V.\] This verifies Condition 1 (I). We consider the stabilized semi-implicit Euler scheme \[\begin{align} \label{eq:CH32euler} \frac{u_N^{j+1}-u_N^j}{\tau} =-\nu \Delta^2 u^{j+1}_N +\gamma\sqrt{\nu}\,\Delta (u_N^{j+1}-u_N^j) +\Delta \mathcal{P}_N\Bigl((u_N^j)^3-u_N^j\Bigr), \end{align}\tag{19}\] for \(\gamma\geq r(U^2+\nu^{-1}|\log(\nu)|^2)\), where \(\mathcal{P}_N\) is the \(L^2\)-orthogonal projector onto the space of mean-zero degree-\(N\) trigonometric polynomials. This scheme can be written in the form of 8 with

1. \(\mathcal{B}^{\mathrm{CH}}:L^2(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) given by \[\begin{align} \mathcal{B}^{\text{CH}}&=\Bigl(\frac{1}{\tau} I+\gamma \bigl(\mathcal{A}^{\mathrm{CH}}\bigr)^{1/2}+ \mathcal{A}^{\text{CH}}\Bigr)^{-1}\mathcal{P}_N, \end{align}\] which implies \(\eta=1/2\), where

   1. \(\mathcal{P}_N:L^2(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) is given by \[\begin{align} (\mathcal{P}_Nu)(x)=\sum_{k\in \mathsf{K}_N\backslash\{0\}}\hat{u}(k)e^{i\langle k,x\rangle}, \end{align}\] for any \(u\in L^2(\mathbb{T}^2,\mathbb{R});\)

   2. \(\mathcal{A}^{\mathrm{CH}}:L^2_N(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) is given by \[\begin{align} (\mathcal{A}^{\mathrm{CH}}u_N)(x)=(\nu\Delta^2 u_N)(x)=\nu \sum_{ k\in \mathsf{K}_N \backslash\{0\} }\|k\|_2^4\hat{u}_N(k)e^{i\langle k,x\rangle}, \end{align}\] for any \(u_N\in L^2_N(\mathbb{T}^2,\mathbb{R})\), which implies \(\beta=4\) and \(c=\nu.\)

   Consequently, \(\mathcal{B}^{\mathrm{CH}}\in \mathsf{B}(4,r(U^2+\nu^{-1}|\log(\nu)|^2),\frac{1}{2},\nu)\).

2. \(\mathcal{D}_s^{\mathrm{CH}}:L^2_{3N}(\mathbb{T}^2,\mathbb{R})\to L^2_{3N}(\mathbb{T}^2,\mathbb{R})\) given by \[\begin{align} \mathcal{D}_s^{\mathrm{CH}}(v_{3N})=-\Delta v_{3N}=-\sum_{j=1}^2\partial^2_{x_j}v_{3N}, \end{align}\] for any \(v_{3N}\in L^2_{3N}(\mathbb{T}^2,\mathbb{R}),\) which implies \(s=2\) and \(C_{\mathcal{D}}=1\), so \(\mathcal{D}_s^{\mathrm{CH}}\in \mathsf{D}(2,1).\)

3. \(\mathcal{G}^{\mathrm{CH}}:L^2_{N}(\mathbb{T}^2,\mathbb{R})\to L^2_{3N}(\mathbb{T}^2,\mathbb{R})\) given by \[\begin{align} \mathcal{G}^{\mathrm{CH}}(u_N)=-u_N+u_N^3, \end{align}\] for any \(u_N\in L^2_{N}(\mathbb{T}^2,\mathbb{R}),\) which implies \(p=3\) and \(C_{\mathcal{G}}=1\). Consequently, \(\mathcal{G}^{\mathrm{CH}}\in \mathsf{G}(3,1)\).

It is shown in [41] that with this choice of \(\gamma\), the scheme is energy-stable, so for any \(j\geq 1,\) \[\begin{align} E\Bigl( \bigl(\mathcal{S}_{\tau,N}^{\text{CH}}\bigr)^j(\mathcal{I}_Nu_0)\Bigr)\leq E\Bigl( \bigl(\mathcal{S}_{\tau,N}^{\text{CH}}\bigr)^{j-1}(\mathcal{I}_Nu_0)\Bigr), \end{align}\] for \(\|u_0\|_{H^{4}}\leq U\) if \(\tau\) is chosen small enough relative to a constant depending only on \(U\), \(\gamma\), and \(\nu\). Here \(E\) is the energy functional defined in 17 . This implies a uniform stability bound in \(L^2\), \[\begin{align} \|\bigl(\mathcal{S}_{\tau,N}^{\mathrm{CH}}\bigr)^j( \mathcal{I}_N \normalcolor u_0 )\|_{L^2}\leq \sqrt{4E(\mathcal{I}_Nu_0)+5\pi^2}\leq \sqrt{ 2\nu U^2+r^4U^4+9\pi^2}=V, \end{align}\] verifying the first part of Condition 1 (II). As shown in [47], this scheme admits an error bound given by \[\begin{align} \|u(j\tau)-\bigl(\mathcal{S}_{\tau,N}^{\text{CH}}\bigr)^j(\mathcal{I}_Nu_0)\|_{L^2}\leq r_1e^{r_2 j \tau}\Bigl(N^{-4}+\tau\Bigr) \end{align}\] for \(\|u_0\|_{H^4}\leq U\) with \(\tau\) small enough. This implies that \(a=4,\) and verifies the second part of Condition 1 (II). Consequently, the Cahn–Hilliard solution operator belongs to \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\), and therefore satisfies the approximation and learning guarantees established in Theorem 3. We thus obtain the following corollary:

As in the previous examples, the rate \(M^{-1/7}\) is determined by the spatial smoothness parameter and the dimension of the physical space, which are \(a=4\) and \(d=2\), respectively. The constant \(r\) is increasing in \(\nu^{-1}\) and \(U.\)

4 Dissipative Models with Smooth Nonlinearity↩︎

In this section, we develop a theory for approximating and learning the time-\(T\) solution operator for evolution equations of the form \[\begin{align} \label{eq:evolution32PDEsmooth} \begin{aligned} \frac{du}{dt}+\mathcal{A}u+\mathcal{D}_s\mathcal{G}( u) &=f, \qquad t>0, \end{aligned} \end{align}\tag{20}\] using FNOs under non-parametric smoothness assumptions on the nonlinearity \(\mathcal{G}.\) This contrasts with Section 3, where \(\mathcal{G}\) is assumed to have polynomial structure. Under these general smoothness assumptions, we still obtain polynomial sample complexity for learning the time-\(T\) solution operator; however, the rate is slower than in the polynomial case. In direct analogue with Section 3, we begin by defining a family of spectral methods and then define a corresponding class of operators that can be approximated by these methods. For this operator class, which includes more general reaction-diffusion equations, we establish FNO approximation bounds in Section 4.2 and learning bounds in Section 4.3. In Section 4.4, we apply this theory to show that reaction-diffusion equations can be efficiently learned by FNOs.

4.1 Spectral Methods and Class of Operators↩︎

4.1.1 Background on Spectral Methods↩︎

Following the ideas in Subsection 3.1.1, we consider first-order-in-time semi-implicit spectral methods of the form \[\begin{align} \label{eq:generic32spectral32method32smooth} \frac{1}{\tau}\bigl(u_N^{j+1}-u_N^j\bigr)+\gamma\mathcal{A}^{\eta}\bigl(u_N^{j+1}-u_N^j\bigr)+\mathcal{A}u_N^{j+1}+ \mathcal{P}_N\mathcal{D}_s\bigl(\mathcal{G}(u_N^j)\bigr)-f_N=0, \end{align}\tag{21}\] where \(\gamma\geq0\) and \(\eta\in [0,1)\) control the strength of the stabilization term. As in the previous section, the mapping \(u_N^{j}\mapsto u_N^{j+1}\) is uniquely defined, and we write a single step of the spectral method \(\mathcal{S}_{\tau,N}:L^2_N(\mathbb{T}^d,\mathbb{R}^{d'})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{d'})\) as \[\begin{align} \label{eq:spec32smooth} \begin{aligned} \mathcal{S}_{\tau,N}(u_N) :=& \, \, \Bigl(\frac{1}{\tau}I+\gamma\mathcal{A}^{\eta}+\mathcal{A}\Bigr)^{-1}\Bigl((\frac{1}{\tau}I+\gamma\mathcal{A}^{\eta})u_N-\mathcal{P}_N\mathcal{D}_s\bigl(\mathcal{G}(u_N)\bigr)+f_N\Bigr)\\ =& \, \, \mathcal{B}\Bigl(\bigl(\frac{1}{\tau}I+\gamma \mathcal{A}^{\eta}\bigr)u_N- \mathcal{D}_s\bigl(\mathcal{G}(u_N)\bigr)+f_N\Bigr), \end{aligned} \end{align}\tag{22}\] where \(\mathcal{B}:= \Bigl(\frac{1}{\tau}I +\gamma\mathcal{A}^{\eta}+ \mathcal{A}\Bigr)^{-1} \mathcal{P}_N.\) Once again, to approximate the time-\(T\) solution operator of the PDE, one can choose \(J \in \mathbb{N}\) and \(\tau>0\) such that \(J=T/\tau,\) and then iterate the single-step spectral method: \[\begin{align} u(T) = \mathcal{S}(u) \approx \mathcal{S}_{\tau,N}^J(\mathcal{P}_Nu(0))=\underbrace{\mathcal{S}_{\tau,N}\circ \cdots \circ \mathcal{S}_{\tau,N}}_{J \text{ times}}(u_N^0). \end{align}\] In the next subsection, we introduce a family of single-step spectral methods with general smooth nonlinearities. Then, in the following subsection, we introduce a class of operators that can be well approximated by these single-step spectral methods.

4.1.2 Family of Single-Step Spectral Methods↩︎

Recall the generic single-step spectral method \(\mathcal{S}_{\tau,N}\) in 22 . For timestep \(\tau>0\) and degree \(N \in \mathbb{N},\) we introduce the parameter \(\vartheta:=( \beta,\gamma,\eta,c,s,C_{\mathcal{D}},\alpha, C_{\mathcal{G}}, F)\) and define the family of single-step spectral methods

\[\begin{align} \mathsf{S}_{\tau,N}^{\text{smooth}}(\vartheta) :=\Biggl\{\mathcal{S}_{\tau,N}:L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\to L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})& \, \bigg| \, \mathcal{S}_{\tau,N}(u_N)=\mathcal{B}\Bigl(\bigl(\frac{1}{\tau}I+\gamma\mathcal{A}^{\eta} \bigr)u_N-\mathcal{D}_s\mathcal{G}(u_N)+f_N\Bigr), \\ &\mathcal{B}\in \mathsf{B}(\beta,\gamma,\eta,c), \, \mathcal{D}_s\in \mathsf{D}(s,C_{\mathcal{D}}), \,\mathcal{G}\in \mathsf{G}(\alpha,C_{\mathcal{G}}), \, f_N\in \mathsf{F}(F) \Biggr\}. \end{align}\] Here, the classes \(\mathsf{B}(\beta,\gamma,\eta,c)\), \(\mathsf{D}(s,C_{\mathcal{D}})\), and \(\mathsf{F}(F)\) are exactly as defined in Subsection 3.1.2, while \(\mathsf{G}(\alpha,C_{\mathcal{G}})\) now denotes the class of smooth nonlinearities (applied component-wise and point-wise) with \(\alpha\) derivatives bounded by \(C_{\mathcal{G}}.\) That is, \[\begin{align} \mathsf{G}(\alpha,C_{\mathcal{G}})&:=\Biggl\{\mathcal{G}:L^2_{N}(\mathbb{T}^d,\mathbb{R}^{{d_*}})\to H^{\alpha-1}(\mathbb{T}^d,\mathbb{R}^{\tilde{d}}) \, \bigg\vert \\\, & \mathcal{G}(v_N)(x)=\begin{bmatrix} g_{1}(v_N(x))\\ \vdots \\ g_{\tilde{d}}(v_N(x)) \end{bmatrix} \,\,\, \forall x\in \mathbb{T}^d, \, g_{j}\in W^{\alpha,\infty}(C_{\mathcal{G}}), 1 \le j \le \tilde{d}\Biggr\}, \end{align}\] where \[\begin{align} W^{\alpha,\infty}(C_{\mathcal{G}}):=\biggl\{g\in L^{\infty}(\mathbb{R}^{{d_*}}):\max_{\omega :|\omega|\leq \alpha}\|\partial^{\omega} g\|_{L^{\infty}(\mathbb{R}^{{d_*}})}\leq C_{\mathcal{G}}\biggr\}. \end{align}\] For simplicity, we assume throughout this section that \(\mathcal{P}_N\) projects onto a subspace of \(L^2_N(\mathbb{T}^d,\mathbb{R}^{{d_*}})\) consisting only of mean-zero functions, so that \[\begin{align} (\mathcal{P}_Nu)(x)=\sum_{k\in \mathsf{K}_N \backslash \{0\}}\hat{P}_{\mathcal{P}_N}(k)\hat{u}(k)e^{i\langle k,x\rangle}, \end{align}\] where \(\hat{P}_{\mathcal{P}_N}(k)\in \mathbb{C}^{{d_*}\times {d_*}}\) are orthogonal projections satisfying \(\hat{P}_{\mathcal{P}_N}(k)=\mkern 1.5mu\overline{\mkern-1.5mu\hat{P}_{\mathcal{P}_N}(-k)\mkern-1.5mu}\mkern 1.5mu\) for all \(k\in \mathsf{K}_N\backslash \{0\},\) and \(\hat{P}_{\mathcal{P}_N}(0)=0.\) This implies that for any \(u_N\in \text{Range}(\mathcal{P}_N)\)—in particular, for any iterate of the spectral method, \(\mathcal{S}_{\tau,N}^j(u_N)\)—we have the Poincaré inequality \[\begin{align} c\|u_N\|_{L^2}\leq \|\mathcal{A}u_N\|_{L^2}. \end{align}\] We will assume throughout that \(\alpha>\frac{d}{2}\) and \(s\leq\alpha-1,\) so that \(\mathcal{D}_s\mathcal{G}(u_N)\) is defined as a classical derivative.

4.1.3 Class of Operators↩︎

Next, we introduce the class of operators that can be approximated by spectral methods in \(\mathsf{S}_{\tau,N}^{\text{smooth}}(\vartheta)\), defined analogously to the class considered in Subsection 3.1.3. We consider operators acting on a Sobolev ball of radius \(U,\) denoted by \(\mathsf{U}_a=\{u\in L^2(\mathbb{T}^d,\mathbb{R}^{{d_*}}):\|u\|_{H^a}\leq U\}\). For fixed \(T>0\) and parameter \(\bar{\vartheta}:=(\vartheta,a,U,V),\) with \(a>d/2,\) we define \(\mathsf{S}^{\text{smooth}}(\bar{\vartheta})\) to be the class of operators \(\mathcal{S}\) for which the following holds:

As in Subsection 3.1.3, we assume that the dynamics are bounded, which is natural for dissipative equations, and that the solution operator \(\mathcal{S}\) can be approximated by a stable and convergent spectral method, now belonging to \(\mathsf{S}_{\tau,N}^{\text{smooth}}(\vartheta)\). The smoothness parameter \(a\) again plays a fundamental role in the approximation and learning rates in the following sections. However, we will see shortly that, under the weaker assumptions on the nonlinearity considered here, the rates also inherit a dependence on the smoothness of the nonlinearity, \(\alpha\), and the rate of dissipation, \(\beta\).

In general, polynomial nonlinearities do not have uniformly bounded derivatives, and thus the class \(\mathsf{S}^{\mathrm{poly}}\) is not a subset of \(\mathsf{S}^{\text{smooth}}.\) However, if the PDE solution and spectral method iterates satisfy an a priori \(L^{\infty}\) stability bound, then polynomial nonlinearities restricted to these bounded inputs have uniformly bounded derivatives of all orders and are contained in the class considered in this section. This is the case for the Allen–Cahn example in Section 3.4.

4.2 Approximation Bound↩︎

We now state the first main result of this section: a quantitative estimate on the FNO size required to approximate any operator in \(\mathsf{S}^{\text{smooth}}(\bar{\vartheta})\) to any desired accuracy. We show that for any \(\mathcal{S}\in \mathsf{S}^{\text{smooth}}(\bar{\vartheta})\), there exists an FNO in the following class that can approximate \(\mathcal{S}\) to accuracy \(\varepsilon\): \[\begin{align} \label{eq:sigmaclass32II} \begin{aligned} \Sigma^{\text{smooth}}(\varepsilon,\bar{\vartheta})=\biggl\{\Psi\in \mathsf{FNO}(\mathscr{D},\mathscr{W},\mathscr{S},\mathscr{C}): \mathscr{D}&=r \varepsilon^{-1}\log(\varepsilon^{-1}),\\ \mathscr{W}&=r\varepsilon^{-\frac{{d_*}d}{2a}-\frac{{d_*}\max\{d+s-\beta/2,0\}}{\alpha a}-\frac{{d_*}}{\alpha}}{\log(\varepsilon^{-1})^{\frac{{d_*}}{\alpha}}},\\ \mathscr{S}&=r\varepsilon^{-\frac{\max\{d+s-\beta/2,0\}}{\alpha a}-\frac{1+s}{\alpha}-\frac{(d+2)s}{2a}}{\log(\varepsilon^{-1})^{\frac{1}{\alpha}}}\, ,\\ \mathscr{C}&=r\varepsilon^{-\frac{1}{\alpha}-\frac{2+d}{2a}},\\ &\sup_{u\in \mathsf{U}_a}\|\Psi(u)\|_{L^2}\leq 2V\biggr\} \, , \end{aligned} \end{align}\tag{23}\] where \(r\) is a positive constant increasing in \(\beta^{-1}\), \(\gamma\), \(\eta\), \(c^{-1},\) \(\alpha,\) \(C_{\mathcal{G}}\), \(F\), \(U\), and \(V,\) which can be traced through the proof of Theorem 5. Precisely, we have the following approximation result.

Theorem 5 provides a uniform approximation guarantee over the class of operators \(\mathsf{S}^{\text{smooth}}(\bar{\vartheta})\) with general smooth nonlinearities, and the required size of the FNO architecture is controlled quantitatively in terms of the desired accuracy \(\varepsilon\), as specified in the definition of \(\Sigma^{\text{smooth}}(\varepsilon,\bar{\vartheta})\). These quantitative estimates on the approximation error will be essential to our learning-theoretic result in Section 4.3.

The proof of Theorem 5 is presented in Section 9.1. The construction of the FNO achieving \(\varepsilon\) accuracy mirrors that in Theorem 2; however, the proof of the error bound differs in several substantial ways. The ReLU neural network approximation rates for functions in \(W^{\alpha,\infty}(C_{\mathcal{G}})\) are slower than those for polynomials, and these bounds cannot be improved in general [48]. As such, we require a more delicate Gronwall argument than in the proof of Theorem 2. Further, it is much more difficult to control the aliasing error incurred by approximating general smooth functions of trigonometric polynomials than it is to control polynomial functions of trigonometric polynomials.

4.3 Learning Bound↩︎

In this section, we show that operators \(\mathcal{S}\in\mathsf{S}^{\text{smooth}}(\bar{\vartheta})\) can be efficiently learned by FNOs trained via empirical risk minimization. As in Section 3.3, we assume access to training data of the form \(\{u^{(m)},\mathcal{S}(u^{(m)})\}_{m=1}^M,\) where the input functions \(u^{(1)},\ldots,u^{(M)}\) are sampled independently from a measure \(\mu\) with \(\text{supp}(\mu)\subseteq \mathsf{U}_a\). As before, we consider the ERM estimator \[\begin{align} \label{eq:empirical32risk32II} \begin{aligned} \hat{\Psi}_{\mathcal{S}}&\in \mathop{\mathrm{arg\,min}}_{ \Psi \in \Sigma^{\mathrm{smooth}}(\varepsilon,\bar{\vartheta})}\hat{\mathcal{L}}(\Psi,\mathcal{S}), \\ \hat{\mathcal{L}}(\Psi,\mathcal{S})&:=\frac{1}{M}\sum_{m=1}^{M}\|\Psi( u^{(m)} )-\mathcal{S}( u^{(m)} )\|_{L^2}^2, \end{aligned} \end{align}\tag{24}\] for a choice of \(\varepsilon\) specified in Theorem 6. We measure error using the population risk \[\begin{align} \label{eq:population32risk32II} \mathcal{L}(\Psi,\mathcal{S}):=\mathbb{E}_{u\sim \mu}\|\Psi(u)-\mathcal{S}(u)\|_{L^2}^2. \end{align}\tag{25}\] The main result of this section is the following:

Theorem 6 demonstrates that the operator class \(\mathsf{S}^{\text{smooth}}(\bar{\vartheta})\) can be learned with polynomial sample complexity using FNOs. As before, we contrast this result with the hardness result in [21], which shows that generic Fréchet differentiable or Lipschitz operators cannot be learned from only polynomially many samples. The structure imposed in \(\mathsf{S}^{\text{smooth}}(\bar{\vartheta})\) is enough to overcome this hardness, even though the spectral methods that approximate our target operators make only generic smoothness assumptions. The proof of Theorem 6 is given in Section 9.2 and proceeds similarly to the proof of Theorem 3.

4.4 Example: Learning the Solution Map of a Nonlinear Gradient Flow↩︎

We consider the Cahn–Hilliard equation with a logarithmic Flory–Huggins potential on the two-dimensional periodic torus \(\mathbb{T}^2\), \[\begin{align} \label{eq:FH32PDE} \frac{\partial u}{\partial t}- \Delta \bigl(-\nu \Delta u +\mathcal{W}'(u)\bigr)=0, \end{align}\tag{26}\] where \(\mathcal{W}:(-1,1)\to \mathbb{R}\) is referred to as the thermodynamic potential or homogeneous free energy, given by \[\begin{align} \mathcal{W}(u)=\frac{\theta}{2}\Bigl((1+u)\log(1+u)+(1-u)\log(1-u)\Bigr)-\frac{\theta_c}{2}u^2. \end{align}\] Here \(\theta\) represents the absolute temperature, and \(\theta_c\) represents the critical temperature. Under the standard condition that \(0<\theta<\theta_c\), the potential \(\mathcal{W}\) has a double-well shape and phase separation can occur. The derivative of \(\mathcal{W}\) is given by \[\begin{align} \mathcal{W}'(u)=w(u)=\frac{\theta}{2}\log((1+u)/(1-u))-\theta_cu. \end{align}\] The logarithmic potential considered here arises from first principles to model the entropy of mixing [49]–[51]. The usual Cahn–Hilliard equation, as considered in Section 3.4, can be derived as a quartic approximation to \(\mathcal{W}\), \[\begin{align} \mathcal{W}(u)\approx \mathcal{W}_{\text{quartic}}(u)=\frac{\theta}{2}\frac{u^4}{6}+\Bigl(\frac{\theta}{2}-\frac{\theta_c}{2}\Bigr)u^2, \end{align}\] with the specific choice of \(\frac{\theta}{\theta_c}=\frac{3}{4}\). From the perspective of physical modeling, the logarithmic potential considered here is more natural despite presenting additional mathematical and computational challenges due to the possibility of singularities in the potential. We let \(\mathcal{S}^{\mathrm{FH}}:L^2(\mathbb{T}^2,\mathbb{R})\to L^2(\mathbb{T}^2,\mathbb{R})\) denote the time-\(T\) solution map of the Cahn–Hilliard equation with Flory–Huggins potential, mapping \(u_0\mapsto u(T)=\mathcal{S}^{\text{FH}}(u_0).\) This map is only well-defined and physically meaningful for initial conditions satisfying \(\|u_0\|_{L^{\infty}}\leq 1-\delta_0,\) for \(\delta_0\in (0,1]\). As such, we restrict our attention to such inputs. We now verify that \(\mathcal{S}^{\mathrm{FH}}\in \mathsf{S}^{\text{smooth}}(\bar{\vartheta})\) for \[\vartheta=(\beta,\gamma,\eta,c, s,C_{\mathcal{D}},\alpha, {C_{\mathcal{G}}}, F)=\Bigl(4,0,0,\nu,2,1,\alpha,C_{\alpha }(\theta_c + \theta(\alpha-1)!\delta_1^{-\alpha}),0 \Bigr)\] and \(\bar{\vartheta}=(\vartheta,a,U,V)=(\vartheta,5,U,2\pi(1-\delta_1))\) for any \(U,\delta_0,\theta,\theta_c>0\) such that \(0<\theta<\theta_c,\) \(\nu\in (0,1]\), and \(\alpha\geq 3.\) First we verify boundedness. From [38], there exists \(\delta_1\in (0,1)\) depending on \(U, \theta,\theta_c,\) and \(\delta_0\) such that \[\begin{align} \sup_{u:\|u\|_{H^5}\leq U}\|\mathcal{S}^{\mathrm{FH}}(u)\|_{L^{\infty}}\leq 1-\delta_1, \end{align}\] implying that \[\begin{align} \sup_{u:\|u\|_{H^5}\leq U}\|\mathcal{S}^{\mathrm{FH}}(u)\|_{L^{2}}\leq 2\pi(1-\delta_1):=V. \end{align}\] This verifies Condition 4 (I). We consider the semi-implicit Euler scheme \[\begin{align} \label{eq:FH32scheme} \frac{u^{j+1}_N-u^{j}_N}{\tau}=-\nu \Delta^2u^{j+1}_N-\theta_c\Delta u_N^{j}+\Delta \mathcal{P}_N\frac{\theta}{2}\log\biggl(\frac{1+u_N^j}{1-u_N^j}\biggr), \end{align}\tag{27}\] where \(\mathcal{P}_N:L^2(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) is the \(L^2\)-orthogonal projection onto the space of mean-zero degree-\(N\) trigonometric polynomials. Theorem 8 demonstrates that the iterates given by 27 remain bounded pointwise by a constant strictly less than \(1\), such that \(\|u_N^j\|_{L^{\infty}}\leq 1-\delta_1\) for all \(j\geq 0.\) To ensure the nonlinearity is globally well-defined on the entire domain space, as is required to apply Theorem 6, we equivalently consider the extended scheme \[\begin{align} \label{eq:FH32scheme322} \frac{u^{j+1}_N-u^{j}_N}{\tau}=-\nu \Delta^2u^{j+1}_N+\Delta \mathcal{P}_N \tilde{w}(u_N^j), \end{align}\tag{28}\] where \(\tilde{w}:\mathbb{R}\rightarrow\mathbb{R}\) is a smooth extension of \(w:{[-1+\delta_1, 1-\delta_1]}\to \mathbb{R}\) to the entire real line. By classical extension results [52], \(\tilde{w}\) can be constructed to satisfy \(\tilde{w}(u) = w(u)\) for all \(u\in[-1+\delta_1, 1-\delta_1]\) and obey the uniform derivative bound \[\begin{align} \label{eq:extension32derivative32bound} \max_{\omega:|\omega|\leq \alpha}\|\partial^{\omega} \tilde{w}\|_{L^{\infty}(\mathbb{R})}\leq C_{\alpha}\max_{\omega:|\omega|\leq \alpha}\|\partial^{\omega} w\|_{L^{\infty}([-1+\delta_1,1-\delta_1])}, \end{align}\tag{29}\] where \(C_{\alpha}\) is a universal constant depending solely on the smoothness \(\alpha\). The iterates given by 27 and 28 coincide exactly provided \(\|u_N^0\|_{L^{\infty}}\leq 1-\delta_0\) for \(\delta_0<1\) and \(u_N^0\in \mathsf{U}_5.\) We now demonstrate that the scheme 28 can be written in the form of 22 with

1. \(\mathcal{B}^{\mathrm{FH}}:L^2(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) given by \[\begin{align} \mathcal{B}^{\text{FH}}&=\Bigl(\frac{1}{\tau} I+ \mathcal{A}^{\text{FH}}\Bigr)^{-1}\mathcal{P}_N, \end{align}\] which corresponds to \(\gamma=\eta=0\), where

   1. \(\mathcal{P}_N:L^2(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) is given by \[\begin{align} (\mathcal{P}_Nu)(x)=\sum_{ k\in \mathsf{K}_N\backslash\{0\}}\hat{u}(k)e^{i\langle k,x\rangle}, \end{align}\] for any \(u\in L^2(\mathbb{T}^2,\mathbb{R}).\)

   2. \(\mathcal{A}^{\mathrm{FH}}:L^2_N(\mathbb{T}^2,\mathbb{R})\to L^2_N(\mathbb{T}^2,\mathbb{R})\) is given by \[\begin{align} (\mathcal{A}^{\mathrm{FH}}u_N)(x)=(\nu\Delta^2 u_N)(x)=\nu\sum_{ k\in \mathsf{K}_N\backslash\{0\}}\|k\|_2^4\:\hat{u}_N(k)e^{i\langle k,x\rangle}, \end{align}\] for any \(u_N\in L^2_N(\mathbb{T}^2,\mathbb{R})\), which implies \(\beta=4\) and \(c=\nu.\)

   Consequently, \(\mathcal{B}^{\mathrm{FH}}\in \mathsf{B}(4,0,0,\nu)\).

2. \(\mathcal{D}_s^{\mathrm{FH}}:H^{\alpha-1}(\mathbb{T}^{2},\mathbb{R}) \to L^2(\mathbb{T}^2,\mathbb{R})\) given by \[\begin{align} \mathcal{D}_s^{\mathrm{FH}}(u)=-\Delta u=-\sum_{j=1}^2\partial_{x_j}^2u, \end{align}\] for any \(u\in H^{\alpha-1}(\mathbb{T}^2,\mathbb{R}),\) which implies \(s=2\) and \(C_{\mathcal{D}}=1\), so \(\mathcal{D}_s^{\mathrm{FH}}\in \mathsf{D}(2,1).\)

3. \(\mathcal{G}^{\mathrm{FH}}:L^2_N(\mathbb{T}^2,\mathbb{R})\to H^{\alpha-1}(\mathbb{T}^2,\mathbb{R})\) given by \[\begin{align} \mathcal{G}^{\mathrm{FH}}(u_N)=\tilde{w}(u_N), \end{align}\] where \(\tilde{w}:\mathbb{R}\to \mathbb{R}\) is the smooth global extension of \(w\) introduced above. Since \[\max_{\omega :|\omega|\leq \alpha}\|\partial^{\omega} w\|_{L^{\infty}([-1+\delta_1,1-\delta_1])}\leq\theta_c + \theta(\alpha-1)!\frac{1}{\delta_1^\alpha},\] it follows from 29 that \[\max_{\omega:|\omega|\leq \alpha}\|\partial^{\omega}\tilde{w}\|_{L^{\infty}}\leq C_{\alpha}\bigl(\theta_c + \theta(\alpha-1)!\frac{1}{\delta_1^\alpha}\bigr),\] and thus \(\mathcal{G}^{\mathrm{FH}}\in \mathsf{G}(\alpha, C_{\alpha }(\theta_c + \theta(\alpha-1)!\delta_{1}^{-\alpha}))\) for any \(\alpha\geq 3.\)

In Theorem 8, we show that for \(u_0\in \mathsf{U}_5\) with \(\|u_0\|_{L^\infty}\leq 1-\delta_0\), \[\sup_{j\geq 1} \mathopen{}\mathclose{\left\|(\mathcal{S}_{\tau,N}^{\mathrm{FH}})^j(\mathcal{I}_N u_0) }\right\|_{L^{\infty}}\leq 1-\delta_1,\] provided \(\tau\) is chosen small enough and \(N\) is chosen large enough relative to \(U,\delta_0,\nu,\theta,\theta_c.\) Thus \[\sup_{j\geq 1}\|(\mathcal{S}_{\tau,N}^{\mathrm{FH}})^j(\mathcal{I}_N u_0)\|_{L^{2}}\leq 2\pi(1-\delta_1)=V.\] This verifies the first part of Condition 4 (II). Then in Theorem 9, we show that this scheme admits an error bound given by \[\begin{align} \mathopen{}\mathclose{\left\|u(\tau j)-(\mathcal{S}_{\tau,N}^{\mathrm{FH}})^j(\mathcal{I}_Nu_0) }\right\|_{L^2}\leq r_1e^{r_2j\tau}\Bigl(N^{-5}+\tau\Bigr) \end{align}\] for \(\|u_0\|_{H^5}\leq U\) with \(\|u_0\|_{L^{\infty}}\leq 1-\delta_0\) for \(\tau\) appropriately small and \(N\) appropriately large relative to \(U, \delta_0, \nu, \theta, \theta_c\). This implies that \(a=5\) and verifies the second part of Condition 4 (II). Consequently, the logarithmic Cahn-Hilliard solution operator belongs to \(\mathsf{S}^{\text{smooth}}(\bar{\vartheta})\), and therefore satisfies the approximation and learning guarantees established in Theorem 6. We obtain the following corollary:

Here the rate \(M^{\frac{-1}{7+\frac{24}{5\alpha}}}\) is determined by the spatial smoothness parameter \(a=5\), the dimension of the physical space \(d=2\), and the smoothness of the nonlinearity \(\alpha\). While \(\alpha\) could be taken to be arbitrarily large to decrease the rate, the constant \(r\) rapidly blows up with \(\alpha\) as \(r\gtrsim \exp\Bigl(\frac{T}{\nu} (\theta_c + \theta(\alpha-1)!\frac{1}{\delta_1^\alpha})^2\Bigr)\) because of the singularities of the logarithmic potential at \(\pm1\). Consequently, taking \(\alpha\) to be too large quickly leads to a vacuous bound. Given a constrained number of samples, one can in principle optimize this bound over \(\alpha\); however we do not pursue this here.

5 Conclusions↩︎

We developed approximation and learning guarantees for FNOs applied to time-\(T\) solution operators of dissipative evolution equations. The main principle underlying our results is that solution operators admitting stable and accurate Fourier spectral discretizations can be efficiently represented by FNOs and learned from data with polynomial sample complexity.

We formalized this principle by introducing classes of evolution operators defined through spectral methods and proving quantitative FNO approximation and learning bounds for these classes. For polynomial nonlinearities, the rates depend on the Sobolev smoothness of the input space and the dimension of the physical domain. We verified that the framework applies to the Navier–Stokes, Allen–Cahn, and Cahn–Hilliard equations, and extended the theory to general smooth nonlinearities, where the rates also depend on the smoothness of the nonlinearity and the strength of dissipation.

These results identify stable spectral approximability as a mechanism leading to parameter-efficient FNO representations and data-efficient learning. Natural directions for future work include non-periodic domains, more general boundary conditions, and broader classes of numerical discretizations and operator-learning architectures.

Our learning guarantees for dissipative evolution equations are in part motivated by the use of data-driven surrogate forecasting models within data assimilation algorithms [53], [54]. Current theory suggests that to fully couple learning theoretic guarantees for surrogate models with accuracy guarantees for data assimilation, the learning bounds should be in a higher-order Sobolev space than \(L^2\) [55], [56]. Pursuing learning guarantees in higher-order norms is therefore a natural direction for future work.

Acknowledgments↩︎

DSA and NW were partly funded by NSF CAREER DMS-2237628. The authors are thankful to Jiaheng Chen for helpful feedback and discussions.

6 Background on Multilayer Perceptrons and FNOs↩︎

This appendix overviews existing approximation results for multilayer perceptrons and FNOs.

6.1 Multilayer Perceptrons↩︎

We characterize the size of an MLP through the following quantities:

• The depth \(\mathscr{D}(\Phi):=L\);

• The width \(\mathscr{W}(\Phi):= \max_{\ell=0,\ldots,L}d_{\ell}\);

• The weight magnitude \(\mathscr{S}(\Phi):=\max_{\ell=1,\ldots,L}\max \bigl\{\| {\boldsymbol{A}}_{\ell}\|_{\infty}, \|b_{\ell}\|_{\infty} \bigr\}\).

In the sequel, we restrict our attention to the rectified linear unit (ReLU) activation function, \(\sigma(z)=\max\{0,z\}.\) We let \(\mathcal{N}_{d_{\rm in},d_{\rm out}}\) denote the set of all ReLU networks with input dimension \(d_{\rm in}\) and output dimension \(d_{\rm out}\). The following lemmas from [57] will be useful for constructing such networks.

Consider the set of polynomials of \(n\) variables of degree at most \(p\), \[\begin{align} h(z)=c_0+\sum_{i_1}c_{i_1}z_{i_1}+\sum_{i_1,i_2}c_{i_1,i_2}z_{i_1}z_{i_2}+\cdots+\sum_{i_1,\ldots,i_p}c_{i_1,\ldots, i_p}z_{i_1}z_{i_2}\ldots z_{i_p}, \end{align}\] where the coefficient vector \(c=(c_0,c_1,\ldots, c_n, c_{1,1},\ldots, c_{n,n},\ldots)\) is of length \({p+n\choose n}\). In the following lemma, we quantify the size of a ReLU network sufficient to approximate a generic polynomial of \(n\) variables of degree at most \(p\).

For any domain \(\Omega \subseteq \mathbb{R}^d\), we define \[\begin{align} W^{\alpha,\infty}_{\Omega}(C):=\biggl\{f\in L^{\infty}(\Omega):\max_{\kappa :|\kappa|\leq \alpha}\|D^{\kappa} f\|_{L^{\infty}(\Omega)}\leq C\biggr\}. \end{align}\] When \(\Omega=\mathbb{R}^d\), we typically write \(W^{\alpha,\infty}_{\mathbb{R}^d}(C)=W^{\alpha,\infty}(C).\)

6.2 FNO Background Lemmas↩︎

We have the following FNO analogs of the preceding lemmas for constructing MLPs.

The following result [13] shows that pointwise MLPs can be realized as FNOs without increasing their depth, width, or weight magnitude.

7 FNO Approximation of the Building-Block Operators↩︎

This appendix contains FNO approximation bounds for the operators involved in our classes of single-step spectral methods.

In the following results, \(\mathcal{A}\) and \(\mathcal{P}_N\) are as in Subsection 3.1.2, and in particular, \(\mathcal{A}\) satisfies that for all \(k\in \mathsf{K}_N\backslash\{0\},\) \[\begin{align} 0 \prec \|k\|_{2}^{\beta}cI \preceq \hat{P}_{\mathcal{A}}(k)\preceq \|k\|_{2}^{\beta}I, \end{align}\] for some \(c\in (0,1].\)

In the following lemma, we construct an FNO to approximate the operator \(\mathcal{D}_s\bigl(\mathcal{G}(\cdot)\bigl)\). We quantify the mean-zero part of the approximation error in the \(L^2\) norm weighted by \(\mathcal{A}^{-1/2}\), which is most convenient for the multi-step approximation bound in Subsection 8.1 and allows for faster approximation by exploiting the dissipation due to the operator \(\mathcal{A}\). We also quantify the size of the mean in the \(L^2\) norm, identifying the mean with the constant function taking that value on \(\mathbb{T}^d\).

8 Supplementary Materials Section 3↩︎

8.1 Proof of Approximation Bound↩︎

8.2 Proof of Learning Bound↩︎

This section contains the proof of Theorem 3. Recall that we train FNOs by minimizing the empirical squared \(L^2\) prediction risk \[\begin{align} \hat{\mathcal{L}}(\Psi,\mathcal{S})=\frac{1}{M}\sum_{m=1}^{M}\|\Psi(u^{(m)})-\mathcal{S}(u^{(m)})\|_{L^2}^2. \end{align}\] We want to show that a global minimizer of the empirical risk over \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\), which we denote by \(\hat{\Psi}_{\mathcal{S}}\) has small population risk, \[\begin{align} \mathcal{L}(\Psi,\mathcal{S})=\mathbb{E}_{u\sim\mu}\|\Psi(u)-\mathcal{S}(u)\|_{L^2}^2, \end{align}\] with high probability. Our proof is adapted from [21]. First, in Lemma 20 we bound the population risk of the ERM by an approximation error term, which uses the approximation result in Theorem 2 and an estimation error term which we need to control. To control this estimation error term, we prove a concentration inequality for the population risk \(\mathcal{L}(\hat{\Psi}_{\mathcal{S}},\mathcal{S})\) in Proposition 7 in terms of the metric entropy of \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\) and \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\). Lemmas 25, 26, and 27 are used to prove this proposition. The major difference between our proof and that in [21] is that they require \(\|\mathcal{S}(u)\|_{L^2}\leq 1\), and our result holds for a generic upper bound \(V\). Throughout this section, we occasionally write \(\mathsf{S}\) and \(\Sigma\) for \(\mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) and \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\) to streamline the notation.

Using Lemmas 25, 26, and 27 from Appendix 8.3 along with a union-bound argument, we will show a concentration inequality for the risk uniformly over \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\times \mathsf{S}.\) To do so, we first recall the definitions of covering number and metric entropy.

Note that the covering number depends on both the set and metric. In what follows, we will consider covering numbers with respect to the supremum norm \[\begin{align} \mathsf{d}(\Psi,\Psi'):=\sup_{u\in \mathsf{U}_a}\|\Psi(u)-\Psi'(u)\|_{L^2(D)}. \end{align}\]

Consequently, we will drop the dependence on the metric in our notation for covering numbers in what follows. Throughout the analysis we allow external covers, meaning that the covering centers need not belong to the set being covered, but only to the ambient metric space. This convention is harmless here: any external \(\delta\)-cover can be converted into an internal \(2\delta\)-cover by choosing, from each nonempty covering ball, one point of the covered set.

The next result characterizes the concentration of the empirical risk over \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta})\times \mathsf{S}.\)

We now apply the previous lemma to prove a high probability bound on the population risk of an empirical risk minimizer, assuming we have enough samples relative to the metric entropy of the class of FNOs and operators. These metric entropies are controlled in the three subsequent results.

It remains to characterize the metric entropy of \(\Sigma\) and \(\mathsf{S}\), which we do in the following three results. The first is a restatement of a result from [21].

As an immediate corollary of Lemma 22, we can characterize the metric entropy of \(\Sigma^{\mathrm{poly}}(\varepsilon,\bar{\vartheta}).\)

The dependence of the metric entropy on the physical dimension \(d\) and the smoothness parameter \(a\) is natural. The factor \(\varepsilon^{-d/a}\) comes from the Fourier cutoff required to approximate functions in the \(H^a\) input class; equivalently, it reflects the number of retained Fourier modes, and hence the number of parameters needed to describe the Fourier multipliers. This number grows rapidly with the physical dimension \(d\). The smoothness parameter \(a\) dictates the rate at which the Fourier cutoff \(\mathscr{C}\) must grow in order to achieve a prescribed approximation accuracy. The remaining factor \(\varepsilon^{-4}\) comes from the depth of the time-stepping construction and the quadratic dependence on depth in the FNO metric entropy bound.

The next lemma transfers a cover of the approximating FNO class to a cover of the PDE solution-operator class. The condition \(\varepsilon\leq\delta\) ensures that the approximation error is no larger than the covering scale.

The proof of Theorem 3 now follows in a straightforward manner from the above lemmas.

8.3 Auxiliary Lemmas↩︎

In the proof of Theorem 2 we use the following classical Gronwall lemma.

In the proof of Theorem 3 in Appendix 8.2, we use the following lemmas. The first of these lemmas relies on the boundedness of any \(\mathcal{S}\in \mathsf{S}^{\mathrm{poly}}(\bar{\vartheta})\) and Bernstein’s inequality to obtain a concentration inequality for the excess risk.

The following lemma gives a Lipschitz bound for the risk functionals.

9 Supplementary Materials Section 4↩︎

9.1 Proof of Approximation Bound↩︎

9.2 Proof of Learning Bound↩︎

9.3 Auxiliary Lemmas↩︎

In the proof of the approximation bound, we used the following discrete Gronwall inequality [59].

9.4 Logarithmic Cahn–Hilliard Results↩︎

References↩︎

---

1. In [35], they show the stability and error bound for the spectral method applied to initial condition \(\mathcal{P}_N u_0,\) not \(\mathcal{I}_Nu_0,\) however the proofs carry through exactly as written with this change, provided that \(u_0\in \mathsf{U}_a\) for \(a>\frac{d}{2}\). This is also the case for the methods for the Allen–Cahn and Cahn–Hilliard equations.↩︎