PaperlayerView exact revision on arXiv

Source-linked AI summary

Numerical Neural Operator: Connection Between Numerical Analysis and the Operator Learning and Building Discretized Neural Operators from Finite Element Methods

Zecheng Zhang

arXiv:2610.11120v1math.NAcs.NA

TL;DR

The paper asks whether neural operators trained on numerical data approximate the continuous PDE operator or the discretized solver map, a connection that remains poorly understood. It derives neural operators from numerical schemes and analyzes their convergence. In numerical tests, the global-basis model has the lowest mean relative errors across three PDEs while using fewer parameters than DeepONet.

  • Problem

    It remains unclear whether neural operators trained on solver-generated data approximate the continuous PDE solution operator or the numerical scheme’s discrete map.

    Source: §1 · p. 2

  • Method

    The paper derives network blocks from numerical discretizations, using neural representations of FEM basis functions and time-unrolled coefficient updates.

    Source: §1 · p. 3, §4.5 · p. 14

  • Results

    Across three PDEs, the proposed global-basis model achieves the lowest mean relative errors—0.08%, 0.28%, and 0.54%—while using approximately 16% fewer parameters than DeepONet.

    Source: §5 · p. 23

  • Takeaways & Limitations

    The approach provides a numerical-scheme-based network design that supports different basis representations and offers a practical tradeoff between accuracy and model size.

    Source: §5 · p. 23

    The reported width bound is not a parameter-count bound, since dense affine matrices may contain many entries.

    Source: §4.5 · p. 15

Abstract

from arXiv · show

Operator learning aims to learn mappings between function spaces and has been widely used for PDE-related problems. This raises two fundamental questions. First, is there a connection between standard numerical methods and neural operators, given that both are designed to approximate PDE solution operators? Second, when a neural operator is trained on numerically generated data, does it approximate the underlying infinite dimension operator, the training data themselves, or the discretized numerical scheme used to generate them? A key hypothesis in the paper is that the trained neural operator may mimic the discretized numerical schemes that generate the training data. The hypothesis has been partially and preliminary verified by the linear approximation principle proposed in Chen Chen 1995 paper and the following scaling law papers. However, neural operator approximation from the perspective of numerical discretization remains unexplored. In this work, we consider a class of PDEs and derive neural operator architectures from the numerical schemes. Firstly, the numerical basis functions, such as finite element basis functions, can be approximated or directly constructed by learnable neural networks. Then the finite difference schemes for temporal discretization can be unrolled into iterative neural architectures that generate the coefficients of the basis of the resulting neural operator. By leveraging numerical analysis, we establish convergence error estimates for these numerically derived neural operators. Notably, the network size depends polynomial on the approximation error in $L_2$, improving the nested exponential dependence derived in the general neural operator approximation and aligning with the numerical observations. Finally, numerical experiments demonstrate reduced errors and improved parameter efficiency.

1 Introduction

The introduction asks whether neural operators trained on solver-generated data approximate the continuous PDE operator or the discretized numerical scheme. It proposes deriving operator architectures from numerical discretizations and analyzing their convergence.

1 Introduction

  • Operator learning also serves PDE prediction and multi-fidelity tasks, including mapping initial conditions to solutions and reconstructing fine-scale solutions from coarse inputs.The introduction describes prediction for new initial conditions and reconstruction using limited fine-scale observations.Source: §1 · pp. 1–2
  • Most scientific training data come from numerical solvers, leaving unclear whether trained neural operators approximate the continuous solution operator or the solver’s discrete input-output map.The paper suggests the latter perspective may help explain observed approximation-rate behavior, alongside low-dimensional structure.Source: §1 · p. 2
  • The connection between classical FEM and neural operators remains incomplete because FEM solves for coefficients algebraically while neural operators generate them with trainable networks.How spatial and temporal discretizations jointly shape neural-operator architectures is also largely unexplored.Source: §1 · p. 2
  • The paper reverses the usual design direction by deriving neural operators from numerical discretizations, then develops convergence estimates relating L2 error to component sizes.The introduction says the resulting rates recover convergence behavior commonly seen in numerical computations.Source: §1 · p. 3

