Multi-stage neural operator learning with application for convolutions

Multi-stage neural operator learning with application for convolutions Zhiping Maoa , Zhenye Wenb , Yong Zhangc,d , Xiaofei Zhaob,e a School of Mathematical Sciences, Eastern Institute of Technology, Ningbo, 315200 Ningbo, China b School of Mathematics and Statistics, Wuhan University, 430072 Wuhan, China c Center for Applied Mathematics and KL-AAGDM, Tianjin University, 300072 Tianjin, China d State Key Laboratory of Synthetic Biology, Tianjin University, Tianjin, 300072, China

arXiv:2608.18851v1 [cs.LG] 19 Aug 2026

e Computational Sciences Hubei Key Laboratory, Wuhan University, 430072 Wuhan, China

Abstract Convolution integrals widely exist in applications, and to enable fast and accurate computations, this paper introduces two general multi-stage neural operator learning frameworks. The first, Deep Collocation Neural Operator (DCNO), is a supervised approach that iteratively refines the operator approximation by learning residuals from input-output data pairs. The second, Deep Galerkin Neural Operator (DGNO), is an unsupervised framework applicable when the target operator can be represented by a PDE, leveraging the weak form of the PDE residual for training. Both methods progressively construct basis operators through multiple training stages to enrich the approximation space, leading to significantly improved accuracy over standard one-shot operator learning. We provide theoretical analysis for their approximation capabilities and implement them for learning convolutions. Extensive numerical experiments demonstrate that both DCNO and DGNO achieve high accuracy, approaching machine precision under single float for convolution problems, and offer substantial efficiency gains for numerous queries or parametric variations compared to traditional solvers. We also extend these frameworks to handle multi-input operator learning scenarios involving variations in both the density and kernel of a convolution. Keywords: Operator learning, Convolution, Residual-driven approximation, Multiple stages, Deep operator network, High accuracy

1. Introduction Convolution computations arise naturally in a wide range of scientific and engineering applications, such as potential evaluation in electrostatics, gravitational particle simulations and signal processing [19, 48, 38, 17, 23], playing a crucial role in efficiently capturing long-range interactions. This paper focuses on the accurate and efficient evaluation of convolution. Let G denote the target convolution operator mapping the density ρ(x) to the potential ϕ(x), i.e., Z ϕ(x) = G(ρ)(x) = K(x, x̃)ρ(x̃) dx̃, x ∈ R ⊂ Rd (1.1) R

where d ∈ N+ is the spatial dimension, R denotes a region, and K(x, x̃) denotes a kernel. Despite the linearity of the mapping, the nonlocal nature of convolution (1.1) renders its evaluation computationally demanding under direct discretizations. Fast solvers such as fast multipole methods, hierarchical matrix techniques, Gaussian-Summation method, and kernel truncation method have therefore been extensively developed to enable efficient computation [19, 48, 20, 15, 42] of convolution for an input ρ. They are able to produce predicted ϕ on N ∈ N+ prescribed grids at cost O(N log N ). Although effective, Email addresses: zmao@eitech.edu.cn (Zhiping Mao), wenzhy@whu.edu.cn (Zhenye Wen), Zhang_Yong@tju.edu.cn (Yong Zhang), matzhxf@whu.edu.cn (Xiaofei Zhao)

repeated calls of the traditional solver for new density inputs in scenarios like parametric variations or manyquery settings may still incur substantial costs. This forms a natural circumstance for operator learning to apply, where neural networks are set to learn the mapping from the input function space (e.g., density and/or kernel) to the output solution space (e.g., potential). The trained operator network will be able to offer a fast prediction of the solution function for a wide range of input functions at arbitrary points of the region [34, 9, 43, 28]. In contrast, traditional solvers produce the predicted solution only on the prescribed grid points, and values elsewhere will need to be obtained posteriorly via proper interpolations to maintain accuracy, e.g., [18]. In the literature, a wide range of neural operator learning methods have been developed, with representative approaches including Fourier neural operators [28, 27], deep operator networks (DeepONet) [34, 24, 35, 44], graph kernel networks [5], Orthogonal Polynomial Neural Operator [32], and random feature-based models [37]. Successful applications have already been made in a broad range of scientific computing problems [10, 11, 14, 26, 29, 36], and for the particular studies on the convolution problem (1.1), we refer to [34, 8]. Despite the remarkable progress, the existing operator learning works are based on direct one-shot training, whose accuracy is often quite limited. The resulting operator network for the convolution problem (1.1) would therefore be less competitive compared to traditional fast solvers. In fact, such a one-shot training accuracy limitation was noted earlier in deep learning of fixed-configuration PDE problems [2, 46, 45], and to address this issue, several residual-based approaches reaching high accuracy have been developed. They progressively seek new bases to improve solution accuracy by minimizing the current residual of the PDE. The Galerkin Neural Networks (GNN) [2] employs a Galerkin framework based on the weak formulation of PDEs to iteratively enrich the solution approximation, particularly for the problems involving singularities. Its extension, Extended Galerkin Neural Networks [3], further generalizes this framework to a broader class of boundary value problems. Under the strong formulation, the Multi-level Neural Networks [4] and the Multi-stage Neural Networks [45] progressively improve the approximation by sequentially introducing new neural network components to approximate the residual, whereas the Multigrade Neural Networks [47] achieves iterative refinement through progressive network deepening. More recently, the Deep Collocation Method [46] adopts a collocation-based least-squares framework to improve training efficiency while maintaining high accuracy. These high-accuracy PDE solvers offer valuable insights for developing high-accuracy operator learning frameworks, which will be one of our goals. To this end, in this work we shall first propose two high-accuracy operator learning frameworks that enrich the approximation space by progressively constructing neural operator bases as DeepONets through multiple residual-based training stages guided by the residuals. Specifically, a deep collocation neural operator (DCNO) method is designed for settings where input-output data pairs are available, falling within a supervised learning paradigm. On the other hand, when the task operator can be represented as a PDE, a deep Galerkin neural operator (DGNO) method is proposed, whose training uses the residual of the PDE in the weak form, enabling fully unsupervised operator learning. Theoretical analysis is done for each of them to support their approximation ability and convergence of the iterative frameworks. Afterwards, we shall apply the two proposed methods to the convolution operator learning task (1.1), and demonstrate through a series of numerical experiments and comparisons the accuracy and efficiency of the proposed algorithms. Extensions to the case with multiple inputs are implemented in the end using MIONet [24]. The main contributions of this work are summarized as follows. 1) Two general multi-stage operator learning frameworks are proposed, including a supervised one based on data pairs and an unsupervised one based on PDE. They are implemented using DeepONet and analyzed for accuracy and convergence. 2) The proposed algorithms are applied to convolution operator learning with extensions made using MIONet. Experimental results demonstrate efficient predictions with accuracy that approaches machine precision under single-precision floating-point arithmetic, making the proposed methods a competitive addition to the state-of-the-art solvers for convolution. The rest of this paper is organized as follows. Section 2 provides the necessary preliminaries on residualbased neural network methods and operator networks. Sections 3 and 4 present and analyze the DCNO and 2

