Numerical Operator Learning Methods

• Numerical operator learning is the approximation of mappings between infinite-dimensional function spaces, especially for PDE solution operators, using trainable surrogates.
• It employs discretization techniques, physics-informed constraints, and statistical methods to create efficient and reliable models for many-query tasks like optimization and uncertainty quantification.
• Key architectures include encoder–decoder frameworks, neural operators with nonlocal layers, and graph-based models that integrate numerical structure with data-driven learning.

Numerical operator learning is the approximation of mappings between function spaces—most prominently PDE solution operators—by trainable surrogates that are defined, trained, and evaluated through finite data and numerical discretizations. In one formulation it is a function-to-function regression problem for operators $G:U\to V$; in another, it is the direct learning of a discretized map from coefficients, initial data, boundary data, or design variables to discrete solution fields on grids, meshes, or sensor sets. Across these formulations, the central aim is to construct surrogates for many-query tasks such as parametric solution, optimization, and uncertainty quantification, while retaining enough numerical, physical, and statistical structure to remain faithful to the underlying operator (Subedi et al., 4 Apr 2025, Boullé et al., 2023, Kovachki et al., 2024).

1. Formal problem setting

At the most general level, operator learning is a regression problem in which the input and output are functions rather than finite-dimensional vectors. A typical target is an operator $G:V\to W$, with $V$ and $W$ separable Banach or Hilbert spaces of functions, learned from samples $(v_i,G(v_i))$ by empirical risk minimization over an operator class $\mathcal{F}$ (Subedi et al., 4 Apr 2025). In PDE applications this operator is often a solution map, such as coefficient-to-solution, forcing-to-solution, or initial-condition-to-solution.

Numerically, the infinite-dimensional problem is always coupled to a finite-dimensional representation. Functions are observed on grids, through basis coefficients, or via bounded linear observation operators such as $\phi:\mathcal{U}\to\mathbb{R}^n$ and $\varphi:\mathcal{V}\to\mathbb{R}^m$, frequently implemented as point-evaluation maps (Batlle et al., 2023, Mora et al., 2024). This is why the term “numerical operator learning” is used both for the numerical approximation of mappings between infinite-dimensional function spaces and for the learning of discretized PDE solution operators directly from numerical data.

A recurring distinction is between continuous operator learning and discrete operator learning. In the continuous view, one seeks an operator $G:U\to V$ that is meaningful independently of a specific mesh. In the discrete view, one learns a map such as $\mathcal{M}:a\mapsto u$ between tensors or coefficient vectors on a fixed grid or mesh, often with the same spatial dimensions in input and output (Zhu et al., 2023). This suggests that numerical operator learning is best understood as a spectrum: some methods target mesh-independent approximations of continuous operators, while others deliberately learn a discretization-specific numerical operator.

2. Core architectural families

Two architectural lineages organize much of the subject. The first is the encoder–decoder family. DeepONet uses a branch network on sampled input functions and a trunk network on output coordinates, yielding

$G:V\to W$0

which behaves like a learned low-rank factorization of the operator range (Boullé et al., 2023, Kovachki et al., 2024). PCA-Net has the same broad structure, but replaces learned encoders and decoders by PCA bases, so that the finite-dimensional neural network acts on reduced coordinates.

The second lineage is the neural-operator family, in which hidden layers are themselves nonlocal operators on feature fields. In a general form, layers apply

$G:V\to W$1

with $G:V\to W$2 an integral operator and $G:V\to W$3 a pointwise linear map (Kovachki et al., 2024). The Fourier Neural Operator specializes this to spectral convolution: Fourier transform, mode truncation, learned spectral multipliers, inverse Fourier transform, and pointwise mixing. In numerical terms, this is a nonlinear spectral method with global receptive field (Boullé et al., 2023, Enyeart et al., 2024).