2 Preliminary

The preliminary section sets out ReLU approximation tools and the general operator-learning framework. It also contrasts general-theory rates with faster convergence commonly observed in numerical experiments.

2 Preliminary

  • ReLU networks can represent clipping exactly, approximate multiplication and squares, and approximate bounded Sobolev functions with stated depth, width, and size bounds.For the Sobolev approximation result, the size bound scales as ε^-d/n(1 + log(1/ε)).Source: §2 · pp. 3–4

2.1 Neural operator

  • The general operator-learning framework approximates continuous operators using neural networks, but its theoretical rates are slower than the power-law scaling commonly observed in numerical experiments.The background notes that specific PDE structure and low-dimensional input parametrizations may contribute to the observed faster rates.Source: §2.1 · p. 4

3 Methodology

The methodology derives a neural operator by discretizing a semilinear parabolic PDE in space and time, then expressing each numerical update as a network block. The resulting architecture combines FEM basis functions with time-unrolled coefficient maps.

3 Methodology

  • The methodology uses a weak-form discretization and demonstrates temporal discretization with forward Euler, while noting that other time schemes can be used similarly.Forward Euler is chosen because its conditional stability links mesh size h and time step Δt for the later error analysis.Source: §3 · p. 4

3.1 General semilinear parabolic problem and weak form

  • The target problem is a semilinear parabolic PDE with zero boundary data, and its weak form follows by testing with a function in H1_0(Ω) and integrating by parts.The stated assumptions include f ∈ L2(0,T; C(Ω)) and locally Lipschitz g.Source: §3.1 · p. 5

3.2 Space and temporal discretizations.

  • The spatial discretization uses conforming finite elements, while forward Euler advances coefficient vectors under a CFL stability condition for the analysis.The analysis chooses Δt ∼ h^2, giving Nt ∼ h^-2; the text says the CFL condition is unnecessary for network design with an unconditionally stable scheme.Source: §3.2 · pp. 5–6

3.3 Quadrature representation for the nonlinear terms and FEM solution

  • Quadrature expresses a time step as an affine evaluation map, componentwise nonlinearity, affine recombination, residual linear branch, and source branch.In the described block, B evaluates the finite-element function at quadrature points and C_h recombines the nonlinear values.Source: §3.3 · pp. 6–7

3.3.1 Network structure

  • The finite-element basis supplies output basis functions, while the time-unrolled numerical scheme supplies coefficient maps realized with learnable layers and neural approximations.P1 basis functions can be represented exactly by ReLU networks; higher-order basis functions are approximated, and the shared coefficient network generates all coefficient outputs.Source: §3.3.1 · pp. 7–9

4 Fisher–KPP Study and Analysis

The Fisher–KPP analysis specializes the update to the reaction term v − v^2 and uses quadrature to express its nonlinear contribution. It assumes sufficiently exact quadrature so the mass matrix and square reaction term are represented consistently.

4 Fisher–KPP Study and Analysis

  • For Fisher–KPP, g(v) = v − v² and f = 0; the method separates the linear term v, defines u_h^n as the fully discrete coefficient vector at t_n = nΔt, and gives a forward Euler update.Source: §4 · p. 9

4.1 Quadrature in the Fisher–KPP update

  • Degree-3r exact quadrature makes the mass matrix and Fisher–KPP square term exact; the analysis does not pursue the non-mass-exact case.Without mass exactness, the linear reaction contribution uses M_Q u rather than M u and adds a quadrature consistency term.Source: §4.1 · pp. 9–10

4.2 Preparation for analysis: norms and operator bounds

The analysis uses mass- and quadrature-weighted norms to measure coefficient and quadrature-value vectors. It then bounds the operators used in one time step under the stated stability conditions.