DGNO methods, respectively. Applications to convolution operator learning problems are given in Section 5. Some conclusions are drawn in Section 6. 2. Preliminary In this section, we briefly review some preliminary knowledge that is fundamental for proposing our method. The first part presents a class of residual-based neural network methods for high-accuracy PDE solving, and the second part presents the neural network for operator learning. 2.1. Residual-based neural network methods for high-accuracy PDE solving Although deep neural networks are capable of representing intricate solution structures of PDEs [40], direct one-shot training often faces limitations in accuracy. Residual-based neural network methods in the literature have addressed this challenge by iteratively constructing adaptive approximation spaces using multiple neural networks, thereby progressively improving approximations for PDEs. For illustration here, we briefly present a representative approach, namely the Galerkin neural network (GNN) [2]. Consider the following PDE problem in a bounded domain Ω ⊂ Rd with boundary ∂Ω: ( Lu(x) = f (x), x ∈ Ω, (2.1) boundary conditions, x ∈ ∂Ω, where L is a linear differential operator and f (x) is a source term. The bounded setup here is natural for computations. Assume that (2.1) can be reformulated into a variational (weak) form: find u∈X

such that

a(u, v) = F (f, v),

∀v ∈ X.

(2.2)

Here, X denotes the solution space, F (f, ·) : X → R is a bounded linear functional on X induced by the source term f , and a(·, ·) : X × X → R is a bounded, symmetric, bilinear form defined on X satisfying a(v, v) ≥ 0, and a(v, v) = 0 ⇐⇒ v = 0. Under this p assumption, a(·, ·) induces an inner product on X, with the corresponding norm defined by ∥v∥X = a(v, v). For example, in (2.1), when L = −∆ with homogeneous Dirichlet boundary R conditions, 1 one takes X = H (Ω), and the operator a and F associated with (2.1) read a(u, v) = ∇u · ∇v dx, 0 Ω R F (f, v) = Ω f v dx. The Galerkin subspace within GNN is constructed as follows. Starting from an initial basis function ψ0 = 0, at stage s ≥ 1, let Ψs−1 = span{ψ0 , . . . , ψs−1 } ⊂ X denote the current approximation space spanned by the basis functions {ψj }s−1 j=0 , and let us−1 ∈ Ψs−1 be the corresponding Galerkin approximation. In principle, the most effective enrichment direction is given by the error es−1 := u − us−1 , since adding this direction to the approximation subspace would recover the exact solution in a single step. Thus, GNN seeks a neural-network basis function ψs that approximates the normalized error φs−1 = (u − us−1 )/∥u − us−1 ∥X . Since the exact solution u is not accessible, the residual functional Rs−1 (f, v) := F (f, v) − a(us−1 , v) can be used as a substitute due to the fact that the normalized error φs−1 mathematically maximizes Rs−1 (f, v) over the unit ball in X [2]: φs−1 ∈ arg max Rs−1 (f, v),

Rs−1 (f, φs−1 ) = ∥u − us−1 ∥X .

(2.3)

∥v∥X =1

This characterization suggests the construction of the new basis function: ψs = arg max v̸=0

Rs−1 (f, v) , ∥v∥X

(2.4)

and with the enriched finite-dimensional subspace Ψs = span{ψ0 , . . . , ψs } ⊂ X, the Galerkin approximation us ∈ Ψs is determined by solving a (s + 1) × (s + 1) linear system us ∈ Ψs :

a(us , ψj ) = F (vj ) 3

0 ≤ j ≤ s,

(2.5)

for the linear expansion coefficients of us on the basis of Ψs . In addition to GNN, residual-based neural network methods based on the strong form of PDE (2.1) have also been proposed, such as Multi-level Neural Networks [4] and the Deep Collocation Method [46]. We shall develop some high-accuracy neural operator learning methods by leveraging the residual-based framework. 2.2. Deep operator network In this work, we adopt DeepONet as the neural architecture for operator approximation, designed to learn mappings between infinite-dimensional function spaces [34]. Consider an operator G : W → X, mapping an input function ρ ∈ W to an output function ϕ ∈ X, where W and X both are (for simplicity) function spaces defined on Ω ⊂ Rd with d ∈ N+ . For any x ∈ Ω, DeepONet approximates G(ρ) via a trunk-branch decomposition: p X G(ρ)(x) ≈ Gθ (ρ)(x) = bk (ρ(s1 ), . . . , ρ(sm ); θb ) tk (x; θt ), (2.6) k=1

where θ = {θb , θt } denotes the trainable parameters, with θb and θt representing the parameters of the branch and trunk networks, respectively. Here, {sl }m l=1 ⊂ Ω are sensor locations, and latent dimension p ∈ N+ is the output width of branch and trunk networks. The branch network is responsible for extracting a p-dimensional latent representation b = (b1 , b2 , . . . , bp ) ∈ Rp of the input function ρ and is typically implemented as a multilayer perceptron (MLP) [41]. In the standard setting, the branch input consists of the evaluations of ρ at fixed sensors {sl }m l=1 . More generally, the branch input can also be taken as the expansion coefficients of ρ with respect to a prescribed basis of W [12, 34, 35]. The trunk network captures the dependence of the output function on the coordinate x ∈ Ω. It is implemented as another MLP that maps x to a latent vector t = (t1 , t2 , . . . , tp ) ∈ Rp , which can be interpreted as a set of learned basis functions over the output domain. The operator approximation is then constructed through the inner product of b and t, resulting in a continuous representation of G. Since its introduction, DeepONet has been successfully applied to a wide range of problems [10, 11, 26, 29]. However, its standard supervised training paradigm, which relies on labeled data and one-shot optimization, often suffers from limited approximation accuracy. The aim of this work, on the one hand, is to develop a general multi-stage framework for high-accuracy operator learning by progressively constructing basis operators guided by residuals and improving the accuracy of operator approximation. On the other hand, the generation of high-fidelity training data for labels typically requires repeated calls to traditional numerical solvers, leading to substantial computational cost. Thus, in addition to the multi-stage framework, we also work for an unsupervised operator learning manner that eliminates the need for labeled data. 3. Deep collocation neural operator (DCNO) method In this section, we introduce the deep collocation neural operator (DCNO) method, a supervised framework for high-accuracy operator learning. Instead of training a single, monolithic neural operator to approximate a target operator G(ρ) = ϕ directly, DCNO constructs a sequence of basis operators and then approximates the target operator as a linear combination of trained basis operator networks. The following first presents the method and then provides some analysis to support its approximation validity. 3.1. Neural network basis operators and collocation scheme Let us first outline the basic framework. We construct the basis operators through a greedy iterative strategy similarly as DCM and GNN [2, 46] (see also [4, 45, 47]), where each basis operator is implemented as DeepONet (2.6) to approximate the solution residual generated at the current iteration. Beginning with an initial operator G0 := Gθ0 , which is chosen as the zero operator in this work, the initial residual is defined as e0 = G − Gθ0 = G. The first basis operator Gθ1 is then trained to approximate the initial residual mapping, i.e., Gθ1 ≈ e0 . The resulting approximation for G is then updated to G1 = Gθ0 + Gθ1 . Assume that we Ps−1 are now at some stage s ≥ 2 with the approximation from the previous stage Gs−1 = j=0 Gθj , then the corresponding residual operator is defined as es−1 = G − Gs−1 . By training a new basis operator Gθs to 4

Ps approximate the residual mapping ρ 7→ es−1 , the update Gs = j=0 Gθj thereby incrementally improves the overall approximation for G. Now, we detail the implementation. We consider here for DCNO in the particular case that solution data (b) b {(ρ(b) , ϕ(b) )}N are sampled inputs with Nb ∈ N+ and ϕ(b) = G(ρ(b) ) are possibly b=1 are available, where ρ obtained from numerical solvers. Then, each stage can be trained in a supervised manner, as detailed below. Let Ω denote the computational (bounded) domain in practice. The trunk network in (2.6) takes the spatial x coordinate x ∈ Ω as input. A set of collocation points is then generated through discretization: {xk }N k=1 ⊂ Ω, where Nx ∈ N+ denotes the number of collocation points. Accordingly, for each input sample ρ(b) , b = 1, . . . , Nb , with corresponding output ϕ(b) , the training data are constructed as {(ρ(b) (xk ), ϕ(b) (xk )) : b k = 1, . . . , Nx }N b=1 . Let G0 = 0. For s ≥ 1, the residual corresponding to the b-th training sample at stage s − 1 is given by es−1 (ρ(b) )(x) := ϕ(b) (x) − Gs−1 (ρ(b) )(x), b = 1, . . . , Nb . (3.1) We measure it under a discrete norm: v u Nb X Nx u X 1 |es−1 (ρ(b) )(xk )|2 =: βs−1 , ∥es−1 ∥ = t Nb × N x

s ≥ 1.

(3.2)

b=1 k=1

As long as βs−1 > 0, DCNO could in principle proceed to the next stage and the DeepONet Gθs is determined by minimizing the empirical loss N

Losss (θ) :=

N

2

b X x X 1 1 es−1 (ρ(b) )(xk ) − Gθ (ρ(b) )(xk ) , Nb × N x βs−1

b=1 k=1

and

θs ∈ arg min {Losss (θ)} . (3.3) θ

Here, the residual is scaled to order one for practical training efficiency [4], and the approximate operator for G at stage s reads: s X Gs (ρ) := βj−1 Gθj (ρ), s ≥ 1. (3.4) j=1

In our numerical experiments later, the unconstrained optimization problem in (3.3) for every s ≥ 1 will be solved by Adam optimizer [25] with the Xavier method [16] for parameter initialization. It is common for the optimizer to produce a practical minimizer for (3.3) with Losss (θs ) > 0, which is in fact fine here since DCNO will pass the error as some residual to the next stage for refinement. This process can be terminated once βs = ∥es ∥ ≤ ϵ for some s, where ϵ > 0 is a prescribed tolerance. This means that the target operator has been approximated to the desired accuracy on the training set. If the training data are sampled sufficiently large and effective, the final DCNO approximation G(ρ) ≈ Gs (ρ) can therefore become valid for any input function ρ (inside or outside the training set) from a desired functional space. Note that by construction, es = G − Gs = es−1 − βs−1 Gθs , and so the stopping criterion ∥es ∥ ≤ ϵ can be replaced equivalently as p βs−1 Losss (θs ) ≤ ϵ. The whole process of the DCNO method is summarized in Algorithm 1. Note that DCNO itself does not require the target operator G to be linear, despite that we will apply it later to convolutions. 3.2. Numerical analysis Here we present some numerical analysis for the proposed DCNO method. Let G̃(ρ) and G(ρ) denote the exact operator and the reference operator, respectively, where the latter serves as the training label in our algorithm P and may deviate from the exact one. They both map W to X ⊂ L2 (Ω). With the numerical s solution Gs = j=1 βj−1 Gθj at the s-th stage, the stage-wise error is defined by es (ρ) := G(ρ) − Gs (ρ), 5

(3.5)

Algorithm 1 DCNO method Parameters: tolerance ϵ; maximum step Ne for optimizer Nx b Input: Training dataset {(ρ(b) , ϕ(b) )}N b=1 ; collocation points {xk }k=1 Output: Basis operator set {Gθj }sj=1 ; Solution operator Gs (ρ) ≈ G(ρ) 1: Initialize: s = 0; G0 = 0; 2: Part I: Construction of neural operator basis 3: do 4: Set s ← s + 1 5: Compute es−1 (ρ(b) ) (3.1) for every b and compute (3.2) for βs−1 6: Initialize DeepONet parameters θs 7: for epoch = 1 to epoch = Ne do 8: compute the loss function by (3.3) 9: update θs via the optimizer 10: end for Ps 11: Update approximation Gs = j=1 βj−1 Gθj p Losss (θs ) > ϵ 12: while βs−1 13: Part II: Inference for new input ρ 14: Given ρ, evaluate {Gθj (ρ)}sj=1 Ps 15: Form the prediction Gs (ρ) = j=1 βj−1 Gθj (ρ)

and by taking the data error ẽ(ρ) := G̃(ρ) − G(ρ) into account, we further define ẽs (ρ) := G̃(ρ) − Gs (ρ).

(3.6)

We consider the DeepONets for stage-wise construction of DCNO to be from D := {Gθ (ρ)(x) = b(ρ) · t(x) | b ∈ Nbranch , t ∈ Ntrunk } ,

(3.7)

where Nbranch and Ntrunk are classes of continuous neural network functions (e.g., the MLP class), and each picked b, t are two vectors of the same size so that their inner product makes sense. To quantify the approximation performance over the input space W, we consider a finite measure µ supported on W and introduce the Bochner space L2 (W; L2 (Ω)) equipped with the norm [1] Z ∥v∥2L2 (W;L2 (Ω)) := ∥v(ρ)∥2L2 (Ω) dµ(ρ), v : W → L2 (Ω). (3.8) W

Note that (3.2) can be interpreted as a realization/discretization of (3.8), so we will still denote βs = ∥es ∥L2 (W;L2 (Ω)) for simplicity in the following. For the error estimate of the DCNO construction process, we make the following basic assumption. Assumption 3.1. Suppose that G : W → X is a Rcontinuous operator, where W is a compact subset of C(Ω) and X = C(Ω) with 0 < |Ω| < ∞ and 0 < C := W dµ(ρ) < ∞. Before the construction process of DCNO in Algorithm 1 stops at some S ∈ N+ , the residual operator is non-zero, i.e., βs = ∥es ∥L2 (W;L2 (Ω)) > 0,

s = 0, 1, . . . S − 1.

We assume that the training and sampling within DCNO are sufficient at each stage so that the resulting DeepONet basis p Gθs for s = 1, . . . S do not pollute the universal approximation ability of D up to a prescribed threshold τ / C|Ω|. Under Assumption 3.1, we can directly have the following error estimate result. 6

Theorem 3.2. Fix one 0 < τ < 1 and let Assumption 3.1 hold. Then, the error of DCNO is strictly decreasing during construction of each stage, i.e., ∥es ∥L2 (W;L2 (Ω)) < ∥es−1 ∥L2 (W;L2 (Ω)) , and the error at the s-stage satisfies the bound ∥es ∥L2 (W;L2 (Ω)) ≤ τ s ∥e0 ∥L2 (W;L2 (Ω)) ,

s = 1, . . . , S.

(3.9)

Taking into account the data error ẽ, we further have ∥ẽs ∥L2 (W;L2 (Ω)) ≤ ∥ẽ∥L2 (W;L2 (Ω)) + τ s ∥e0 ∥L2 (W;L2 (Ω)) ,

s = 1, . . . , S.

(3.10)

Proof. Since each Gθs is taken from D, we have Gs as defined in (3.4) is a continuous operator and so does es (ρ) given by (3.5) for each s. By the universal approximation theorem for DeepONet [34], the basis operator network generated within DCNO at stage s, i.e., Gθs ∈ D for 1 ≤ s ≤ S, can approximate the last scaled residue es−1 /βs−1 to a range: τ 1 , es−1 (ρ)(x) − Gθs (ρ)(x) ≤ p βs−1 C|Ω|

∀ρ ∈ W, x ∈ Ω.

(3.11)

By taking the square on both sides of (3.11) and then integrating with respect to x ∈ Ω, we have 1 βs−1

2

es−1 (ρ) − Gθs (ρ)

Z =

L2 (Ω)

Ω

1 βs−1

2

es−1 (ρ)(x) − Gθs (ρ)(x)

dx ≤

τ2 , C

∀ρ ∈ W.

Further integrating the above with respect to ρ ∈ W under the measure µ and then taking the square root, 1 we are led to ∥ βs−1 es−1 − Gθs ∥L2 (W;L2 (Ω)) ≤ τ. By the construction of DCNO and the definition (3.5), we have es = es−1 − βs−1 Gθs , and so ∥es ∥L2 (W;L2 (Ω)) ≤ τ ∥es−1 ∥L2 (W;L2 (Ω)) ,

s = 1, . . . , S.

(3.12)

As τ < 1, the error strictly decreases. Applying the estimate (3.12) recursively, we obtain (3.9). Then, by subtracting (3.6) from (3.5), we see that ẽs (ρ) = ẽ(ρ) + es (ρ). Applying the triangle inequality to it yields ∥ẽs ∥L2 (W;L2 (Ω)) ≤ ∥es ∥L2 (W;L2 (Ω)) + ∥ẽ∥L2 (W;L2 (Ω)) . Using (3.9) further gives (3.10). 4. Deep Galerkin neural operator (DGNO) method Alternative to DCNO, in this section, we develop an unsupervised high-accuracy operator learning framework, termed the deep Galerkin neural operator (DGNO) method, for approximating the solution operator of a PDE. For illustration, the target PDE here is set as (2.1) with G : f 7→ u, the mapping of the input source term to the corresponding solution. Similarly as DCNO, DGNO also builds a finite sum of basis operator networks for approximation, but the training is fully based on equation residue without requiring labeled data. 4.1. Constructing neural network basis operator in an unsupervised manner As we are working with a PDE, to reduce the dependence of regularity, we choose to build DGNO for (2.1) upon its variational formulation (2.2), via a framework analogous to GNN. The progressive construction procedure now is to find DeepONets as basis neural operators based on residuals to enrich the approximation space. The iterative procedure starts with a zero initial operator G0 = Gθ0 = 0, then the initial normalized error under an input f is φ0 (f ) = (G(f ) − G0 (f ))/∥G(f ) − G0 (f )∥X = G(f )/∥G(f )∥X . The first basis operator 7

Gθ1 is trained to approximate the initial error mapping, i.e., Gθ1 ≈ φ0 , and the resulting approximation for G is then updated to G1 := β1 Gθ1 , where the coefficient β1 will be determined via a Galerkin projection detailed later.PAssume that we are now at some stage s ≥ 2 with the approximation from the previous s−1 stage Gs−1 = j=1 βj Gθj , the new basis operator Gθs is trained to approximate the current optimal update direction, namely the normalized error operator read under an input f as φs−1 (f ) :=

G(f ) − Gs−1 (f ) . ∥G(f ) − Gs−1 (f )∥X

(4.1)

For fixed f , by (2.3), φs−1 (f ) can be achieved by maximizing the normalized PDE-residual functional: φs−1 (f ) ∈ arg max Rs−1 (f, v),

Rs−1 (f, v) := F (f, v) − a(Gs−1 (f ), v).

where

∥v∥X =1

To extend this pointwise characterization to the operator level, we consider the expectation of the normalized residual functional with respect to a measure µ for input f , and further motivated from the GNN strategy (2.4), we are led to the following operator-level optimization problem:   Rs−1 (f, Gθ (f )) , s ≥ 1. (4.2) Gθs ∈ arg max Ef ∼µ ∥Gθ (f )∥X Gθ ∈D Here, D denotes the class of DeepONets defined in (3.7). The target function of the optimization problem (4.2) is purely based on the residue of PDE, without calling for labeled data, and thus provides an unsupervised way for training the basis operators. Now we give some details on implementation. At stage s ≥ 1, given the set of constructed operators Ps−1 (b) (b) Nb {Gθj }s−1 }b=1 , and the approximate solution Gs−1 (f (b) ) = j=0 βj Gθj (f (b) ). j=0 , the training set of input {f We look for the basis operator Gθs ≈ φs−1 based on (4.2). Here, the expectation is discretized as the empirical average over the training samples in practice, and then we are motivated to train θs of Gθs by minimizing the loss N

Losss (θ) = −

b Rs−1 (f (b) , Gθ (f (b) )) 1 X , Nb ∥Gθ (f (b) )∥X

b=1

θs ∈ arg min {Losss (θ)} .

(4.3)

θ

To ensure high accuracy, both the residual functional Rs−1 (·, ·) and the energy norm ∥·∥X are computed using x a Gauss–Legendre quadrature rule on the quadrature points {xk }N k=1 . With the obtained basis operator set s {Gθj }j=1 and an input f , we seek a Galerkin approximation Gs (f ) for the solution of (2.2) in span{Gθj (f )}sj=1 , Ps i.e., Gs (f )(x) = j=1 βj (f ) Gθj (f )(x), such that a(Gs (f ), Gθj (f )) = F (f, Gθj (f )),

1 ≤ j ≤ s.

The coefficient vector β(f ) = [β0 (f ), . . . , βs (f )]T is thus obtained by solving the linear system K(f ) β(f ) = F(f ),

(4.4)

where K = (Kk,j )s×s and F = (Fj )s×1 with entries
Kk,j (f ) = a Gθk (f ), Gθj (f ) ,

Fj (f ) = F f, Gθj (f ) .

The Galerkin projection scheme (4.4) ensures that for each input f , it provides the optimal approximation to the solution of the variational problem within the current approximation subspace given by the basis operators. The construction process of DGNO can be stopped when a prescribed maximum stage number S ≥ 1 is reached. The overall framework of DGNO is summarized in Algorithm 2. Remark 4.1. The current DGNO is designed based on the bilinear form (2.2). Nevertheless, by incorporating linearization techniques such as fixed-point iteration, the proposed framework could be extended to general nonlinear case, which will be investigated in our future work. 8

Algorithm 2 DGNO method Parameters: maximum stage number S; maximum step Ne for optimizer Nx b Input: Training dataset {f (b) }N b=1 ; Gaussian quadrature points {xk }k=1 Output: Basis operator set {Gθj }Sj=1 ; Solution GS ≈ G(f ) 1: Initialize: G0 ← 0; Ψs ← ∅ 2: Part I: Greedy construction of neural operator basis 3: for s = 1 to S do 4: Initialize DeepONet parameters θs 5: for epoch = 1 to epoch = Ne do 6: compute the loss function by (4.3) 7: update θs via the optimizer 8: end for 9: Update basis set: Ψs+1 ← Ψs ∪ {Gθs } 10: Compute coefficients β(f (b) ) for every Ps b via (4.4) 11: Update approximation Gs (f (b) ) = j=1 βj (f (b) ) Gθj (f (b) ) for b = 1, · · · , Nb 12: end for 13: Part II: Inference for new input f 14: Given f , evaluate {Gθj (f )}S j=1 15: Compute coefficient β(f ) via solving linear system (4.4) PS 16: Form the prediction GS (f )(x) = j=1 βj (f )Gθj (f )(x)

4.2. Numerical analysis Here we perform some numerical analysis for the proposed DGNO method. Our analysis considers a continuous projection neglecting the quadrature errors. To quantify the overall approximation performance over the input space W, we consider a probability measure µ supported on W and introduce the Bochner space L1 (W; X) equipped with the norm [1] Z ∥v∥L1 (W;X) := ∥v(f )∥X dµ(f ), v ∈ L1 (W; X). (4.5) W

In the context of solving PDE, it is natural to consider the solution space X as some Sobolev space and we make the following basic assumption. Assumption 4.2. For some 0 ≤ q1 , q2 < ∞, let the target operator G : W → X be a continuous operator, where W is a compact subset of H q1 (Ω) and X is subspace of H q2 (Ω) with the energy norm ||·||X continuously controlled by ∥ · ∥H q2 (Ω) and Ω a bounded open domain. We first establish a universal approximation result for DeepONet in the Sobolev space. To do this, we begin with the following lemma. Lemma 4.3 ([22]). Suppose 0 ≤ q < ∞ and Ω ⊂ Rd is a bounded open domain. If σ ∈ C q (R) is non-constant and bounded, then the space   ∞  N  [ X Vσ = fN : Rd → R fN (x) = αj σ(wj · x + zj ), ∀αj , zj ∈ R, wj ∈ Rd .   N =1

j=1

is dense in the Sobolev space H q (Ω). Based on the neural-network approximation for functions in Sobolev spaces in Lemma 4.3, we establish the following universal approximation in Sobolev space for operator networks. 9

Lemma 4.4 (Universal approximation for stacked DeepONet in Sobolev norm). Let Assumption 4.2 hold. If σ ∈ C q2 (R) is non-constant and bounded, then for any ϵ > 0, there exists a stacked DeepONet Gθ (f )(x) := P p l=1 bl,θb (cK (f )) tl,θt (x), such that: sup ∥G(f ) − Gθ (f )∥H q2 (Ω) < ϵ. f ∈W

PNl l PMl l Here, x ∈ Ω, bl,θb (cK (f )) := k=1 αk σ(wkl · cK (f ) + zkl ), tl,θt (x) := j=1 αj σ(wjl · x + zjl ). The coefficient ∞ K vector cK (f ) = (⟨f, ζ1 ⟩H q1 (Ω) , . . . , ⟨f, ζK ⟩H q1 (Ω) ) ∈ R , where {ζk }k=1 is an orthonormal basis of H q1 (Ω). PK Proof. Define the projection operator PK : H q1 (Ω) → H q1 (Ω) by PK f := k=1 ⟨f, ζk ⟩H q1 ζk . Since W ⊂ H q1 (Ω) is compact and {PK }K∈N is a sequence of uniformly bounded projections, we have supf ∈W ∥PK f − f ∥H q1 (Ω) → 0 as K → ∞. Moreover, since G is continuous on the compact set W, it is uniformly continuous. Therefore, for any ϵ > 0, there exists δ > 0 such that ∥f −PK f ∥H q1 (Ω) < δ =⇒ ∥G(f )−G(PK f )∥H q2 (Ω) < ϵ/4. Choosing K sufficiently large so that supf ∈W ∥PK f − f ∥H q1 (Ω) < δ, we obtain sup ∥G(f ) − G(PK f )∥H q2 (Ω) < f ∈W

ϵ . 4

q2 q1 Let {ηl }∞ l=1 be an unit orthonormal basis of H (Ω). Since PK (W) is compact in H (Ω) and G is q2 continuous, the set G(PK (W)) is compact in H (Ω). Similarly, there exists a positive integer p such that

sup G(PK f ) − f ∈W

p X ⟨G(PK f ), ηl ⟩H q2 ηl l=1

H q2 (Ω)

<

ϵ . 4

For each l = 1, . . . , p, define Bl (cK (f
)) := ⟨G(PK f ), ηl ⟩H q2 , f ∈ W. Since PK f is uniquely determined by cK (f ) = ⟨f, ζ1 ⟩H q1 , . . . , ⟨f, ζK ⟩H q1 , the above definition is well posed. Moreover, Bl is continuous on cK (W) as it is the composition of continuous mappings. Since W is compact in H q1 (Ω), cK (W) is a compact subset of RK . By the universal approximation theorem for continuous functions [13], for each l = 1, . . . , p, PNl l there exists a neural network bl,θb (cK (f )) = k=1 αk σ(wkl · cK (f ) + zkl ), such that ϵ sup |Bl (cK (f )) − bl,θb (cK (f ))| < √ . 4 p

f ∈W

Since p is finite, these approximations can be chosen simultaneously for all l = 1, . . . , p. Then, sup f ∈W

p X

Bl (cK (f )) − bl,θb (cK (f )) ηl

l=1

H q2 (Ω)

= sup f ∈W

≤

p X l=1

p X

!1/2 2

|Bl (cK (f )) − bl,θb (cK (f ))|

l=1

!1/2 sup |Bl (cK (f )) − bl,θb (cK (f ))|2

<

f ∈W

ϵ . 4

PMl l Now, define the network tl,θt (x) := j=1 αj σ(wjl · x + zjl ), 1 ≤ l ≤ p. Since p is finite, for each 1 ≤ l ≤ p, by the universal approximation from Lemma 4.3, there exists Ml such that ∥ηl −tl,θt ∥H q2 (Ω) < δ0 . Moreover, since cK (W) is compact and bl,θb is continuous, exists a constant C0 > 0 such that |bl,θb (cK (f ))| ≤ Pthere p C0 , ∀f ∈ W, 1 ≤ l ≤ p. Defining Gθ (f )(x) := l=1 bl,θb (cK (f ))tl,θt (x) and by the triangle inequality, we have sup f ∈W

p X l=1

bl,θb (cK (f ))ηl − Gθ (f )

H q2 (Ω)

≤ sup

p X

f ∈W

≤ C0

p X l=1

10

|bl,θb (cK (f ))| · ∥ηl − tl,θt ∥H q2 (Ω)

l=1

∥ηl − tl,θt ∥H q2 (Ω) ≤ C0 pδ0 .



By choosing δ0 = ϵ/(4C0 p), we obtain

Pp

l=1 bl,θb (cK (f ))ηl − Gθ (f )

H q2 (Ω)

< ϵ/4. Then we have

p  X Bl (cK (f ))ηl sup ∥G(f ) − Gθ (f )∥H q2 (Ω) ≤ sup ∥G(f ) − G(PK f )∥H q2 (Ω) + G(PK f ) − f ∈W

f ∈W

+

p X

l=1

[Bl (cK (f )) − bl,θb (cK (f ))]ηl

l=1

H q2 (Ω)

+

p X

H q2 (Ω)

bl,θb (cK (f ))ηl − Gθ (f )

l=1

 H q2 (Ω)

< ϵ.

Lemma 4.4 only considers DeepONet constructed from stacked shallow networks. The following extends the result by allowing for more general branch and trunk networks. Lemma 4.5 (Generalized universal approximation for DeepONet in Sobolev norms). Under Assumption 4.2, for any ϵ > 0, there exist positive integers K, p, continuous functionals c1 , c2 , . . . , cK : H q1 (Ω) → R, and continuous vector functions b : RK → Rp , t : Rd → Rp , such that t ∥H q2 (Ω) < ϵ. sup ∥G(f ) − b(c1 (f ), c2 (f ), . . . , cK (f )) · |{z} {z } |

f ∈W

branch

(4.6)

trunk

The functions b and t can be chosen as diverse classes of neural networks which satisfy the classical universal approximation theorem of functions, including (stacked/unstacked) fully connected neural networks and convolutional neural networks. PNl l Proof. Let b = (b1 , · · · , bp )T and t = (t1 , · · · , tp )T . We choose bl (cK (f )) = k=1 αk σ(wkl · cK (f ) + zkl ), and PMl l l l d tl (x) = j=1 αj σ(wj · x + zj ) for x ∈ R , where all the notations involved are defined in Lemma 4.4. Then, this result holds immediately following Lemma 4.4. Note that the domains of b and t can be restricted on compact sets. In fact, t is defined on Ω, and b can be defined on [−Z, Z]K , where Z > 0 is a uniform bound of the functions in W. Therefore, b and t can be approximated by neural networks according to the universal approximation of functions [13]. Now, we are ready to present a convergence result for the proposed DGNO method. Theorem 4.6. Under Assumption 4.2, fix one 0 < τ < 1 and assume that the training and sampling within DGNO are sufficient at each stage so that the resulting DeepONet basis Gθs for s = 1, 2, . . . S can meet (4.6) with ϵ = τ . Then, the following error estimates hold. 1) The error of DGNO is strictly decreasing during construction of each stage, i.e., ∥es ∥L1 (W;X) < ∥es−1 ∥L1 (W;X) , and the error at the s-stage satisfies the bound ∥es ∥L1 (W;X) < τ s ∥e0 ∥L1 (W;X) . 2) Define Js−1 (Gθs ) := bound

