Self-Explainable Operator Learning

• Self-explainable operator learning is a design strategy where prediction models inherently decompose outputs into localized, interpretable contributions.
• It employs frameworks like generalized functional linear models, GP-based schemes, and spectral decompositions to reveal kernels, modes, and uncertainty.
• This approach enables transparent mapping in PDEs and dynamical systems, balancing interpretability with predictive performance trade-offs.

Searching arXiv for papers directly related to self-explainable operator learning and closely adjacent operator-learning paradigms. Self-explainable operator learning concerns operator-learning schemes in which the mathematical objects used to compute a prediction also expose the basis of that prediction. In the most explicit formulation, the operator is written as a linear combination of generalized functional linear models expressed through integral equations, and the additivity of those integrals yields exact spatially localized contributions without invoking a separate explainer (Alishiri et al., 2 Jul 2026). A broader reading is suggested by recent work that represents operators through explicit kernels, bilinear forms, orthonormal modal bases, spectral multipliers, eigenfunctions, selector weights, or uncertainty models, so that explanation is tied to the operator’s native representation rather than added post hoc (Boullé et al., 2023, Kovachki et al., 2024, Mora et al., 2024, Turri et al., 24 May 2025).

1. Conceptual scope

Operator learning studies maps between function spaces, such as $\mathcal{A}:\mathcal{U}\to\mathcal{V}$, $G:U\to V$, or $\Psi^\dagger:U\to V$, where inputs and outputs are functions rather than finite-dimensional vectors (Boullé et al., 2023, Kovachki et al., 2024, Mora et al., 2024). In PDE settings, the operator is often a solution operator: for Poisson, Darcy, Burgers, advection, or heat-conduction problems, it maps forcing terms, coefficient fields, boundary data, or initial conditions to solution fields (Boullé et al., 2023, Kovachki et al., 2024, Huang et al., 12 Oct 2025). In dynamical systems, the relevant object may be an evolution operator or Koopman/transfer operator, acting on observables by conditional expectation,

$$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$

with spectral structure encoding decay rates, oscillation frequencies, coherent structures, and metastable sets (Turri et al., 24 May 2025).

Within this setting, “self-explainable” refers to by-design interpretability. The explicit contrast is with post-hoc XAI methods such as SHAP, Grad-CAM, saliency, or occlusion, which analyze a trained black box from the outside; by contrast, the self-explainable formulation makes the explanation an algebraic consequence of the operator definition itself (Alishiri et al., 2 Jul 2026). This suggests a criterion broader than any single architecture: an operator learner is self-explainable to the extent that its prediction can be decomposed into inspectable basis functions, kernels, spectral factors, selector votes, or uncertainty contributions that retain physical or semantic meaning (Mora et al., 2024, Bistrian, 5 Aug 2025, Bayani, 2022).

2. Mathematical forms of embedded explanation

Several recent formulations make explanation-native operator learning mathematically explicit. One is the generalized functional linear model, where a function-to-function operator on a 2D domain is written as

$$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$

and the function-to-scalar version replaces the output field by a scalar response (Alishiri et al., 2 Jul 2026). Because the mapping is linear in the integral terms and additive across the domain, partitioning the input domain as $\Omega_T=\cup_{i=1}^k\Omega_i$ yields an exact decomposition $\mathbf{u}=\sum_i \mathbf{u}_i$, where each $\mathbf{u}_i$ is the localized contribution of subdomain $\Omega_i$ (Alishiri et al., 2 Jul 2026).

A second formulation represents the operator through an associated real-valued bilinear form,

$$\widetilde{G}(u,\varphi):=[\varphi,G(u)], \qquad G(u)(y)=\widetilde{G}(u,\delta_y),$$

and places a Gaussian Process prior on the scalar map $G:U\to V$0 (Mora et al., 2024). Here explanation is tied to the mean function, product-kernel structure, posterior covariance, and, in hybrid versions, the correction term applied by the GP to a base neural operator (Mora et al., 2024).

A third family uses explicit spectral decompositions. In self-supervised evolution operator learning, the learned latent operator $G:U\to V$1 satisfies