4.2 Preparation for analysis: norms and operator bounds

  • The mass-matrix norm of a coefficient vector equals the L2(Ω) norm of its finite element function.The quadrature-weight norm is the analogous L2-type norm for values at quadrature points.Source: §4.2 · p. 10
  • The analysis defines induced norms for Lh, Ch, and B using the mass and quadrature norms.These norms match the two maps used in one time step.Source: §4.2 · p. 11
  • The bounds are ∥B∥M→W = 1 and ∥Ch∥W→M ≤ ∆t; under the explicit diffusion CFL condition, ∥Lh∥M→M ≤ 1 + ∆t.The CFL condition is ∆tλmax(M−1S) ≤ 2.Source: §4.2 · p. 11

4.3 Predicted coefficient trajectory error

The coefficient-error analysis compares network iterates with a reference finite element trajectory. Under numerical boundedness and the stated scheme assumptions, it controls the accumulated error from approximating the clipped-square update.

4.3 Predicted coefficient trajectory error

  • The perturbation argument assumes the reference solution is bounded at quadrature points, but does not require a pointwise bound on network iterates.The clipped-square map is globally Lipschitz, allowing the perturbation argument to remain valid even when network iterates are unbounded.Source: §4.3 · p. 11
  • The clipping level M0 is chosen at least M*, independently of h, ∆t, and the initial condition.This ensures clipping does not change reference solution values at quadrature points.Source: §4.3 · pp. 11–12
  • Reference boundedness is an explicit assumption because it is not automatic for consistent-mass Galerkin schemes under the diffusion CFL condition alone.A method preserving [0, 1] at quadrature points yields M* = 1 for admissible initial data.Source: §4.3 · p. 11
  • The accumulated coefficient error is bounded by CTδ under the stated quadrature, CFL, and reference-boundedness assumptions.CT is independent of δ and ∆t; the estimate transfers to the finite element function through mass-norm and L2 equivalence.Source: §4.3 · pp. 12–13

4.4 Numerical error

The numerical-error analysis separates finite element approximation, initial encoding, and quadrature contributions. Under standard regularity, stability, and reaction-control assumptions, it invokes a semilinear parabolic finite element estimate.

4.4 Numerical error

  • The numerical estimate assumes the initial finite element encoding has L2 error at most Cinithm, uniformly over admissible initial data.This initial approximation error is included in the spatial error term.Source: §4.4 · p. 13
  • For Pr elements and sufficiently smooth solutions, the spatial order m is typically r + 1.The estimate also includes quadrature consistency error ηQ.Source: §4.4 · p. 13
  • For the Fisher–KPP square case on affine elements, quadrature exact through degree 3r makes ηQ = 0.With exact square quadrature, the text then specializes the estimate under the CFL scaling ∆t ≲ h2.Source: §4.4 · p. 13

4.5 Semi-network approximation error

The semi-network combines a time-unrolled coefficient network with finite element basis networks. Under the theorem’s assumptions, its approximation error and network scaling are controlled by the spatial and time discretization.

4.5 Semi-network approximation error

  • For continuous P1 basis networks, the construction introduces no basis-realization error, so the neural operator is the finite element reconstruction of network-generated coefficients.The stated error estimate combines classical finite element error with coefficient-network error.Source: §4.5 · pp. 13, 15
  • Continuous P1 nodal basis functions on locally convex simplicial meshes have exact ReLU representations with depth at most C(1 + log2 kh) and width at most Ckh.The theorem’s P1 setting uses affine, shape-regular, quasi-uniform, locally convex simplicial meshes.Source: §4.5 · pp. 13–14
  • The exact P1 basis realization does not directly apply to Q1 elements in dimensions d ≥ 2, whose bases are generally bilinear or multilinear.Q1 bases therefore require an additional approximation argument, except in one dimension where Q1 = P1.Source: §4.5 · p. 14
  • With h ∼ ε^(1/η), ∆t ∼ h2, and δ ∼ ε, N ∼ ε^(−d/η), Nt ∼ h^−2, and the L2(Ω) error is ≲ ε.Here η = min{m, 2}; the theorem assumes a classical fully discrete P1 finite element solution with the stated uniform accuracy.Source: §4.5 · p. 14
  • Each time-step block applies Lh u + ChRδ(Bu), with the scalar clipped-square approximation as the coefficient network’s only internal approximation.Lh, Ch, and B are exactly representable as affine ReLU layers, and stacking Nt blocks produces the coefficient vector.Source: §4.5 · p. 14
  • The coefficient network’s depth scales as ≲ h^−2(1 + log(M2/δ)), while its width scales as ≲ h^−d.The depth reflects Nt ∼ h^−2 time-step blocks; the width estimate counts maximum hidden-layer width.Source: §4.5 · p. 15
  • The width bound is not a parameter-count bound: dense affine matrices may contain many entries unless sparsity or structured implementations are exploited.The text lists locality, mass lumping, iterative solvers, and structured implementations as possible ways to exploit structure.Source: §4.5 · p. 15

