Property-Preserving Kernel Operator Learning

• Property-preserving kernel-based operator learning is a framework that uses kernel and RKHS methods to embed critical physical or analytical constraints into mappings between function spaces.
• The methodology incorporates operator-valued kernels, divergence-free bases, and structured recovery maps to maintain PDE operators, boundary conditions, and other invariants directly in the learning process.
• Recent studies show that these approaches yield highly accurate surrogates with explicit error bounds and interpretability, enhancing applications in PDE solving, fluid dynamics, and dynamical systems.

Property-preserving kernel-based operator learning denotes a class of kernel and RKHS methods that learn maps between function spaces, Hilbert spaces, or Banach spaces while embedding structural constraints into the hypothesis space, the kernel, the basis, or the learning objective. In the recent arXiv literature, the preserved structure may be the governing PDE operator and boundary conditions, the measurement and recovery geometry of partial observations, analytic incompressibility and periodicity, temporal differentiability and Koopman spectral structure, or dissipativity certificates for unknown nonlinear systems (Hu et al., 10 May 2026, Batlle et al., 2023, Sharma et al., 17 Feb 2026, Withanachchi, 23 Aug 2025, Ye et al., 31 Oct 2025). The unifying premise is that the learned operator should respect salient properties of the target problem by construction or through the variational formulation, rather than recover them only approximately from end-to-end data fitting.

1. Scope and meaning of property preservation

Across the literature, “property-preserving” is not a single formal definition but a family of design principles. In PDE operator learning, the preserved object is the differential structure itself: the PDE operator $P_s$ appears directly in the loss, boundary data are incorporated into the formulation, and the resulting estimator defines a solution operator $T:h\mapsto u$ instead of fitting a single solution field (Hu et al., 10 May 2026). In partial-observation operator learning, the emphasis is different: the approximation is built as

$\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$

so that measurement maps $\phi,\varphi$ and optimal recovery maps $\psi,\chi$ preserve the observation model and the reconstruction geometry (Batlle et al., 2023).

A distinct form appears in incompressible-flow surrogates, where property preservation is transferred into the representation itself. Output velocity fields are expanded in a matrix-valued kernel basis that is divergence-free, periodic when required, and multiscale when turbulence-related scaling is to be encoded. Because prediction proceeds only through this basis, incompressibility, periodicity, and turbulence-related scaling are inherited analytically rather than approximately (Sharma et al., 17 Feb 2026). In spatio-temporal dynamics, operator-valued kernels and time-regularized RKHS norms preserve vector-valued spatial structure, temporal differentiability, Sobolev regularity, and spectral characteristics of kernel Koopman approximations (Withanachchi, 23 Aug 2025). In dissipativity learning, the preserved object is not a state-to-state or input-to-output map, but a storage/supply-rate certificate represented by Hilbert–Schmidt operators in RKHS feature space (Ye et al., 31 Oct 2025).

This variation matters conceptually. A common misconception is to identify property preservation exclusively with hard physics constraints. The partial-observation framework shows that preservation can also be functional-analytic: exact consistency with measurements, mesh invariance, and structured reconstruction can be the primary invariant even when no PDE is imposed explicitly (Batlle et al., 2023). Conversely, the PDE-discovery framework shows that some kernel operator learners preserve stability, differential consistency, and boundary residual structure without guaranteeing conservation laws, symmetry, positivity, monotonicity, or the maximum principle (Long et al., 2022).

2. Functional-analytic architecture and RKHS formulations

The common substrate of these methods is reproducing kernel Hilbert space theory. In scalar-valued learning, regularized empirical risk minimization with square loss yields kernel ridge regression. In vector-valued or operator learning, the output space becomes a Hilbert space $\mathcal Y$, the kernel becomes operator-valued,

$$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$

and the representer form becomes

$$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$

This places learned maps inside a vvRKHS whose geometry is determined by the choice of $\Gamma$, for example linear, separable, or more general operator-valued kernels (Rosasco, 22 Sep 2025).

The partial-observation operator-learning framework makes this structure explicit by separating observation, reduced learning, and recovery. The target operator is

$\mathcal G^\dagger:\mathcal U\to \mathcal V,$

but only partial linear measurements $T:h\mapsto u$0 and $T:h\mapsto u$1 are observed. Optimal recovery maps

$T:h\mapsto u$2

reconstruct minimum-norm elements consistent with the measurements, and the reduced operator