$G:U\to V$2

so eigenvalues determine timescales and frequencies, while eigenfunctions $G:U\to V$3 act as dynamical coordinates (Turri et al., 24 May 2025). In reduced-order Koopman twin models, orthonormal Koopman modes $G:U\to V$4 and temporal coefficients $G:U\to V$5 give

$G:U\to V$6

making the spatial operator linear and directly analyzable (Bistrian, 5 Aug 2025).

More generally, the mathematical literature on operator learning already identifies basis expansions, integral kernels, graph-local kernels, and Fourier multipliers as native operator representations. DeepONet writes the output as $G:U\to V$7, while FNO layers apply spectral convolutions through learned multipliers in Fourier space (Boullé et al., 2023, Kovachki et al., 2024).

Family	Explanatory object	Representative source
Generalized functional linear model	Localized integral contributions $G:U\to V$8 and coefficients $G:U\to V$9	(Alishiri et al., 2 Jul 2026)
GP bilinear operator learning	Bilinear form, kernel correlation structure, posterior uncertainty	(Mora et al., 2024)
Basis/spectral neural operators	Trunk bases, branch coefficients, Fourier multipliers, kernels	(Boullé et al., 2023, Kovachki et al., 2024)
Evolution-operator models	Eigenvalues, eigenfunctions, latent operator matrix $\Psi^\dagger:U\to V$0	(Turri et al., 24 May 2025)
Koopman-ROM twin models	Orthonormal Koopman modes, MAC matrices, modal amplitudes	(Bistrian, 5 Aug 2025)
Selector-based explanation operators	Selector weights, state distances, UCB statistics	(Bayani, 2022)

3. Architectures and learning procedures

The most direct self-explainable architecture replaces deep neural operators with a library of analytic kernels and lifting maps. In that formulation, $\Psi^\dagger:U\to V$1 may be polynomial, exponential, or tanh, the kernels include Gaussian/RBF-like exponentials and polynomial forms, and training reduces to a single convex ridge-regression problem,

$\Psi^\dagger:U\to V$2

with normal equations solved, for example, by GMRES; no backpropagation and no epochs are required (Alishiri et al., 2 Jul 2026). The deliberate removal of nonlinear activations after the integral terms preserves exact additivity, which is the essential source of its built-in explanation mechanism (Alishiri et al., 2 Jul 2026).

GP-based operator learning preserves a probabilistic operator semantics. The GP may use a zero mean or a neural operator mean function, with separable kernels $\Psi^\dagger:U\to V$3 and $\Psi^\dagger:U\to V$4 producing Kronecker-structured covariances and robust maximum-likelihood training. Physics-informed variants further add PDE, boundary-condition, and initial-condition residual terms to the loss, so that the learned operator is constrained simultaneously by data fit and explicit physical equations (Mora et al., 2024).

Self-supervised evolution operator learning uses an encoder-only architecture, a linear predictor $\Psi^\dagger:U\to V$5, and a contrastive $\Psi^\dagger:U\to V$6 density-ratio objective,

$\Psi^\dagger:U\to V$7

which is equal to the Hilbert–Schmidt distance between the true evolution operator and its projected approximation up to a constant under the stated assumptions (Turri et al., 24 May 2025). With optimal $\Psi^\dagger:U\to V$8, minimizing this loss is equivalent to maximizing VAMP-2, so the learned feature space is aligned with dominant singular functions of the operator (Turri et al., 24 May 2025).

Koopman-based self-explainable reduced operators follow a different route. KROD uses randomized SVD on time-shifted snapshot matrices, constructs a projected Koopman matrix $\Psi^\dagger:U\to V$9, forms the Gram matrix $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$0, and defines orthonormal Koopman modes by $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$1 (Bistrian, 5 Aug 2025). Model rank is then selected by Pareto optimization balancing error and correlation, and temporal evolution is modeled by an explainable NLARX network operating on the reduced modal coefficients (Bistrian, 5 Aug 2025).