4.6 Full approximation with nonlocal approximation

This section derives a full operator approximation using neural coefficient and basis networks in a finite element setting. Its error analysis combines numerical-scheme and network approximation errors under the stated assumptions.

4.6 Full approximation with nonlocal approximation

  • Nonlocalized basis networks approximate each finite element basis to accuracy γ, but need not preserve its support.The analysis bounds the basis-network contribution without requiring support preservation.Source: §4.6 · pp. 16–17
  • The coefficient network iterates a learned block, while basis networks reconstruct the finite element output from the resulting coefficients.The finite element encoder supplies the same initial coefficient vector to the numerical and network schemes.Source: §4.6 · pp. 16–18
  • The coefficient network is formed by iterating the prescribed block, and its depth and width are bounded in terms of δ and η.The section also notes that parallel concatenation of the basis networks dominates the resulting width.Source: §4.6 · pp. 17–19
  • With the stated CFL condition and parameter choices, the approximation error is O(ε) in L2.The choices include h ∼ ε^1/η and γ = ε^1+d/η, with η = min{m, 2}.Source: §4.6 · p. 18

4.7 Full approximation case with localized basis networks

This section constructs localized neural approximations of finite element bases and uses them in a full operator approximation. Localization preserves the final L2 error order while reducing the basis-approximation burden and, under the stated construction, total basis-network size.

4.7.1 Strict support-preserving approximation

  • The strict-support approximation uses a local mask and a zero-preserving product subnet to enforce support while retaining the usual network-size order.The localization lemma applies when the multiplier is accurate over a bounded interval containing the raw approximant’s range.Source: §4.7.1 · p. 19

4.7.2 Localized neural approximations of finite element bases

  • Lemma 10 constructs ReLU basis networks with γ approximation error and support contained in Dj,h, although some affine weights may depend on h.The stated size bound is independent of h, j, and γ.Source: §4.7.2 · p. 20

4.7.3 Full operator approximation

  • Theorem 3 gives an L2 approximation error at most Cε for the full operator under its stated finite element, quadrature, boundedness, CFL, and numerical-error assumptions.The coefficient network is obtained by iterating a ReLU block, and the basis networks are localized.Source: §4.7.3 · p. 21
  • Localization retains final L2 error order ε, removes the h^-d factor from the basis-error bound, and saves a factor ε^-d2/η in total basis-network size under the common construction.The required basis accuracy changes from γ ∼ εh^d to γ ∼ ε.Source: §4.7.3 · p. 22

5 Numerical validation

Across three PDE examples, the proposed global-basis model achieves the lowest reported mean relative errors while using fewer parameters than DeepONet. The experiments also compare localized and FEM-based basis choices, with the localized version trading higher errors for a smaller model.