$T:h\mapsto u$3

is then learned in a vector-valued RKHS. The structured approximation

$T:h\mapsto u$4

belongs exactly to the operator-valued RKHS

$T:h\mapsto u$5

which is the cleanest statement that the learned operator is confined to a function class already compatible with the observation and recovery pipeline (Batlle et al., 2023).

This architecture is one source of the field’s interpretability. The kernel is not only a similarity function; it specifies regularity, admissible coupling between outputs, and, in many formulations, the admissible operator class itself. A plausible implication is that “property preservation” in kernel operator learning is often best understood as a statement about the admissible hypothesis space before optimization begins.

3. Differential operators and solution-operator learning for PDEs

A central recent development is the direct kernel learning of PDE solution operators. For a linear inhomogeneous boundary-value problem,

$T:h\mapsto u$6

with

$T:h\mapsto u$7

the Dirichlet case is rewritten as

$T:h\mapsto u$8

where $T:h\mapsto u$9 bundles the interior forcing and boundary data. The learning goal is then the solution operator

$\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$0

so PDE solving is reinterpreted as learning a map between function spaces (Hu et al., 10 May 2026).

The estimator is defined by regularized empirical risk minimization in an RKHS $\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$1: $\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$2 with

$\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$3

A key result is the closed-form operator expression

$\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$4

where $\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$5 and $\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$6. The corresponding differential representer theorem gives

$\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$7

with

$\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$8

Here the generalized Gram matrix $\bar{\mathcal G}=\chi\circ \bar f\circ \phi,$9 is a PDE-aware analogue of the usual kernel matrix, because the differential operator acts on the kernel itself (Hu et al., 10 May 2026).

The property-preserving mechanism is the differential reproducing property. If $\phi,\varphi$0, then $\phi,\varphi$1 is a bounded linear operator from $\phi,\varphi$2 to $\phi,\varphi$3, and

$\phi,\varphi$4

Thus the governing operator is not merely penalized; it is built into the RKHS geometry and into the basis functions used by the estimator. Boundary conditions are incorporated through the formulation, and the paper explicitly characterizes the method as shifting from a PDE solver to a solution-operator solver (Hu et al., 10 May 2026).

The theory separates total error into approximation and estimation components,

$\phi,\varphi$5

Under the source condition

$\phi,\varphi$6

the approximation error satisfies

$\phi,\varphi$7

With $\phi,\varphi$8, the paper proves a uniform convergence rate of the form

$\phi,\varphi$9

up to constants and logarithmic factors (Hu et al., 10 May 2026).

4. PDE discovery and surrogate operator construction

A different kernel route to operator learning begins not from a known PDE operator, but from its reconstruction. The three-step framework for PDE discovery and operator learning first denoises noisy solution data by kernel smoothing, then learns the PDE operator by kernel regression on solution-derived features, and finally solves the learned surrogate PDE for new source or boundary data (Long et al., 2022).

In the smoothing stage, each noisy training solution is replaced by the RKHS estimator

$\psi,\chi$0

and derivatives are obtained analytically by differentiating the kernel expansion,

$\psi,\chi$1

The learned PDE operator then has representer form

$\psi,\chi$2

where $\psi,\chi$3 collects local features such as $\psi,\chi$4, $\psi,\chi$5, and derivatives of $\psi,\chi$6 (Long et al., 2022).

For a new forcing term, the method solves a regularized residual problem over an RKHS,

$\psi,\chi$7

This framework preserves several forms of structure: stability through RKHS norm penalties, consistency with differential operators through kernel differentiation, and boundary conditions through boundary residual terms. At the same time, it explicitly does not guarantee by construction symmetry, positivity, conservation laws, monotonicity, the maximum principle, or exact well-posedness of the learned PDE (Long et al., 2022).

This contrast with direct PDE-operator learning is instructive. When the differential operator is known, kernels can embed it directly in the loss and basis. When the operator must be discovered, preservation shifts from exact physics encoding to regularized surrogate consistency. A plausible implication is that the degree of preservation depends strongly on whether the structural prior is known a priori or has to be inferred from data.

5. Analytic preservation of incompressibility, periodicity, and turbulence

For incompressible Navier–Stokes surrogates, property preservation is implemented at the level of the output basis rather than through residual penalties. The learning problem is an operator map

$\psi,\chi$8