Graph-based variants encode locality and multiscale structure more explicitly. Graph Neural Operators restrict interactions to local neighborhoods, while Multipole Graph Neural Operators are designed to reflect multi-scale and hierarchical low-rank structure in Green’s functions (Boullé et al., 2023). Multigrid-inspired architectures go further: MgNet makes the link between multigrid and CNNs explicit, and the enhanced V-cycle MgNet adds low-frequency correction to residuals because low-frequency errors decay slowly under smoothing iterations (Zhu et al., 2023).

A separate but related architectural current treats neural networks themselves as iterated operators. Fixed-point formulations of the form $G:V\to W$4 and relaxed Picard iterations are used to reinterpret diffusion models, AlphaFold recycling, transformers, neural integral equations, and graph networks as iterative operator-learning schemes, with convergence established under Lipschitz or contraction conditions (Zappala et al., 2023). This places numerical operator learning in direct dialogue with classical iterative methods.

3. Discretization-aware and numerically structured models

A defining feature of numerical operator learning is that many successful models are not merely expressive approximators; they are structured by the numerical representation of the PDE. EV-MgNet is emblematic: it begins from a multigrid-based convolutional neural network architecture known as MgNet and introduces a low-frequency correction structure for residuals. The resulting model captures low-frequency features considerably better than the standard V-cycle MgNet and is more robust in case of low- and high-resolution data during training and testing, respectively (Zhu et al., 2023). Its design explicitly mirrors multigrid frequency separation.

DimOL introduces a different kind of structure, drawn from dimensional analysis. Its ProdLayer can be seamlessly integrated into FNO-based and Transformer-based PDE solvers and is designed to capture sum-of-products structures inherent in many physical systems (Song et al., 2024). Empirically, DimOL models achieve up to 48% performance gain within the PDE datasets, and the analysis of Fourier components’ weights can symbolically discern the physical significance of each term (Song et al., 2024). The broader implication is that inductive bias need not be restricted to locality or spectral sparsity; it can also encode algebraic structure of PDE terms.

Finite Operator Learning makes the discretization dependence fully explicit. It uses an uncomplicated feed-forward neural network to directly map the discrete design space to the discrete solution space at a finite number of sensor points in an arbitrary shape domain, with the discretized governing equations, as well as the design and solution spaces, derived from well-established numerical techniques (Rezaei et al., 2024). In the reported realization, the Finite Element Method is used to approximate fields and their spatial derivatives, and the target problem is the steady-state heat equation within heterogeneous materials that exhibits significant phase contrast and possibly temperature-dependent conductivity (Rezaei et al., 2024). Rather than attempting mesh independence, this approach treats the learned operator as a numerical surrogate for a specific discretization.

4. Physics-informed and structure-preserving training

Physics-informed operator learning augments data fitting with explicit numerical or physical constraints. In statistical form, it can be written as regularized empirical risk minimization with additional penalties for PDE residuals, boundary conditions, conservation laws, or variational identities (Subedi et al., 4 Apr 2025). This unifies pure data-driven regression with physics-based regularization.

Finite Operator Learning is explicitly data-free. It can parametrically solve partial differential equations in a data-free manner and provide accurate sensitivities, meaning the derivatives of the solution space with respect to the design space (Rezaei et al., 2024). Its Sobolev training minimizes a multi-objective loss function that includes the discretized weak form of the energy functional, boundary conditions violations, and the stationarity of the residuals with respect to the design variables (Rezaei et al., 2024). Because the network’s tangent matrix is directly used for gradient-based optimization, the learned operator is tied not only to forward prediction but also to sensitivities and design updates.

PENCO extends the same philosophy to stiff time-dependent phase-field models. It introduces a hybrid operator-learning framework that integrates physical laws and numerical structure within a data-driven architecture. Its objective combines an enhanced $G:V\to W$5 Gauss-Lobatto collocation residual around the temporal midpoint, a Fourier-space numerical consistency term that captures the balanced behavior of semi-implicit discretizations, an energy-dissipation constraint that ensures thermodynamic consistency, and additional low-frequency spectral anchoring and teacher-consistency mechanisms (Bamdad et al., 4 Dec 2025). The resulting operator is not simply trained to match trajectories; it is constrained to behave like a stable, structure-preserving time integrator.