5 Numerical validation

  • The experiments compare proposed global-basis, exact FEM-basis, hybrid, and stacked DeepONet models, using 32 basis functions across the models.The proposed network has a coefficient network and a basis network; training uses Adam for 4,000 epochs.Source: §5 · pp. 22–23
  • Across the three PDEs, the proposed global-basis model has the lowest mean relative errors (0.08%, 0.28%, and 0.54%) and approximately 16% fewer parameters than DeepONet.Tables 1–3 cover the Allen–Cahn, KPP, and aggregation–diffusion examples.Source: §5 · p. 23, Table 1 · p. 23, Table 2 · p. 25
  • The localized version uses approximately 80–85% fewer parameters than the global version but has higher approximation errors and does not consistently outperform baselines in accuracy.It uses the fewest parameters in every experiment.Source: §5 · p. 23

5.1 Allen Cahn equation

  • Allen–Cahn maps initial conditions with independently sampled a1,a2 in [-0.5, 0.5] to solutions at T = 2 under homogeneous Dirichlet boundaries; Figure 3 shows two examples.The example initial conditions use tanh combinations of the first two sine modes.Source: §5.1 · p. 23, Figure 3 · p. 23

5.2 KPP example

  • KPP uses homogeneous Dirichlet boundaries and initial parameters sampled in [-1, 1]; its exact-FEM-basis model differs from FERN chiefly in the coefficient-network design.The proposed model uses a coefficient network derived from the numerical scheme, whereas FERN uses the standard DeepONet design.Source: §5.2 · pp. 24–25

5.3 Aggregation-diffusion equation example

  • The aggregation–diffusion example maps initial conditions parameterized by c0 in [1] [5] to terminal solutions on x in [-6, 6] and t in [0, 200].Figure 6 shows two predictions, and Table 3 reports numerical results.Source: §5.3 · p. 25, Figure 6 · p. 25, Table 3 · p. 25

6 Conclusion

The paper studies whether neural operators trained for forward PDE problems may approximate the numerical schemes that generate their data. It derives neural operators from numerical discretizations and reports an error analysis with a power scaling law observed in experiments.

6 Conclusion

  • For forward PDE solving, a trained neural operator may approximate the numerical scheme generating its training data.The paper frames this as a hypothesis about the relationship between operator learning and numerical methods.Source: §6 · pp. 25–26
  • The paper derives neural operators from numerical discretizations, yielding a network structure that returns to the classical Chen-Chen neural operator form.The derivation is based on the hypothesis that trained operators may mimic the schemes generating the data.Source: §6 · p. 26
  • The paper provides error analysis and derives a power scaling law that is widely observed in experiments.Source: §6 · p. 26

7 Appendix

The appendix derives the diffusion CFL scaling and stability bounds for the discretization, then gives constructions and complexity estimates for neural approximations of finite-element basis functions and locally supported functions.

7.1 CFL.

  • For fixed-degree Lagrange elements on shape-regular, quasi-uniform meshes, λmax(M^-1S) scales as h^-2, so choosing Δt ~ h^2 gives Nt ~ h^-2 for fixed T.An explicit reaction term may impose an additional restriction, but it is not a parabolic mesh restriction; the diffusion restriction controls the time-unrolled network depth.Source: §7.1 · p. 28

7.2 Proof of Lemma 3

  • Under the stated assumptions and diffusion CFL condition, the stability bounds include ∥B∥M→W = 1, ∥Ch∥W→M ≤ Δt, and ∥Lh∥M→M ≤ 1 + Δt.The heat update is non-expansive in the M-norm exactly when max_j |1 − Δtλj| ≤ 1.Source: §7.2 · p. 28, §7.1 · p. 28

7.3 Proof of the approximating to FEM basis

  • Each fixed-degree finite-element basis function has a ReLU approximation with L∞ error at most γ and size O(γ^-d(1 + log(1/γ))), uniformly in h and the basis index.The construction rescales the basis function to a reference cube, approximates it there, then composes with affine maps and clipping.Source: §7.3 · p. 30

7.4 Proof of the approximation with local support

  • A locally supported approximation combines a ReLU approximation with a mask and multiplier, vanishes outside D, and has error at most μ + ρ on Dκ.When μ ~ ρ ~ ε, the size of the function-approximation network dominates the mask and multiplier networks.Source: §7.4 · p. 31
Loading 2610.11120v1…