A related but conceptually distinct strand concerns learning operators over explanation states rather than physical fields. In Fanoos, explanations are stored in states $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$2, operators transform one state into another, and operator selection is learned through a weighted ensemble of selectors combined with UCB. The implementation uses a library of $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$3 operators and $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$4 selectors, of which $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$5 are higher-order selectors (Bayani, 2022). Since selectors are explicit modules—uniform, applicability-based, similarity-based, or multiplicative combinations of other selectors—and weights are updated by a simple rule, the operator policy remains inspectable during learning and deployment (Bayani, 2022).

Attention-enhanced neural operators occupy an intermediate position. SAU-FNO combines FNO, U-Net, and a final self-attention block in order to capture global spectral structure, local high-frequency detail, and long-range spatial dependencies in 3D-IC thermal simulation (Huang et al., 12 Oct 2025). The paper does not present the method as self-explainable, but it explicitly identifies attention maps, multi-scale U-Net features, and spectral filters as natural explanatory artifacts that could be surfaced (Huang et al., 12 Oct 2025). This suggests a distinction between intrinsically explainable operator classes and explanation-friendly neural operators.

4. Intrinsic attribution and diagnostic objects

The defining mechanism of explicit self-explainable operator learning is exact additive attribution. After partitioning the domain, the integral model produces zonal fields $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$6, mean-centered contributions $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$7, and exact reconstruction of the mean-centered full output by summation over zones (Alishiri et al., 2 Jul 2026). Importance is then quantified by Pearson correlation between each zonal contribution and the full output, either for scalar responses or per-sample field responses (Alishiri et al., 2 Jul 2026). Because these quantities are computed from the same operator used for prediction, they are not surrogate attributions.

In selector-based operator learning, the diagnostic objects are different but equally explicit. Fanoos makes the state features themselves interpretable—numbers of predicates, conjunctions, box volumes, and abstraction statistics—and bases selector outputs on distances in that feature space, decayed success rates over similar past states, applicability rules, and UCB exploration bonuses (Bayani, 2022). The autouser adds a transparent notion of what counts as “more abstract” or “less abstract” by evaluating named-predicate volume, box-range volume, number of predicates, number of conjunctions, and number of box-range predicates through a ternary step function and history-adaptive thresholds (Bayani, 2022). The resulting policy can therefore be narrated in terms of similar past explanation states and explicit abstraction criteria.

Evolution-operator methods use spectral diagnostics rather than spatial partitions. The learned eigenfunctions in the self-supervised framework are the primary explanatory objects: they encode slow collective variables, committor-like structure, and coherent modes, while eigenvalues determine implied timescales $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$8 (Turri et al., 24 May 2025). In Koopman-ROM models, orthogonality is verified through MAC matrices, modal amplitudes show when a spatial pattern grows or decays, and input-output validation is performed per mode, making failure modes interpretable at the level of specific coefficients (Bistrian, 5 Aug 2025).

GP-based operator learning contributes a different explanatory layer: calibrated posterior uncertainty. Because the predictor is defined on the bilinear form and uses product kernels over input-function features and output locations, length scales and covariance structure describe how similarity in function space and physical space is translated into confidence, interpolation, and smoothness (Mora et al., 2024). For hybrid GP/NN models, the residual GP correction is itself interpretable as covariance-weighted interpolation of errors left unexplained by the mean operator (Mora et al., 2024).

In neural-operator systems such as SAU-FNO, explanation is structurally available but not intrinsic in the same strict sense. Spatial attention matrices $$(\mathcal{K}f)(x_t)=\mathbb{E}_{y\sim X_{t+1}\mid X_t=x_t}[f(y)],$$9, channel attention maps $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$0, U-Net feature hierarchies, and Fourier-space kernels $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$1 can be visualized and interpreted as long-range couplings, scale-specific thermal features, and learned transfer functions; however, the paper explicitly notes that these artifacts are not yet used as formal explanations and that no PDE-residual loss is imposed (Huang et al., 12 Oct 2025).

5. Empirical domains and observed explanatory patterns