This numerics-aware view differs from generic physics-informed residual minimization. The aim is not only to enforce the PDE pointwise, but to align the learned operator with a chosen discretization, variational principle, or Lyapunov structure. This suggests a convergence of numerical analysis and machine learning in which discretizations, residuals, and energy laws become architectural and loss-design primitives.

5. Statistical, kernel, Gaussian-process, and noisy-data perspectives

A parallel literature formulates numerical operator learning in explicitly statistical terms. Operator learning becomes function-to-function regression with approximation, statistical, truncation, and discretization error all contributing to the total risk (Subedi et al., 4 Apr 2025). This viewpoint is important because many numerical claims—generalization across meshes, sample efficiency, or uncertainty quantification—are ultimately statistical claims about finite data.

Kernel methods provide one of the clearest non-neural formulations. A general framework approximates $G:V\to W$6 by

$G:V\to W$7

where $G:V\to W$8 and $G:V\to W$9 are linear observation operators, $V$0 and $V$1 are optimal recovery maps, and $V$2 is a kernel approximation on the finite-dimensional measurement space (Batlle et al., 2023). Even when using vanilla kernels such as linear, Matérn, or rational quadratic, this framework is competitive in terms of cost-accuracy trade-off and offers simplicity, interpretability, convergence guarantees, a priori error estimates, and Bayesian uncertainty quantification (Batlle et al., 2023).

Gaussian-process operator learning pushes this further by placing a GP prior on the bilinear form $V$3, so that $V$4. Numerically, the GP is defined on $V$5, can use separable product kernels and Kronecker product matrix representations, and can take a neural operator as its mean function (Mora et al., 2024). This yields a hybrid GP/NN framework in which the GP interpolates residual structure around a parametric mean, and it enables zero-shot data-driven models for accurate predictions without prior training (Mora et al., 2024).

Noise in both variables changes the statistical picture substantially. Error-in-variables modelling for operator learning shows that ordinary least squares has an operator analogue of attenuation bias: in the nonlinear setting there is underprediction of the action of the Burgers operator in the presence of noise in the independent variable (Patel et al., 2022). EiV formulations for MOR-Physics and DeepONet reduce this bias and robustly recover operators in high-noise regimes that defeat OLS operator learning (Patel et al., 2022). This is especially relevant because both independent and dependent variables in operator learning are signals and are therefore naturally subject to measurement error.

At a more fundamental level, data complexity can itself be intractable. Lower bounds on $V$6-widths for general classes of Lipschitz and Fréchet differentiable operators demonstrate a curse of data-complexity: learning on such general classes requires a sample size exponential in the inverse of the desired accuracy $V$7 (Kovachki et al., 2024). The same paper shows, using FNO as a case study, that parametric efficiency implies data efficiency on narrower operator classes (Kovachki et al., 2024). This identifies a central principle: data efficiency in numerical operator learning is inseparable from structural assumptions on the target operator class.

6. Empirical landscape and application domains

Quantitative benchmark suites have made numerical operator learning unusually comparative. In EV-MgNet, representative relative $V$8 errors at the highest reported resolutions include approximately $V$9 on 1D Burgers at $W$0, $W$1 on 1D KdV at $W$2, $W$3 on 2D Darcy at $W$4, and $W$5 on 2D Navier–Stokes with $W$6, $W$7, and $W$8 (Zhu et al., 2023). The same study shows that the enhanced model can be trained at low resolution and tested at much higher resolution with small degradation in error (Zhu et al., 2023).

Phase-field dynamics furnish a demanding long-horizon testbed. PENCO reports final-time errors around 0.02 on Allen–Cahn, 0.04 on Cahn–Hilliard, 0.05 on Swift–Hohenberg, 0.04 on Phase Field Crystal, and 0.03 on Molecular Beam Epitaxy when $W$9, while the corresponding FNO-4D and MHNO baselines are substantially larger in the hardest regimes (Bamdad et al., 4 Dec 2025). On Molecular Beam Epitaxy, FNO-4D is approximately 0.80, MHNO approximately 0.73, PENCO approximately 0.03, and pure physics approximately 0.15 (Bamdad et al., 4 Dec 2025). These results are tied to long-term stability, suppression of temporal error accumulation, and preservation of physically consistent evolution.

