Solution-Operator Learning Methods

• Solution-operator learning methods are techniques that approximate mappings between infinite-dimensional function spaces using neural and statistical architectures to efficiently solve parameterized PDEs.
• They leverage designs like DeepONet, Fourier neural operators, and RKHS-based methods to produce mesh-independent surrogates and enable uncertainty quantification.
• Practical implementations blend empirical loss minimization with physics-informed regularization to ensure stability and convergence in complex PDE simulations.

Solution-operator learning methods seek to approximate maps between infinite-dimensional function spaces—typically the solution mapping of parameterized partial differential equations (PDEs)—by employing neural or statistical architectures that capture operator-theoretic, functional-analytic, and physical properties. These methods have become central in computational mathematics, scientific machine learning, optimal control, and uncertainty quantification, enabling efficient surrogates for PDE solvers, inference in inverse problems, and simulation-driven modeling far beyond the reach of traditional grid-based numerics.

1. Mathematical Formulation and Operator Learning Paradigms

Let $D \subset \mathbb{R}^d$ be a spatial domain, and consider a parameterized PDE

$$\mathcal L_\alpha u = f \quad \text{in } D, \qquad u = g \quad \text{on } \partial D,$$

where $\alpha, g$ are fixed and $f$ varies over a Banach or Hilbert function space $\mathcal{F}$. The solution operator $\mathcal{G}^\dagger : \mathcal{F} \to \mathcal{U}$ maps $f$ to the unique solution $u$. In abstract terms, solution-operator learning aims to learn $\mathcal{G}^\dagger$ (or more generally, an operator $G: V \to W$ between Banach or Hilbert spaces of functions) from input-output samples or, in some frameworks, from the PDE structure itself (Boullé et al., 2023, Subedi et al., 4 Apr 2025).

The goal is to construct a parametric family $G_\theta$ (neural, kernel, or polynomial-based) minimizing empirical or physics-informed risk,

$$L(\theta) = \frac1N \sum_{i=1}^N \|G_\theta(f_i) - u_i\|^2_W + \text{(regularization)},$$

where $(f_i, u_i)$ are (possibly noisy) samples, or by unsupervised minimization of action, residual, or energy functionals.

2. Representative Architectures and Parameterizations

A variety of architectures implement solution-operator learning, differing in inductive bias, computational mechanism, and theoretical guarantees:

• Function-valued RKHS neural operator: Bao et al. (Bao et al., 2022) embed the operator learning problem in a function-valued reproducing kernel Hilbert space, leveraging operator-valued kernels $K$ and neural network parameterizations for Hilbert–Schmidt integral kernels. The learned operator takes the form

$$\hat G_n = \sum_{i=1}^n K(f_i, \cdot) A_i, \quad A_i \in \mathcal{U}.$$

Neural networks encode both the $f_i$ similarity kernel and the spatial integral kernel, enabling data-efficient inference and mesh-independence.

• DeepONet and trunk-branch decompositions: DeepONet (Boullé et al., 2023) uses a branch net encoding input function samples and a trunk net encoding output coordinates, reconstructing the operator as

$$\mathcal G_\theta(f)(y) = \sum_{k=1}^p b_k(f(x_1), \dots, f(x_m))\, t_k(y).$$

• Fourier/pseudo-differential neural operators (FNO, PDNO): FNO (Boullé et al., 2023), and generalizations such as PDNO (Shin et al., 2022), parameterize convolution or pseudo-differential kernels in frequency (and optionally position), implementable as

$$u_{\ell+1}(x) = \sigma\big( W_\ell u_\ell(x) + \mathcal{F}^{-1} [a_\theta(x,\xi) \cdot \widehat{u}_\ell(\xi)](x) \big),$$

systematically enforcing continuity in Sobolev/Hörmander classes and supporting space- or time-dependent coefficients.

• Graph neural operators (GNO, MGNO): For PDEs whose Green's functions have local or hierarchical structure, GNOs (Boullé et al., 2023) perform message passing on mesh or graph representations, efficiently capturing nonlocal interactions in complex geometries.
• Polynomial Chaos and weighted least squares: PCE methods (Sharma et al., 28 Aug 2025) model the unknown solution as an expansion in polynomial chaos bases over stochastic input variables, reducing operator learning to solving explicit regression or constrained systems for expansion coefficients with built-in uncertainty quantification.
• Energy-based and variational learning: Physics-informed MLPs or FNOs trained by minimizing elementwise or global discrete energy (Ritz or Galerkin forms) avoid reliance on solved data (Larson et al., 2024, Xu et al., 2023), supporting label-free operator learning and matrix-free end-to-end training.

3. Regularization, Losses, and Training Strategies

Operator learning methods balance data fidelity with stability and generalization via:

• Empirical $L^2$ or relative-error loss over function evaluations or $L^2$-normed outputs.
• RKHS norm regularization, enforcing boundedness and controlling estimation error; in practice absorbed into neural weight decay (Bao et al., 2022).
• Sobolev (derivative) supervision: Combined $L^2$ and derivative mismatch loss, approximated with moving-least-squares on unstructured meshes, yields 10–30% lower error, noise robustness, and faster local convergence (Cho et al., 2024).
• Physics-informed residuals: PDE and boundary residuals are penalized at collocation points, enabling data-free training (Larson et al., 2024, Bi et al., 2024), often combined with weak-form projections or energy minimization.
• Variance-based or Christoffel weighting: In optimal weighted least squares frameworks, sample weights and measures are adapted to condition the operator-level Gram matrix, yielding $O(N \log N)$ sample complexity and geometric regularization (Turnage et al., 11 Dec 2025).