The blood-flow and airfoil studies provide the clearest demonstration of intrinsic attribution. In the aneurysm function-to-scalar task, partitioning the input velocity slice into four zones showed that Zone 3 had the highest average Pearson correlation for Max WSS, $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$2, while Zone 1 dominated Mean WSS with PC $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$3, Zone 2 reached $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$4 for Mean DNWSS, and function-to-function WSS prediction again emphasized Zone 3 (Alishiri et al., 2 Jul 2026). The stated physical interpretation is that the operator most often prioritizes regions with strong velocity-magnitude gradients, consistent with the relation between near-wall velocity gradients and wall shear stress (Alishiri et al., 2 Jul 2026). In the airfoil problem, low-frequency pitching associated $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$5 primarily with Zones 1 and 4 and $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$6 with Zone 2, whereas at higher pitching frequency both $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$7 and $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$8 were dominated by Zone 4, the wake region (Alishiri et al., 2 Jul 2026).

The same paper compares its built-in decomposition with Kernel SHAP, occlusion sensitivity, and Grad-CAM applied to a U-Net surrogate for the blood-flow function-to-function task. All three post-hoc methods identified Zone 3 as the dominant region on average, agreeing qualitatively with the intrinsic zonal decomposition while differing on lower-ranked zones (Alishiri et al., 2 Jul 2026). The comparison supports the claim that built-in explanations can recover the same dominant structures without retraining or perturbation-based attribution.

In Burgers shock-wave experiments, the Koopman-based twin-model framework used $$\mathbf{u}(x,y) = \sum_{n=1}^{N}\sum_{m=1}^{M} w_{n,m} \iint \psi_n(x-\zeta,y-\eta)\,\mathcal{T}_m \mathbf{f}(\zeta,\eta)\,d\zeta\,d\eta,$$9, $\Omega_T=\cup_{i=1}^k\Omega_i$0, and Pareto selection to choose $\Omega_T=\cup_{i=1}^k\Omega_i$1, $\Omega_T=\cup_{i=1}^k\Omega_i$2, and $\Omega_T=\cup_{i=1}^k\Omega_i$3 modes for three increasingly complex initial conditions (Bistrian, 5 Aug 2025). The resulting reduced models achieved correlation coefficient $\Omega_T=\cup_{i=1}^k\Omega_i$4 for all experiments and mean absolute errors of order $\Omega_T=\cup_{i=1}^k\Omega_i$5, while MAC matrices confirmed numerical orthogonality (Bistrian, 5 Aug 2025). Here interpretability is modal rather than zonal: the number of retained modes quantifies intrinsic complexity, and each orthonormal mode can be visualized as a distinct spatial structure associated with shock formation and propagation (Bistrian, 5 Aug 2025).

In high-dimensional dynamical systems, the self-supervised evolution-operator framework produced explanatory eigenfunctions across several scientific domains. For Trp-cage folding, the leading eigenfunction correlated strongly with RMSD from the folded state and radius of gyration, and sparse LASSO regression linked it to a specific set of hydrogen-bond and side-chain interactions (Turri et al., 24 May 2025). For ligand binding, $\Omega_T=\cup_{i=1}^k\Omega_i$6 described the semi-bound to bound transition and $\Omega_T=\cup_{i=1}^k\Omega_i$7 the unbound to bound transition, with a trapped water molecule identified as a kinetic bottleneck (Turri et al., 24 May 2025). For climate, an eleventh eigenfunction recovered ENSO-like variability, activated in the tropical Pacific, and correlated strongly with the Oceanic Niño Index, including the 2023 El Niño event in the validation period (Turri et al., 24 May 2025).