Kernel and GP methods are also empirically competitive on standard operator-learning suites. On seven PDE benchmarks, a best Matérn or rational quadratic kernel achieved 2.15% on Burgers, 2.75% on Darcy, $(v_i,G(v_i))$0% on Advection I, 11.44% on Advection II, 1.00% on Helmholtz, 5.18% on structural mechanics, and 0.12% on Navier–Stokes, matching or beating DeepONet, FNO, PCA-Net, or PARA-Net on a majority of tasks (Batlle et al., 2023). In the GP/NN hybrid setting, FNO-mean GP reduced Burgers from 1.93% to 0.08%, Darcy from 2.41% to 2.19%, Advection from 0.66% to 0.23%, and structural mechanics from 6.62% to 6.49%, while zero-shot zero-mean GP remained competitive with far fewer parameters (Mora et al., 2024).

Application domains now span steady and unsteady PDEs, elliptic and parabolic operators, coefficient-to-solution maps, sequence-to-sequence fluid prediction, and PDE-constrained optimization. A notable optimization-oriented example is Finite Operator Learning, where the network’s tangent matrix is used directly for gradient-based optimization to improve the microstructure’s heat transfer characteristics (Rezaei et al., 2024). This suggests that forward surrogate accuracy is only one use case; differentiable operator surrogates are increasingly valued for sensitivities and design.

7. Theory, practice, and open directions

Several clarifications follow from the current literature. Numerical operator learning is not synonymous with mesh-independence: some methods emphasize evaluation on different meshes or point sets, while others are explicitly discretization dependent, as in Finite Operator Learning (Boullé et al., 2023, Rezaei et al., 2024). Nor is numerical operator learning exclusively neural: kernel methods, Gaussian processes, and random feature models remain central because they offer interpretability, uncertainty quantification, and convergence guarantees (Batlle et al., 2023, Mora et al., 2024).

The sharp difficulty boundary is still unsettled. For general Lipschitz or Fréchet differentiable operator classes, the existing theory points to exponential dependence on inverse accuracy, described as a curse of parametric complexity and a curse of data-complexity (Kovachki et al., 2024, Kovachki et al., 2024). More favorable algebraic rates appear only on narrower classes, such as linear operators, holomorphic operators, and certain operator families efficiently approximated by FNO. A plausible implication is that successful numerical operator learning will continue to depend less on universal approximation in the abstract and more on discovering operator classes with exploitable structure.

Another active strand concerns iteration and stability. Fixed-point formulations $(v_i,G(v_i))$1, Banach–Caccioppoli convergence, and relaxed Picard schemes have been used to analyze iterative graph neural networks, transformers, neural integral equations, diffusion models, and AlphaFold recycling (Zappala et al., 2023). This suggests that future neural operators may increasingly be designed as convergent iterative solvers rather than single-pass regressors.

Training practice is also beginning to stabilize. Across DeepONets, Fourier Neural Operators, and Koopman autoencoders, GELU was best in every single experiment, any positive dropout rate degraded performance, stochastic weight averaging improved accuracy when its learning rate was equal to or one-tenth the base learning rate, and learning-rate finders were inconsistent (Enyeart et al., 2024). These are implementation-level observations, but in numerical operator learning they matter because hyperparameter choices interact directly with approximation quality and cross-resolution robustness.

Recurring open problems include active data collection, rigorous uncertainty quantification frameworks, theory for multiresolution generalization, frequency-decay analysis in multigrid-informed architectures, non-periodic boundary conditions and complex geometries, multi-physics, time-adaptive operators, and standardized, realistic benchmarks with centralized repositories of datasets and models (Subedi et al., 4 Apr 2025, Zhu et al., 2023, Bamdad et al., 4 Dec 2025). Together these directions indicate that numerical operator learning is evolving from a collection of surrogate architectures into a broader synthesis of numerical analysis, statistical learning, and scientific computing.