R W

(4.7)

Rs−1 (f, Gθs (f )) dµ(f ), then the error at stage s − 1 of DGNO satisfies the

1 1 Js−1 (Gθs ) < ∥es−1 ∥L1 (W;X) < Js−1 (Gθs ). 1+τ 1−τ

(4.8)

Proof. By the definition (4.1), φs−1 is a continuous operator from W ⊂ H q1 (Ω) to X ⊂ H q2 (Ω). Under the assumption and also by Lemma 4.5, the basis operator network generated within DGNO at stage s, i.e., Gθs ∈ D for 1 ≤ s ≤ S, can approximate the φs−1 to a range: ∥Gθs (f ) − φs−1 (f )∥X < τ, 11

∀f ∈ W.

(4.9)

For any fixed f ∈ W, note that Gs−1 + ∥es−1 (f )∥X Gθs (f ) ∈ span{Gθj (f )}sj=1 . Hence, by the property of Galerkin projection (see also [2]), we have ∥es (f )∥X = ∥G(f ) − Gs (f )∥X ≤∥G(f ) − Gs−1 (f ) − ∥es−1 (f )∥X Gθs (f )∥X =∥∥es−1 (f )∥X φs−1 (f ) − ∥es−1 (f )∥X Gθs (f )∥X =∥es−1 (f )∥X ∥φs−1 (f ) − Gθs (f )∥X