SAU-FNO addresses a different trade-off. On Chip1, it reduced RMSE relative to FNO from $\Omega_T=\cup_{i=1}^k\Omega_i$8 to $\Omega_T=\cup_{i=1}^k\Omega_i$9 at $\mathbf{u}=\sum_i \mathbf{u}_i$0 and from $\mathbf{u}=\sum_i \mathbf{u}_i$1 to $\mathbf{u}=\sum_i \mathbf{u}_i$2 at $\mathbf{u}=\sum_i \mathbf{u}_i$3, while reducing maximum temperature error from approximately $\mathbf{u}=\sum_i \mathbf{u}_i$4–$\mathbf{u}=\sum_i \mathbf{u}_i$5 K to less than $\mathbf{u}=\sum_i \mathbf{u}_i$6 K (Huang et al., 12 Oct 2025). It achieved average prediction time $\mathbf{u}=\sum_i \mathbf{u}_i$7 s, compared with $\mathbf{u}=\sum_i \mathbf{u}_i$8 s for MTA, a reported $\mathbf{u}=\sum_i \mathbf{u}_i$9 speedup (Huang et al., 12 Oct 2025). These are performance results rather than intrinsic-attribution results, but the paper explicitly situates attention maps, multi-scale features, and spectral filters as explanation-friendly structures that could be exposed in future operator learners (Huang et al., 12 Oct 2025).

The Fanoos operator-selection study remains primarily methodological: it outlines an experimental plan using two simple testing domains and random question generation, with quantitative results described as ongoing (Bayani, 2022). Its significance for self-explainability lies in showing that operator learning can target explanation-management policies themselves, not only physical PDE maps (Bayani, 2022).

6. Limitations, tensions, and research directions

A recurrent limitation is the trade-off between transparency and predictive flexibility. The generalized functional linear model is intrinsically decomposable, but compared with FNO baselines it is typically about one order of magnitude worse in mean absolute error, and its feature-library construction scales as $\mathbf{u}_i$0 for 2D $\mathbf{u}_i$1 grids (Alishiri et al., 2 Jul 2026). The same framework is most interpretable when inputs have coherent spatial patterns and becomes harder to interpret on highly irregular or chaotic synthetic PDE data (Alishiri et al., 2 Jul 2026).

Spectral and Koopman-based approaches inherit their own assumptions. The self-supervised evolution-operator method relies on a Hilbert–Schmidt assumption for exact identification of the loss with operator regression, even though the paper notes that this may fail for many deterministic systems; climate experiments also showed training instability and possible overfitting (Turri et al., 24 May 2025). The reduced twin-model framework is demonstrated on 1D viscous Burgers’ equation with $\mathbf{u}_i$2 snapshots and identity observables, and it does not yet implement formal XAI tools inside the NLARX component (Bistrian, 5 Aug 2025).

Selector-based operator learning in Fanoos is transparent but heuristic in several respects. The paper explicitly notes that vanilla UCB is used in a nonstationary setting, that the autouser reward shaping may not match all human preferences, that selectors do not yet use domain- and question-type information explicitly, and that there is no current mechanism to avoid or reason explicitly about cycles (Bayani, 2022).

GP-based operator learning offers uncertainty and functional-analytic clarity, but kernel design, scalability, and joint optimization of neural means with kernel hyperparameters remain difficult; the paper reports that deep kernels did not yield consistent gains and that rich mean functions can leave ill-conditioned residual covariance structure (Mora et al., 2024). Attention-enhanced neural operators show strong accuracy and speed, but their explanations remain post-hoc or potential rather than guaranteed faithful, because the papers do not visualize the learned attention maps as formal explanatory objects and do not impose PDE-residual constraints (Huang et al., 12 Oct 2025).

Across these strands, a common future direction is to preserve the explanatory status of kernels, modes, selectors, or uncertainty while increasing expressivity. The explicit suggestions in the literature include more principled contextual or Bayesian bandits for selector policies, richer yet still structured nonlinear operators, deeper integration of explainability techniques into modal temporal models, nontrivial multi-output kernels, physics-aware priors, sparse or block-structured latent operators, symbolic or sparse regression on learned eigenfunctions, and explicit meta-explanation layers that surface case-based or mode-based rationales directly to practitioners (Bayani, 2022, Mora et al., 2024, Turri et al., 24 May 2025, Bistrian, 5 Aug 2025). In that sense, self-explainable operator learning is best understood not as a single architecture, but as a design principle: the operator representation should itself serve as the explanation.