from input functions such as initial conditions, boundary data, or geometry parameters to velocity fields. Instead of predicting the velocity field directly, the method first expresses each output in a property-preserving kernel interpolant,

$\psi,\chi$9

and then learns an operator-valued kernel interpolant from input samples to the coefficients $\mathcal Y$0. The final prediction is reconstructed only from this basis, so the output inherits the basis properties analytically (Sharma et al., 17 Feb 2026).

The central construction is the divergence-free matrix-valued kernel generated from a scalar positive-definite kernel $\mathcal Y$1: $\mathcal Y$2 For radial kernels this is equivalent to

$\mathcal Y$3

Each column of $\mathcal Y$4 is divergence-free in the evaluation coordinate, so any interpolant of the form $\mathcal Y$5 satisfies

$\mathcal Y$6

pointwise everywhere, independent of the coefficients or sample locations (Sharma et al., 17 Feb 2026).

Periodicity is imposed by embedding $\mathcal Y$7 into $\mathcal Y$8 via

$\mathcal Y$9

forming a periodic scalar kernel $$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$0, and then defining the periodic divergence-free kernel

$$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$1

To encode turbulence-related scaling, the paper uses a multiscale additive kernel

$$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$2

with

$$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$3

A combined kernel $$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$4 simultaneously enforces incompressibility, periodicity, and turbulence scaling (Sharma et al., 17 Feb 2026).

The operator stage uses a diagonal operator-valued kernel and yields the two-stage pipeline

$$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$5

The paper emphasizes that physics preservation is guaranteed by the basis, not by optimization. It also reports that the overall pipeline uses only two trainable scalar parameters: one kernel shape parameter for the property-preserving interpolant and one for the operator kernel (Sharma et al., 17 Feb 2026).

6. Spatio-temporal dynamics, Koopman structure, and dissipativity certificates

Kernel operator learning has also been extended to time-dependent vector fields and dynamical systems. In the OV-RKHS and kernel Koopman framework, the target is a time-dependent vector field $$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$6 on $$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$7, learned in an operator-valued RKHS with kernel

$$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$8

A central choice is the separable kernel

$$\Gamma:\mathcal X\times \mathcal X\to \mathcal L(\mathcal Y),$$9

together with the time-regularized norm

$$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$0

The corresponding representer theorem yields a finite expansion involving $$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$1, and the paper proves Sobolev-type approximation bounds for both the field and its time derivative as well as operator-norm and spectral convergence for kernel Koopman approximations (Withanachchi, 23 Aug 2025).

The Koopman perspective supplies the operator-theoretic interpretation. Nonlinear state evolution is recast as linear evolution of observables, and a reduced spectral forecast is written as

$$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$2

The expository treatment of kernel learning for functions, operators, and dynamical systems formulates this as a regularized operator-learning problem in an RKHS of observables, where the empirical objective is

$$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$3

The minimizer approximates the Koopman operator restricted to the chosen RKHS, and self-adjointness appears when the underlying Markov process satisfies detailed balance (Rosasco, 22 Sep 2025).

A further extension replaces state prediction by certificate learning. For the discrete-time nonlinear system

$$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$4

dissipativity is defined by the existence of a storage function $$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$5 and supply rate $$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$6 such that

$$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$7

The RKHS formulation lifts variables into feature spaces and represents

$$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$8

where $$f_\lambda(x)=\sum_{i=1}^n \Gamma(x,x_i)c_i,\qquad c=(\Gamma+n\lambda I)^{-1}y.$$9 and $\Gamma$0 are Hilbert–Schmidt operators. Learning is posed as a one-class SVM-like margin maximization problem over $\Gamma$1, and a representer theorem reduces the infinite-dimensional problem to a finite-dimensional convex program involving kernel Gram matrices (Ye et al., 31 Oct 2025).

This development broadens the meaning of operator learning. The learned object need not be a direct surrogate for dynamics; it may instead be a provable energy-balance certificate. A plausible implication is that property preservation in kernel methods can target analysis and verification objectives, not only forward prediction.

7. Guarantees, empirical performance, and limitations