(4.10)

<τ ∥es−1 (f )∥X . where the final step follows (4.9). Integrating over W with respect f under the measure µ immediately yields ∥es ∥L1 (W;X) < τ ∥es−1 ∥L1 (W;X) . As τ < 1, the error strictly decreases. Applying this estimate recursively, we obtain (4.7). By the definition of Rs−1 , for every f ∈ W, Rs−1 (f, φs−1 (f )) = ∥es−1 (f )∥X . It is then direct to deduce that



Rs−1 f, Gθs (f ) − ∥es−1 (f )∥X = Rs−1 f, Gθs (f ) − Rs−1 f, φs−1 (f ) = Rs−1 f, Gθs (f ) − φs−1 (f )

= F f, Gθs (f ) − φs−1 (f ) − a Gθs−1 (f ), Gθs (f ) − φs−1 (f )


= a G(f ), Gθs (f ) − φs−1 (f ) − a Gθs−1 (f ), Gθs (f ) − φs−1 (f ) = a es−1 (f ), Gθs (f ) − φs−1 (f ) . By the Cauchy–Schwarz inequality and (4.10), we find

Rs−1 f, Gθs (f ) − ∥es−1 (f )∥X = a es−1 (f ), Gθs (f ) − φs−1 (f ) ≤∥es−1 (f )∥X ∥Gθs (f ) − φs−1 (f )∥X < τ ∥es−1 (f )∥X , Thus,
(1 − τ ) ∥es−1 (f )∥X < Rs−1 f, Gθs (f ) < (1 + τ ) ∥es−1 (f )∥X , and integrating over W immediately yields (4.8). 5. Application to convolution problems In this section, we apply the proposed DCNO and DGNO methods to solve the convolution problem (1.1). We first detail the implementation of the methods for the convolution operator, followed by a series of numerical experiments to validate and compare accuracy and efficiency. Extensions to the multi-input case are made in the end. All experiments1 are conducted using a single CPU with 32 GB RAM for CPU-based methods and a single NVIDIA A800-PCIE-80GB GPU for GPU-based methods. 5.1. Implementation of DCNO & DGNO on convolution Focusing on the convolution problem (1.1): ϕ(x) = G(ρ)(x) = (K ∗ ρ)(x), the target operator G maps the input density ρ to the solution potential ϕ. The generality of G can separate the use of DCNO from DGNO. When G does not have a precise (local) PDE formulation, then DGNO might not be applicable. In b such a case, we would employ DCNO for supervised training using data pairs {(ρ(b) , ϕ(b) )}N b=1 that could be generated by traditional numerical convolution solvers. By sampling the input density function ρ(b) , the solver produces accurate enough ϕ(b) = G(ρ(b) ). After the training dataset is constructed, the learning of the convolution operator via DCNO can be performed as in Algorithm 1. When G can be expressed through a proper PDE, i.e., (2.1), DGNO becomes available and its unsupervised training can certainly be rewarded with more efficiency. The task of learning the operator f 7→ u in (2.1) then follows the DGNO Rprocedure in Algorithm 2. If the original convolution problem (1.1) is imposed in the whole space: ϕ(x) = Rd K(x, x̃)ρ(x̃) dx̃, the resulting equivalent whole PDE problem, for computational reasons, usually needs to be truncated onto a bounded domain Ω with appropriate/approximate 1 Code is available at https://github.com/wzhy0777/DCNO-DGNO