Supervised, label-free, and unsupervised regimes appear depending on application and data availability; physics-based architectures can entirely avoid ground-truth solution computation.

4. Theoretical Guarantees, Convergence and Complexity

Key theoretical properties and results include:

• Universal Approximation and RKHS Representer Theorems: Neural operator models such as DeepONet and FNO are universal in the space of continuous operators; RKHS approaches guarantee that minimizers have explicit kernel expansion forms (Bao et al., 2022, Boullé et al., 2023, Subedi et al., 4 Apr 2025).
• Curse of Parametric Complexity: For generic $C^r$- or Lipschitz-regular operators, network size and required discretization grow exponentially in target error, unless problem structure is exploited. This lower bound holds for FNO, DeepONet, and linear architectures (Lanthaler et al., 2023).
• Structure-informed architectures: For Hamilton–Jacobi equations, HJ-Net explicitly encodes characteristic flows and beats the generic curse, achieving polynomial complexity in error (Lanthaler et al., 2023).
• Sample complexity and stability: Operator-level Christoffel function weighting enables $M \sim N\log N$ sample scaling for $N$-dimensional approximation, with uniformly conditioned regression and nonasymptotic stability bounds (Turnage et al., 11 Dec 2025).
• Continuity and regularity: PDNO provides mathematical control in Sobolev spaces via symbol class constraints; time-modulated FNO achieves global Lipschitz (and hence stable) operator learning (Shin et al., 2022, Park et al., 2023).

5. Numerical Benchmarks and Practical Guidelines

Empirical evaluation spans canonical elliptic, parabolic, and hyperbolic PDEs under varying data and mesh regimes:

Method	PDEs/Domains (sample)	Relative $L^2$ Error (typical)	Key Features
RKHS Neural Operator (Bao et al., 2022)	Advection, Burgers, KdV, Poisson	1.5–4.8% (low data)	Mesh-independence, up-sampling, RKHS control
Sobolev Training (Cho et al., 2024)	Darcy2d, NS2d, Heat, Elasticity	10–30% error reduction	Derivative matching, noise-/grid-robust
DeepONet, FNO (Boullé et al., 2023)	1D Burgers, 2D Darcy	1–5%	Efficient with 100–1000 PDE solves
PCE/PC	1D/2D Advection, Burgers, Heat	$10^{-9}$–$10^{-4}$	Closed-form training, UQ, no neural nets
PDNO (Shin et al., 2022)	Darcy, Navier–Stokes	$1.4 \times 10^{-3}$	PDO theory, smooth symbols, Sobolev estimates
Energy-MLP (Larson et al., 2024)	Poisson, nonlinear elasticity	Theor. FEM bounds	Data-free, discrete energy minimization
MeshONet (Xiao et al., 21 Jan 2025)	Mesh generation	$<1\%$ geom. error	Dual-branch, multi-input, 4–5 orders faster
One-shot local operator (Jiao et al., 2021)	1D/2D linear, nonlinear PDEs	$1$–$5\%$	Only one global PDE solve, locality principle

MeshONet demonstrates high efficiency in mesh generation tasks, generalized to variable geometries. CHONKNORIS achieves machine-precision solution operator learning for forward and inverse nonlinear PDEs by regressing the Cholesky factors of Tikhonov-regularized Newton–Kantorovich updates, with convergence guarantees in terms of surrogate accuracy (Bacho et al., 25 Nov 2025).

6. Recent Extensions and Open Directions

Solution-operator learning is a rapidly active area, with open directions including:

• Unsupervised and label-free learning: Trajectory-sampling and amortized-variational frameworks for mean-field games enable mesh-free, unsupervised, and dimension-agnostic operator learning (Huang et al., 2024).
• Adaptive data and active sampling: Statistical theory for linear operators shows that active, data-adaptive sampling can achieve super-parametric convergence rates, far surpassing the classical $n^{-1/2}$ (Subedi et al., 4 Apr 2025).
• Multivariable and manifold-valued operators: Modern architectures can encode multiple input/output fields (MeshONet dual-branch), variable geometries, and learning on manifolds.
• Operator uncertainty quantification (OUQ): PCE and Bayesian/conformal prediction methods provide exact or distribution-free UQ for operator predictions (Sharma et al., 28 Aug 2025, Subedi et al., 4 Apr 2025).
• Scalability and foundation models: FONKNORIS aggregates multiple expert surrogates for cross-PDE generalization; scaling laws and standardized PDE benchmarks are under active investigation (Bacho et al., 25 Nov 2025, Subedi et al., 4 Apr 2025).

7. Theoretical Limits, Practical Caveats, and Recommendations

While operator learning delivers dramatic practical gains for PDE surrogate modeling and inverse problems, inherent statistical, approximation, and computational limitations remain. Universality and stability hinge on function-space properties, architecture design, and data distribution. Physics-informed regularization, derivative supervision, and architecture–PDE matching remain critical for avoiding overfitting, underfitting, and instability.

For practitioners, best practices include adaptation of mesh, basis, and kernel to PDE structure; use of Sobolev-based losses when regularity allows; kernel or Christoffel weighting for efficient sampling; and leveraging physics-informed losses or variational forms for data-scarce or unsupervised scenarios.

---

Comprehensive references and further details can be found in (Bao et al., 2022, Cho et al., 2024, Xu et al., 2023, Subedi et al., 4 Apr 2025, Boullé et al., 2023, Sharma et al., 28 Aug 2025, Turnage et al., 11 Dec 2025, Shin et al., 2022, Xiao et al., 21 Jan 2025, Bacho et al., 25 Nov 2025, Larson et al., 2024, Jiao et al., 2021, Benth et al., 2024, Lanthaler et al., 2023, Li et al., 2022, Huang et al., 2024).