A recurring feature of kernel operator learning is the availability of explicit theory. The partial-observation framework proves a priori error bounds, asymptotic convergence under dense sampling, and identifies the structured operator-valued RKHS in which the learned map lives. It also interprets the regularized estimator as the posterior mean of an operator-valued Gaussian process, with predictive covariance and deterministic error bounds derived from the same kernel geometry (Batlle et al., 2023). The PDE solution-operator framework supplies approximation and estimation error decompositions and uniform convergence rates under source conditions and suitable regularization (Hu et al., 10 May 2026). The spatio-temporal OV-RKHS framework proves convergence in Sobolev norms for both functions and time derivatives, operator-norm convergence of Koopman approximations, and convergence of eigenvalues and eigenfunctions (Withanachchi, 23 Aug 2025). The dissipativity framework supplements representer-theorem reduction with confidence bounds on the dissipation rate and $\Gamma$2-gain (Ye et al., 31 Oct 2025).

Empirically, kernel methods have been reported as competitive or superior on a broad set of operator-learning benchmarks. In the general partial-observation framework, benchmarks include Burgers’ equation, Darcy flow, Advection I and II, Helmholtz, structural mechanics, and Navier–Stokes. Reported examples include Darcy with kernel RQ about $\Gamma$3, Helmholtz with the best kernel about $\Gamma$4, Navier–Stokes with the best kernel about $\Gamma$5, and Advection I with the linear kernel essentially exact, while structural mechanics is identified as the hardest case for the vanilla kernel method (Batlle et al., 2023). In direct PDE solution-operator learning, Darcy and Helmholtz experiments report relative $\Gamma$6 errors around $\Gamma$7, runtime below $\Gamma$8 seconds on Darcy cases, and seconds versus about $\Gamma$9 hours against Green-operator learning in high-frequency Helmholtz regimes (Hu et al., 10 May 2026).

The incompressible-flow results are especially strong where the preserved property is analytically hard to satisfy with neural surrogates. The property-preserving kernel method is reported to achieve up to six orders of magnitude lower relative $\mathcal G^\dagger:\mathcal U\to \mathcal V,$0 errors upon generalization and to train up to five orders of magnitude faster compared to neural operators, while enforcing incompressibility analytically and yielding exactly zero divergence. Reported benchmark highlights include errors as low as $\mathcal G^\dagger:\mathcal U\to \mathcal V,$1 on 2D flow past a cylinder, $\mathcal G^\dagger:\mathcal U\to \mathcal V,$2 on the 2D lid-driven cavity, around $\mathcal G^\dagger:\mathcal U\to \mathcal V,$3 on 2D Taylor–Green vortices, and near $\mathcal G^\dagger:\mathcal U\to \mathcal V,$4 on the 2D backward-facing step (Sharma et al., 17 Feb 2026).

These advantages are accompanied by explicit caveats. Full kernel methods can be computationally heavy because dense linear systems scale cubically with the number of sample points. For the PDE solution-operator framework, an online low-rank kernel regression scheme is introduced to reduce memory from $\mathcal G^\dagger:\mathcal U\to \mathcal V,$5 to $\mathcal G^\dagger:\mathcal U\to \mathcal V,$6 and cost from $\mathcal G^\dagger:\mathcal U\to \mathcal V,$7 to about $\mathcal G^\dagger:\mathcal U\to \mathcal V,$8 (Hu et al., 10 May 2026). For incompressible flows, efficient implementation uses approximate kernel Fekete points, Schur complement reduction, and streaming block construction, with reported reductions from $\mathcal G^\dagger:\mathcal U\to \mathcal V,$9 to $T:h\mapsto u$00 in preprocessing and from $T:h\mapsto u$01 to $T:h\mapsto u$02 in storage (Sharma et al., 17 Feb 2026). On the modeling side, the direct PDE solution-operator theory is developed for linear PDEs, extensions to nonlinear PDEs and manifold settings are future work, and sharper rates may require stronger spectral or source assumptions (Hu et al., 10 May 2026). The PDE-discovery framework warns that the learned surrogate PDE may be ill-posed and does not guarantee exact physical invariants (Long et al., 2022). Even in strongly structured settings, out-of-distribution robustness is not uniformly favorable: on one backward-facing-step OOD test, kernel methods reported error $T:h\mapsto u$03, compared with Geo-FNO $T:h\mapsto u$04 and Transolver $T:h\mapsto u$05 (Sharma et al., 17 Feb 2026).

Taken together, these works define property-preserving kernel-based operator learning as a technically diverse but coherent research direction. Its common strategy is to encode the desired invariants or structural priors into kernels, recovery maps, basis functions, or convex variational principles, so that the learned operator is constrained by the problem’s geometry before generalization is assessed.