12

boundary conditions. Careless truncations may cause loss of accuracy or require a very large domain size. To enlarge computational efficiency in such a case, we here combine with a moment-matching strategy [31] briefly presented below, which can accurately approximate the whole-space convolution as a local PDE with zero boundary conditions in a bounded domain. The moment-matching strategy begins with a decomposition [7, 31] on the density: ρ = ρmm + f with ρmm constructed for moment-matching and f the residual. Here, ρmm is assumed to be ρmm (x) = 2 Pm α 2 −d/2 exp(− |x| |α|=0 γα ∂ G(x), with G(x) = (2πγ ) 2γ 2 ), where α the multi-index and the coefficients {γα } are determined by imposing the moment-matching condition up to order m: Z Z ρmm (x)xα dx = ρ(x)xα dx, for |α| = 0, 1, . . . , m. (5.1) Rd

Rd

R R Pm α β β Conditions in (5.1) read equivalently as a linear system: |α|=0 γα Rd ∂ G(x)x dx = Rd ρ(x)x dx, for |β| = 0, 1, . . . , m, which admits an analytical solution [31]. The advantage of this construction is that ϕmm = K ∗ ρmm can be analytically evaluated and for the residual u = ϕ − ϕmm = K ∗ f , the condition (5.1) ensures a rapid decay of f at infinity. This in turn can lead to a rapid far-field decay also on u for K satisfying suitable far-field asymptotic behavior [31], which is the case for Green’s functions of common elliptic operators like the Laplacian. Then, the unbounded-domain problem can be truncated to a Ω not too large and u can be accurately approximated as the solution of ( Lu(x) = f (x), x ∈ Ω, (5.2) u(x) = 0, x ∈ ∂Ω. b Coming back to the convolution operator learning task, with the sampled training set {ρ(b) }N b=1 , now under (b) b this strategy, we first compute the corresponding {f (b) = ρ(b) − ρmm }N b=1 as inputs associated with (5.2) and then DGNO in Algorithm 2 can be applied for GS (f ) ≈ u. Thus, for any given ρ, the final prediction for the convolution outcome ϕ under our method is ϕmm + GS (f ). By the structure of DeepONet (2.6), the trained GS from DCNO or DGNO can give the predicted solution at N ∈ N+ arbitrary spatial points with computational complexity O(N ). In comparison, most existing traditional solvers such as the hierarchical matrix technique [21], the Gaussian-Summation method [15] and the kernel truncation method [30, 42] cost O(N log N ) for the prescribed N grid points. For a new spatial point outside the prescribed grid, traditional solvers will need to do extra delicate interpolations, e.g., the NUFFT [18], if one wants to maintain high accuracy. For the traditional method SOE [50], which is specially designed to achieve O(N ) cost for one-dimensional convolution, we will provide a numerical comparison below.

5.2. Numerical experiments We now carry out some numerical experiments to illustrate the performance of DCNO and DGNO. The involved DeepONet (2.6) uses identical depth L and width W for both branch and trunk networks, with the same number of neurons in each hidden layer for simplicity. Example 5.1 (Poisson potential). We begin with a 2D Poisson potential in (1.1), i.e., d = 2, x = (x, y) ∈ R2 p 1 ln x2 + y 2 . We consider the density function of form ρ(x, y) = and the convolution kernel K(x, y) = − 2π −(x2 +y 2 )/α2 e with α > 0, and the corresponding analytical solution of this problem is given by   2   α2 x + y2 2 2 ϕ(x, y) = − −Ei − + ln(x + y ) , satifying − ∆ϕ = ρ or ϕ = K ∗ ρ on R2 . 4 α2 Rr 2 Here, Ei(r) = −∞ t−1 et dt. Note that the limiting value at the origin is finite: ϕ(0, 0) = α4 (γE − 2 ln α) with γE ≈ 0.5772156649 the Euler–Mascheroni constant. Our goal is to learn the operator G : ρ 7→ ϕ and the computational domain is set as Ω = (−I, I)2 with I = 12. 13









í



í



í

/RVV



D

6WDJH 6WDJH 6WDJH

6WDJH



6WDJH



6WDJH



6WDJH 6WDJH

6WDJH



,WHUDWLRQV

6WDJH





√ Figure 1: Results of DCNO in Algorithm 1 for Example 5.1: (a) Change of Loss= βs−1 Losss during the training; (b) Mean relative L2 error with respect to the number of stages on test set. Table 1: Relative L2 -errors obtained by DCNO and DGNO under different network architectures.

L=1 L=2 L=3 L=4

W0 = 10 1.75 e−3 1.94 e−5 4.01 e−6 1.61 e−6

DCNO W0 = 30 W0 = 50 1.81 e−3 1.20 e−3 1.80 e−5 1.15 e−5 3.17 e−6 1.48 e−6 1.24 e−6 7.04 e−7

12



W0 = 70 5.37 e−4 9.46 e−6 8.29 e−7 5.47 e−7

W0 = 10 2.99 e−3 5.80 e−5 1.68 e−5 1.43 e−5

12

í

í

y í

−12 −12

í

x

D

12 ×10−5

12





y

W0 = 70 6.70 e−4 1.62 e−5 9.20 e−7 6.17 e−7

12 í

í

í

y

DGNO W0 = 30 W0 = 50 2.15 e−3 1.27 e−3 3.17 e−5 2.38 e−5 6.03 e−6 3.47 e−6 3.10 e−6 1.89 e−6

y

í

í −12 −12

í x

E

12 ×10−6

12

í

−12 −12

x

F

12 ×10−5

12







y

 

y

 

−12 −12

x

D

12

 −12 −12

x

E

12

−12 −12



x

F

12

Figure 2: Numerical solution (1st row) and pointwise error (2nd row) of DCNO for Example 5.1 on three test samples (a) α = 1; (b) α = 3; (c) α = 5.

We first test the performance of DCNO in the supervised scenario by applying Algorithm 1 to Example b 5.1. For training, a dataset {(ρ(b) , ϕ(b) )}N b=1 is generated with samples Nb = 100, where each parameter (b) (b) α of ρ is independently sampled from the uniform distribution U[1, 5]. Motivated by the strategy in [2, 4, 45, 46], we take the width parameters W and p of the DeepONet (2.6) in a stage-dependent manner: W = W0 + 20 (s − 1) and p = 120 + 20 (s − 1), where s denotes the stage index. The base width is selected as W0 ∈ {10, 30, 50, 70}, and the network depth is chosen as L ∈ {1, 2, 3, 4}. The input coordinates of 14

β0 θ1 (β)

β1 θ (β)

β2 θ (β)

2

3

×10

í

í

í

í í

í

í

í

í

í



y





í í









 









í í







−2



y

x

β3 θ (β) 4



 í í

×10





y

x

β4 θ (β) 5

−3



 í í

×10





−5











x



Hí



í



í









y



 í í





y

x





 í í

x

Figure 3: Basis functions generated by the basis operators of DCNO method.

the trunk network are sampled from a uniform spatial grid defined by xk = −I + 2kI/Nx , k = 1, . . . , Nx , yk = −I + 2kI/Ny , k = 1, . . . , Ny , where Nx = Ny = 100 is used during training and a finer grid with Nx = Ny = 150 is adopted for testing error. The activation function is chosen as σ = sin, and the learning rate is set to 5 × 10−4 . With L = 4 and W0 = 70 fixed, the stopping criterion in Algorithm 1 is set to √ ϵ = 2 × 10−6 . Fig. 1(a) shows the evolution of βs−1 Losss during training at each stage, with a decay observed as the number of training iterations increases. For testing, 50 samples disjoint from the training set are randomly selected. The corresponding solutions are predicted by DCNO, and the performance is evaluated by the mean relative L2 error over all samples in the test set. The errors of the numerical solution from DCNO under different depth and width are presented in Table 1. It can be observed that increasing the network width and depth consistently reduces the error, indicating improved approximation of the numerical solution. With L = 4 and W0 = 70 fixed, Fig. 1(b) illustrates a clear decrease in the test error as the number of basis operators increases, reaching eventually 5.47 × 10−7 which is near the machine precision under single-precision floating-point arithmetic. Notably, without the progressive construction of basis operators employed in DCNO, the accuracy at stage 1 corresponds to that of the standard DeepONet, which is only 8.17 × 10−3 . This demonstrates the significant advantage of DCNO in improving the accuracy of operator learning. Fig. 2 shows the DCNO numerical solutions for several test samples, along with the corresponding pointwise error distributions. As observed, DCNO method can accurately predict the potential ϕ across the entire computational domain. For a fixed 2 2 ρ(x, y) = e−(x +y )/4 , Fig. 3 visualizes the basis functions {Gθj (ρ)}5j=1 generated by the basis operators. It is observed that, as the stage index s increases, the basis functions tend to capture higher-frequency components. Also, we test the robustness of DCNO in the presence of noisy training data here to validate the theoretical result (3.10). Independent Gaussian noises with relative noise levels ranging from 10−6 to 10−1 are added to the exact solution data, which serve as training labels. With L = 4 and W0 = 70 fixed and the total number of stages set to S = 5, all other training settings remain the same as before. The left plot of Fig. 4 presents the test error of DCNO versus the noise level. The results indicate that the prediction error increases monotonically with the noise amplitude, while remaining consistently below the noise level across all tested regimes spanning several orders. This demonstrates that DCNO achieves robust operator 15

0HDQUHODWLYHL 2 HUURU

×10−3



í



í





í





í



í



í

12

y



(UURURI'&12 1RLVHOHYHO 

í



í

í

í

  1RLVHOHYHO



í



×10−4

12

 

y



 í

−12 −12

x

12

 −12 −12

x

12

/RVV /RVV

Figure 4: Relative L2 error of DCNO under varying noise levels (left); A Gaussian noise sample with noise level 0.001 (middle); Pointwise error of DCNO for Example 5.1 on a test sample with α = 1 (right).



6WDJH

/RVV

í 

6WDJH

í 

6WDJH

/RVV /RVV

í í í í í

D

6WDJH 6WDJH 







,WHUDWLRQV





Figure 5: Results of DGNO in Algorithm 2 for Example 5.1: (a) Change of Loss during the iterations; (b) Mean relative L2 error with respect to the number of stages. Table 2: Number of stages, epochs per stage, training time and error for DCNO and DGNO.

supervised DCNO unsupervised DGNO

Stages 5 5

Epochs per stage 3000 3000

Training time(s) 251.23s 491.60s

Mean relative L2 error 1.31E-6 6.17E-7

approximation under noisy supervision. For illustration, a Gaussian noise with level 0.001 is shown in the middle plot of Fig. 4, together with the corresponding pointwise error of the predicted solution for the test 2 2 sample ρ(x, y) = e−(x +y ) in the right plot. The outcomes in total support Theorem 3.2. We move on to test DGNO on Example 5.1. The unsupervised DGNO needs only the input functions b as the training set, which is adopted as the same {ρ(b) }N b=1 with Nb = 100 for training in DCNO. We set γ = 2 and m = 4 in the combined moment-matching and zero boundary conditions are encoded Pstrategy, p into DeepONet via Gθs (f, x, y) = (I 2 − x2 )(I 2 − y 2 ) k=1 bk (f ) tk (x, y). The trunk input is discretized on Gauss–Legendre nodes with Nx = Ny = 110 in each spatial direction. The total number of stages set to S = 5, and all other training settings are kept the same as for DCNO above. Fig. 5(a) shows the evolution of the loss function during training. In each stage, the loss decreases rapidly as the training iterations proceed. In the testing phase, the same test set as that used for DCNO is adopted. Table 1 presents the errors obtained by DGNO for various network depths and widths. The results indicate a clear reduction in error with increasing network depth and width. With L = 4 and W0 = 70 fixed, Fig. 5(b) demonstrates the reduction in the test error as the stage increases, with the final error 6.17 × 10−7 also approaching the machine precision of single float. These results support Theorem 4.6. Lastly, we compare the performance of DCNO and DGNO on Example 5.1. The results are summarized in Table 2, which presents the number of stages, epochs per stage, training time, and mean relative L2 error for both approaches. It can be observed that, with the same number of stages and training epochs, DGNO 16

achieves a lower generalization error, albeit at the cost of longer training time, which is due to the additional computational overhead of evaluating derivatives in its loss function. Overall, DGNO benefiting from its unsupervised feature is certainly more efficient. Nevertheless, there are convolution problems that do not enjoy a PDE formulation, where DGNO is not applicable and DCNO becomes the solution. The following example is presented to illustrate this. √ Example 5.2 (Multiquadrics). Consider the multiquadric kernels K(x) = 1/ x2 + a2 with 0 < a ≪ 1, commonly used in electrostatic computation [50, 51, 52], for a 1D convolution: Z 1 K(x − x̃)ρ(x̃) dx̃,

ϕ(x) =

x ∈ Ω = (0, 1).

0

It does not have a PDE formulation [51, 52], and our input is the Gaussian-type density function ρ(x) = 1 2 e−β (x− 2 ) with β > 0. In this example, the parameter a > 0 controls the regularization scale of the kernel, and as a → 0, the 1 2 kernel approaches a weakly singular form. For a fixed ρ(x) = e−2(x− 2 ) , the corresponding solutions ϕ(x) for different a are shown in the left plot of Fig. 6. As a decreases, the kernel becomes more singular at 0, resulting in solutions with steeper variations near the boundary. We then set a = 0.0001 to examine DCNO in this near-singular regime. The training data is provided by the traditional method SOE [50], (b) b generating a training set {(ρ(b) , ϕ(b) )}N is generated by uniformly b=1 containing Nb = 50 samples, where ρ sampling β ∼ U[1, 5]. The trunk network is evaluated at collocation points defined on the uniform grid xk = k/Nx , k = 0, . . . , Nx , with Nx = 200 for training and Nx = 300 for testing. The total number of stages is set to S = 4, with σ = tanh and the depth of DeepONet fixed at L = 4. The network width and p are increased with the stage index s, i.e., W = 70 + 10 (s − 1) and p = 120 + 10 (s − 1). For test, 100 samples that are disjoint from the training set are randomly selected to form the test set. The right plot of Fig. 6 shows the mean relative L2 error on the test set, which decreases steadily as the number of basis operators 1 2 increases with the stage index s, and eventually reaches 1.16 × 10−6 . For the test sample ρ(x) = e−2(x− 2 ) , Fig. 7(a) presents the exact and numerical solutions, while the corresponding pointwise error is shown in Fig. 7(b). The results demonstrate that DCNO can achieve high accuracy even for this near-singular setting. The mean relative L2 errors of the approximate solution ϕ over the range of a values are presented in Table 3 to evaluate the robustness of DCNO. The results show that DCNO maintains consistently high accuracy across all tested values of a, indicating that DCNO is uniformly effective from the near-singular (a = 10−4 ) regime to the regular (a = 10−1 ) regime. 

ϕ = 0.0001 ϕ = 0.001

ϕ = 0.01 ϕ = 0.1

ϕ(x)

   0

x

1

Figure 6: Solutions for different values of the kernel parameter a (left); Mean relative L2 error with respect to the number of stages for a = 0.0001 (right).

By Example 5.2, we assess the computational efficiency of our operator learning by comparing the traditional SOE method with DCNO on a test set of 1000 samples. Both methods are evaluated on the same uniform grid with 64 spatial points for a = 0.0001. The results in Table 4 indicate that while DCNO incurs 17

Hí





D

F

'&12WUDLQLQJHYDOXDWLRQ '&12HYDOXDWLRQ

7LPHV

(UURU

ϕ(x)

 

7UDGLWLRQDOVROYHU

E





([DFWVROXWLRQ 1XPHULFDOVROXWLRQ



0

x

1

í

0

 

1

x





1XPEHURIFRPSXWDWLRQDOFDVHV

Figure 7: Results of DCNO for Example 5.2: (a) Exact and numerical solutions; (b) Error of the numerical solution; (c) Computational time of the traditional solver and DCNO against the number of computational cases. Table 3: Mean relative L2 error of numerical solutions from DCNO for varying a.

a Error

0.0001 1.20 e−6

0.001 3.08 e−6

0.001 2.71 e−6

0.1 3.48 e−6

Table 4: Comparison between DCNO and the traditional method SOE on test set.

Method Traditional method DCNO method

Evaluation time (s) 670.44 0.021

Training time (s) — 202.14

additional training overhead, it drastically reduces the evaluation time compared to the traditional approach, achieving several orders of magnitude speedup during inference and demonstrating clear advantages for large-scale evaluations. Fig. 7(c) further reports the total computational time of SOE and DCNO as the number of test cases increases. DCNO demonstrates a much slower growth rate, highlighting its superior computational efficiency at scale. So will DGNO if applicable. Example 5.3 (Dynamics of Schrödinger-Poisson system). For quantum mechanics taking into account selfgravitating interaction, one considers the nonlinear Schrödinger-Newton/Poisson equations that reads in dimensionless form [39, 49]: 1 i∂t ψ(x, t) = − ∆ψ + V (x)ψ + ϕ(x, t)ψ + |ψ|2 ψ, 2 − ∆ϕ(x, t) = |ψ|2 , x ∈ Ω, t > 0, ψ(x, 0) = ψ0 (x),

x ∈ Ω;

x ∈ Ω ⊂ R2 , t > 0,

ψ(x, t) = 0, ϕ(x, t) = 0,

(5.3a) (5.3b)

x ∈ ∂Ω, t ≥ 0.

(5.3c)

Here, ψ is a complex-valued wave function, ϕ is the real-valued gravitational potential, V (x) = 21 (x2 + y 2 ) denotes an external harmonic oscillator trapping potential. We consider the circular computational domain Ω = {x ∈ R2 : |x| < I} with radius I = 8 and solve its dynamics for t ∈ [0, 1]. The initial condition is given 2 2 by ψ0 (x, y) = √12π e−(x +y )/4 . A popular numerical solver for the dynamics of (5.3) on a circular domain is the time-splitting finite element method, briefly described below as a reference. The spatial domain is discretized by a conforming triangular mesh with Nh interior nodes, where continuous piecewise linear (P1 ) finite elements are employed. h Let {ϑj }N j=1 denote the corresponding nodal basis functions. The finite element approximation of ψ(·, tn ) PNh n n T n n ∈ CNh . Given the time discretization is given by ψh (·, tn ) = j=1 Ψj ϑj , with Ψ = (Ψ1 , . . . , ΨNh )

18

 7UDGLWLRQDOVROYHU

F

E

 '*12



0DVV

7LPHV























 





t







t

Figure 8: Results of Example 5.3: (a) Relative L2 error of ψ over time; (b) Total evaluation time of the traditional solver and DGNO for solving the Poisson equation (5.3b) during the simulation; (c) Evolution of the mass over time. W 

W 

×í

4 0.3

0.3

2

2

4

y

x 4−4

0

ψ|

0.2

0.1

−4



|

ψ|

| 0.2

0.1

−4

4

y

x 4−4

y



0

−4−4

x

4



Figure 9: |ψ(·, t = 0)|2 (left); |ψ(·, t = 1)|2 by DGNO (middle); Pointwise absolute error at t = 1 (right) for Example 5.3.

tn = n∆t, n ≥ 0, the Strang splitting Crank–Nicolson finite element scheme [6] reads     i∆t i∆t ∗ Mh + Kh Ψ = Mh − Kh Ψn , 8 8 
 Ψ∗∗ = exp −i∆t V + Φ∗ + |Ψ∗ |2 Ψ∗ ,     i∆t i∆t n+1 Kh Ψ = Mh − Kh Ψ∗∗ . Mh + 8 8

(5.4a) (5.4b) (5.4c)

Here, Mh and Kh denote the finite element mass and stiffness matrices, respectively, with entries (Mh )jk = ⟨ϑj , ϑk ⟩L2 (Ω) , (Kh )jk = ⟨∇ϑj , ∇ϑk ⟩L2 (Ω) . The numerical potential Φ∗ ∈ RNh is obtained from the finite element discretization of the Poisson equation (5.3b), leading to the linear system Kh Φ∗ = Fh (|Ψ∗ |2 ), where Fh (|Ψ∗ |2 ) denotes the assembled load vector. Since the stiffness matrix Kh is time-independent, its sparse LU factorization is performed only once before time stepping, while each Poisson solve requires only triangular substitutions. Nevertheless, solving the Poisson equation repeatedly at every time step remains computationally demanding, especially for small ∆t or long-time simulations. Here, we employ our DGNO method to replace the evaluation process for Φ∗ in (5.4). The training of DGNO is performed similarly to that in Example 5.1. After training, only a single neural operator inference is required at each time step to efficiently obtain an approximation of the gravitational potential. We set ∆t = 0.001, Nh = 80000 for the experiment. Fig. 8(a) illustrates the relative L2 error of the solution field ψ over time. The linear drift in the solution error with respect to time comes from the natural accumulation of the approximation error. Fig. 8(b) compares the accumulated evaluation time of the Poisson solver based on the finite element method and DGNO throughout the simulation. The evolution of the mass R conventional 2 |ψ(x, t)| dx over time is presented in Fig. 8(c). Fig. 9 further presents |ψ|2 at t = 0 and t = 1, together Ω with the corresponding pointwise absolute error at t = 1. As can be seen, the proposed method remains stable and accurate throughout the evolution, provides an accelerated solver for (5.3), and exhibits excellent 19

mass conservation. 5.3. Extension to multiple inputs We generalize our methods to further study the convolution problem (1.1) in which both the density function and the convolution kernel vary, forming a multi-input operator learning task. In fact, this can be quite straightforward under our approach by adding additional branches to (2.6), replacing DeepONet with the Multiple-Input Operator Network (MIONet) [24]. |x|2

− 2α2 1 ,x∈Ω= Example 5.4 (DCNO with MIONet). Consider the Gaussian kernel family Kα (x) = 2πα 2e (−8, 8)2 , α > 0. Define the convolution operator Z G : (ρ, Kα ) 7→ ϕ, ϕ(x) = (Kα ∗ ρ)(x) = Kα (x − x̃)ρ(x̃) dx̃. Ω

Such a convolution is repeatedly needed in constructing the scale-space of the scale-invariant feature transform 1 −|x|2 /(2τ 2 ) with parameter algorithm [33]. Here we take ρ from a family of Gaussian densities ρ(x) = 2πτ 2e 1 −|x|2 /(2α2 +2τ 2 ) τ > 0, and then the convolution admits the analytical solution ϕ(x) = 2π(α2 +τ 2 ) e as reference. To solve Example 5.4, we apply the DCNO algorithm with the basis operator network constructed using the MIONet architecture [24] to accommodate multiple input functions. For training, a dataset (l) Nb ,Nl {(ρ(b) , Kα , ϕ(b,l) )}b=1,l=1 with Nb = Nl = 20 samples is generated, with each parameter α(b) and τ (l) is independently drawn from U[1, 2]. All other training settings are the same as those in Example 5.1, except that Nx = Ny = 64 are used during training, while a finer grid with Nx = Ny = 100 is adopted for error evaluation. At each stage s, the practical loss function for DCNO in the multiple-input case is extended to N

Losss (θ) =

N

N

2

x b X l X X 1 1 es−1 (ρ(b) , Kα(l) )(xk ) − Gθ (ρ(b) , Kα(l) , xk ) . Nb × N l × N x βs−1

(5.5)

b=1 l=1 k=1

In the testing phase, 200 samples disjoint from the training set are randomly selected. Fig. 10(a) illustrates a clear decrease in the test error as the number of basis operators increases, and ultimately reaches an accuracy of 1.40 × 10−5 . Notably, without the progressive construction of basis operators employed in DCNO, the accuracy at stage 1 corresponds to that of the standard MIONet, which is only 4.84 × 10−2 . This demonstrates the significant advantage of DCNO in improving the accuracy of multi-input operator learning. Fig. 11 presents the numerical solutions on the test set for different combinations of density functions and convolution kernels, together with the corresponding pointwise error distributions. These results demonstrate that DCNO maintains strong performance in the multiple-input setting for generalizing simultaneously over the density and the convolution kernel.

Figure 10: Mean relative L2 error with respect to the number of stages: (a) Example 5.4; (b) Example 5.5 .

20

8



8



y

x

D



8

×10−6

8

y

−8 −8

x

E



8

×10−7

8

−8 −8





x

D

8

 

−8 −8

x

F



8

×10−6

8







y







−8 −8







y

8



y



y

 −8 −8

x

E

8





−8 −8

x

F

8

Figure 11: Numerical solution (1st row) and pointwise error (2nd row) of DCNO for Example 5.4 on three test samples: (a) α = τ = 2; (b) α = τ = 1.5; (c) α = τ = 1.

Example 5.5 (DGNO with MIONet). We turn to the unsupervised scenario by considering the following general elliptic boundary value problem (
−∇ · κ(x)∇ϕ(x) = ρ(x), x ∈ Ω, (5.6) ϕ(x) = 0, x ∈ ∂Ω, where the computational domain is set as Ω = (0, 1)2 . In (5.6), different choices of κ(x) correspond to different integral kernels in the convolution (1.1). Our goal is to learn the mapping G : (ρ, κ) 7→ ϕ with ρ and κ jointly serving as inputs, which constitutes a multiple-input operator learning problem. The differential formulation in which the problem is posed enables an unsupervised training strategy, and we apply the DGNO algorithm with two branch networks to input ρ and κ. To obtain reference solutions for tests, we set κ(x) = 1 + α sin(2πx) sin(2πy) with |α| < 1 and h i h ρ(x) = 2π 2 α κ(x) γ sin(πx) sin(πy) + 4 sin(2πx) sin(2πy) − 2π 2 α2 γ cos(2πx) sin(2πy) cos(πx) sin(πy)

i + sin(2πx) cos(2πy) sin(πx) cos(πy) + 2 cos2 (2πx) sin2 (2πy) + sin2 (2πx) cos2 (2πy) , where ϕ(x) = αγ sin(πx) sin(πy) + α sin(2πx) sin(2πy) is the exact solution to (5.6). For the training data, a b ,Nl set {(κ(b) , ρ(l) )}N b=1,l=1 with Nb = Nl = 20 is generated, resulting in a total of 400 samples. The parameters α(b) and τ (l) are independently sampled from U (0.5, 1) and U[2, 4], respectively. All other training settings are the same as those in Example 5.1, except that Nx = Ny = 110 are used during training, while a finer grid with Nx = Ny = 150 is adopted for error evaluation. In view of (4.3), the extended practical loss for DGNO in the multiple-input case reads N

Losss (θ) =

N

b X l Rs−1 (κ(b) , ρ(l) , Gθs (κ(b) , ρ(l) )) −1 X . Nb × Nl ∥Gθs (κ(b) , ρ(l) )∥X

(5.7)

b=1 l=1

For testing, 200 samples that are disjoint from the training set are randomly selected to form the test set. As shown in Fig. 10(b), the mean relative L2 error on the test set decreases progressively with the stage 21

1



1

1







y





y



0



0

x

D

1 ×10−4

1



0

0

x

E



1

×10−4

1



0

x

D

1

0

0

x

F

1



y

 

y



 0

0

x

E

1





×10−4

1





0







y

y

0

0

x

F

1



Figure 12: Numerical solution (1st row) and pointwise error (2nd row) of DGNO for Example 5.5 on three test samples: (a) α = 0.5, γ = 2; (b) α = 0.6, γ = 3; (c) α = 0.7, γ = 4.

index s, and ultimately reaches an accuracy of 7.08 × 10−5 . We also select some representative samples from the test set to visually assess the approximation capability of the trained basis operators in Fig. 12. These results demonstrate that DGNO as the unsupervised method also performs well for the multiple-input case. 6. Conclusion In this work, we have introduced two novel and general multi-stage neural operator learning frameworks, the Deep Collocation Neural Operator (DCNO) and the Deep Galerkin Neural Operator (DGNO), designed to significantly enhance the accuracy of operator approximation, particularly applied for convolutions. DCNO offers a supervised approach that progressively builds a richer approximation space by training basis operators on the residuals of the target operator, leveraging available input-output data. Complementing this, DGNO provides an unsupervised framework that learns operator bases by minimizing the residual of a governing PDE in its weak form, thus eliminating the need for labeled data. We have demonstrated the efficacy of these methods through extensive numerical experiments. Both DCNO and DGNO achieved high accuracy in learning convolution operators, often approaching machine precision of single-precision floating-point arithmetic. Notably, for convolution problems lacking a straightforward PDE formulation, DCNO proved to be a robust supervised solution. When a PDE formulation is available, DGNO offers an efficient unsupervised alternative. Furthermore, we have extended these frameworks to handle multi-input operator learning scenarios, encompassing variations in both the density and the kernel for convolutions. Acknowledgements Z. Wen and X. Zhao are supported by National Key Research and Development Program of China, National MCF Energy R&D Program (No. 2024YFE03240400), National Natural Science Foundation of China (Nos. 42450275, 12271413). Z. Mao is supported by the Zhejiang Provincial Natural Science Foundation of China under Grant No. LR26A010001 and the National Natural Science Foundation of China under Grant No. 12531015. Y. Zhang is supported by the National Natural Science Foundation of China (No. 12271400), the National Key R&D Program of China (No. 2024YFA1012803) and basic research fund of Tianjin University under Grant 2025XJ21-0010. 22

References [1] R. A. Adams and J. J. Fournier, Sobolev Spaces, vol. 140, Elsevier, 2003. [2] M. Ainsworth and J. Dong, Galerkin neural networks: A framework for approximating variational equations with error control, SIAM J. Sci. Comput., 43 (2021), pp. A2474–A2501. [3] M. Ainsworth and J. Dong, Extended Galerkin neural network approximation of singular variational problems with error control, SIAM J. Sci. Comput., 47 (2025), pp. C738–C768. [4] Z. Aldirany, R. Cottereau, M. Laforest, and S. Prudhomme, Multi-level neural networks for accurate solutions of boundary-value problems, Comput. Methods Appl. Mech. Engrg., 419 (2024), p. 116666. [5] A. Anandkumar, K. Azizzadenesheli, K. Bhattacharya, N. Kovachki, Z. Li, B. Liu, and A. Stuart, Neural operator: Graph kernel network for partial differential equations, ICLR, (2020). [6] W. Auzinger, T. Kassebacher, O. Koch, and M. Thalhammer, Convergence of a Strang splitting finite element discretization for the Schrödinger–Poisson equation, ESAIM Math. Model. Numer. Anal., 51 (2017), pp. 1245–1278. [7] W. Bao, S. Jiang, Q. Tang, and Y. Zhang, Computing the ground state and dynamics of the nonlinear Schrödinger equation with nonlocal interactions via the nonuniform FFT, J. Comput. Phys., 296 (2015), pp. 72–89. [8] H. Bassi, Y. Zhu, S. Liang, J. Yin, C. C. Reeves, V. Vlček, and C. Yang, Learning nonlinear integral operators via recurrent neural networks and its application in solving integro-differential equations, Mach. Learn. Appl., 15 (2024), p. 100524. [9] S. L. Brunton and J. N. Kutz, Promising directions of machine learning for partial differential equations, Nat. Comput. Sci., 4 (2024), pp. 483–494. [10] S. Cai, Z. Wang, L. Lu, T. A. Zaki, and G. E. Karniadakis, DeepM&Mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks, J. Comput. Phys., 436 (2021), p. 110296. [11] Z. Chang, B. Gao, R. Yin, and X. Zhao, DLTM: A deep learning method for tearing mode simulation and prediction, Nucl. Fusion, 65 (2025), p. 086009. [12] Z. Chang, Z. Wen, and X. Zhao, Unsupervised operator learning approach for dissipative equations via Onsager principle, SIAM J. Sci. Comput., 48 (2026), pp. C1060–C1085. [13] T. Chen and H. Chen, Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems, IEEE Trans. Neural Netw., 6 (1995), pp. 911–917. [14] P. C. Di Leoni, L. Lu, C. Meneveau, G. E. Karniadakis, and T. A. Zaki, Neural operator prediction of linear instability waves in high-speed boundary layers, J. Comput. Phys., 474 (2023), p. 111793. [15] L. Exl, N. J. Mauser, and Y. Zhang, Accurate and efficient computation of nonlocal potentials based on Gaussian-sum approximation, J. Comput. Phys., 327 (2016), pp. 629–642. [16] X. Glorot and Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in AISTATS, JMLR Workshop and Conference Proceedings, 2010, pp. 249–256. [17] R. C. Gonzalez, Digital Image Processing, Pearson Education India, 2009. 23

[18] J. Y. Greengard, L.and Lee, Accelerating the nonuniform fast Fourier transform, SIAM Rev., 46 (2004), p. 443–454. [19] L. Greengard, The Rapid Evaluation of Potential Fields in Particle Systems, MIT Press, 1988. [20] W. Hackbusch, A sparse matrix arithmetic based on H-matrices. Part I: Introduction to H-matrices, Computing, 62 (1999), pp. 89–108. [21] W. Hackbusch et al., Hierarchical Matrices: Algorithms and Analysis, vol. 49, Springer, 2015. [22] K. Hornik, Approximation capabilities of multilayer feedforward networks, Neural Netw., 4 (1991), pp. 251–257. [23] Y. Ji, J. Liang, and Z. Xu, Machine-learning interatomic potentials for long-range systems, Phys. Rev. Lett., 135 (2025), p. 178001. [24] P. Jin, S. Meng, and L. Lu, MIONet: Learning multiple-input operators via tensor product, SIAM J. Sci. Comput., 44 (2022), pp. A3490–A3514. [25] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, ICLR, (2015), p. 13. [26] H. Li and M. Shatarah, Operator learning for urban water clarification hydrodynamics and particulate matter transport with physics-informed neural networks, Water Res., 251 (2024), p. 121123. [27] Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar, Fourier neural operator with learned deformations for PDEs on general geometries, J. Mach. Learn. Res., 24 (2023), pp. 1–26. [28] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, ICLR, (2021). [29] C. Lin, Z. Li, L. Lu, S. Cai, M. Maxey, and G. E. Karniadakis, Operator learning for predicting multiscale bubble growth dynamics, J. Chem. Phys., 154 (2021). [30] X. Liu, Q. Tang, S. Zhang, and Y. Zhang, On optimal zero-padding of kernel truncation method, SIAM J. Sci. Comput., 46 (2024), pp. A23–A49. [31] X. Liu, Q. Tang, and Y. Zhang, Fast convolution solvers using moment-matching, arXiv preprint arXiv:2602.12850, (2026). [32] Z. Liu, H. Wang, H. Zhang, K. Bao, X. Qian, and S. Song, Render unto numerics: Orthogonal polynomial neural operator for PDEs with nonperiodic boundary conditions, SIAM J. Sci. Comput., 46 (2024), pp. C323–C348. [33] D. G. Lowe, Distinctive image features from scale-invariant keypoints, Int. J. Comput. Vis., 60 (2004), pp. 91–110. [34] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nat. Mach. Intell., 3 (2021), pp. 218–229. [35] L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, and G. E. Karniadakis, A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data, Comput. Methods Appl. Mech. Engrg., 393 (2022), p. 114778. [36] Z. Mao, L. Lu, O. Marxen, T. A. Zaki, and G. E. Karniadakis, DeepM&Mnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators, J. Comput. Phys., 447 (2021), p. 110698. 24

[37] N. H. Nelsen and A. M. Stuart, The random feature model for input-output maps between Banach spaces, SIAM J. Sci. Comput., 43 (2021), pp. A3212–A3243. [38] A. V. Oppenheim, Discrete-Time Signal Processing, Pearson Education India, 1999. [39] R. Penrose, On the gravitization of quantum mechanics 1: Quantum state reduction, Found. Phys., 44 (2014), pp. 557–575. [40] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys., 378 (2019), pp. 686–707. [41] F. Rosenblat, The perceptron: A probabilistic model for information storage and organization in the brain, Psychol. Rev., 65 (1958), pp. 386–408. [42] F. Vico, L. Greengard, and M. Ferrando, Fast convolution with free-space Green’s functions, J. Comput. Phys., 323 (2016), pp. 191–203. [43] H. Wang, T. Fu, Y. Du, W. Gao, K. Huang, Z. Liu, P. Chandak, S. Liu, P. Van Katwyk, A. Deac, et al., Scientific discovery in the age of artificial intelligence, Nature, 620 (2023), pp. 47–60. [44] S. Wang, H. Wang, and P. Perdikaris, Learning the solution operator of parametric partial differential equations with physics-informed DeepONets, Sci. Adv., 7 (2021), p. eabi8605. [45] Y. Wang and C. Lai, Multi-stage neural networks: Function approximator of machine precision, J. Comput. Phys., 504 (2024), p. 112865. [46] M. Weng, Z. Mao, and J. Shen, Deep collocation method: A framework for solving PDEs using neural networks with error control, SIAM J. Sci. Comput., 48 (2026), pp. C77–C102. [47] Y. Xu, Multi-grade deep learning, Communications on Appl. Math. Comput., 8 (2026), pp. 778–829. [48] L. Ying, G. Biros, and D. Zorin, A kernel-independent adaptive fast multipole algorithm in two and three dimensions, J. Comput. Phys., 196 (2004), pp. 591–626. [49] Y. Zhang and X. Dong, On the computation of ground state and dynamics of Schrödinger–Poisson– Slater system, J. Comput. Phys., 230 (2011), pp. 2660–2676. [50] Y. Zhang, C. Zhuang, and S. Jiang, Fast one-dimensional convolution with general kernels using sum-of-exponential approximation, Commun. Comput. Phys., 29 (2021), pp. 1570–1582. [51] C. Zhuang and R. Zeng, A local discontinuous Galerkin method for 1.5-dimensional streamer discharge simulations, Appl. Math. Comput., 219 (2013), pp. 9925–9934. [52] C. Zhuang, Y. Zhang, X. Zhou, R. Zeng, J. He, and L. Liu, A fast tree algorithm for electric field calculation in electrical discharge simulations, IEEE Trans. Magn., 54 (2018), pp. 1–4.

25

Related documents

• Multiscale Training of Convolutional Neural Networks#625463
• Multiscale Training of Convolutional Neural Networks#623190
• MRI-based brain tumor prediction using convolutional neural network framework.#221076
• Very Deep Convolutional Networks for Large-Scale Image Recognition#6258
• Hybrid Dealiasing and Implicit Packing for Real Convolutions#609658