Graph-Based Neural Operators

• Graph-based neural operators are neural architectures that use graph-supported operators to transform signals on nodes, meshes, or irregular samples.
• They discretize continuum operators via message passing, spectral filtering, and learnable propagation schemes to handle diverse resolutions and topologies.
• Research shows their utility in PDE surrogate modeling, classification, and numerical discretization, while highlighting trade-offs in computational efficiency and design choices.

Graph-based neural operators are neural architectures in which the core transformation is realized by a graph-supported operator acting on signals, fields, or functions defined on nodes, meshes, or irregular point samples. In PDE settings, they are used to approximate mappings between function spaces and to transfer across resolutions, discretizations, and mesh types; in graph representation learning, the propagation operator itself is treated as a design variable that controls smoothing, topology sensitivity, spectral response, and embedding geometry (Li et al., 2020, Jeong et al., 2022). Across this literature, the graph serves variously as a discretization of a continuous domain, a support for integral-kernel message passing, and a parameterized shift or convolution operator, yielding families that include graph kernel networks, multiscale graph neural operators, learnable path operators, centrality-aware shifts, compact operator spaces, and scalable spectral transformer constructions (Migus et al., 2022, Eliasof et al., 2022, Abbahaddou et al., 2024, Rigotti et al., 17 Mar 2026).

1. Operator-theoretic foundations

The operator-learning formulation is explicit in early graph neural operator work on PDEs. Rather than learning a map between fixed-dimensional vectors, the target is an operator $F:A\to U$ between function spaces, with training samples $\{a_j,u_j\}_{j=1}^N$ satisfying $u_j=F(a_j)$. In the graph kernel network formulation, the latent update is written in continuum form as

$$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$

and then discretized by message passing on a neighborhood graph,

$$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$

with edge features $e(x,y)=(x,y,a(x),a(y))$ (Li et al., 2020). The same work emphasizes that the graph is a numerical realization of a continuum operator, which is why one parameter set can be reused across different resolutions and mesh types.

A closely related formulation appears in neural integral operators discretized on graphs. There the generic update is

$$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$

and the graph approximation replaces the integral by a local neighborhood average over $\mathcal{N}(x)$. This produces a message-passing graph neural network whose node features encode both physical state and position, and whose graph approximates local physical interactions on the domain (Migus et al., 2022). In this sense, graph-based neural operators are not merely GNNs applied to scientific data; they are graph discretizations of learned integral operators.

The same operator viewpoint also appears in graph signal processing formulations based on spectral filtering. Polynomial filters of the normalized Laplacian,

$$\mathbf{f}_{\boldsymbol{\theta}} *_G \mathbf{x} = \sum_{i=0}^{k} \theta_i \mathbf{L}^i \mathbf{x},$$

lead to single-layer graph convolution operators such as

$$\mathbf{Y}=softmax(\mathbf{S}^k\mathbf{X}\mathbf{\Theta})$$

and more general linear graph convolutions

$\{a_j,u_j\}_{j=1}^N$0

where the graph operator is an explicit learnable filter rather than a deep stack of nonlinear propagations (Pasa et al., 2021). A still earlier precursor generalized convolution to irregular domains through a local allocation map and weight sharing, writing the operator as

$\{a_j,u_j\}_{j=1}^N$1

and showing that it reduces exactly to standard convolution on regular grids (Vialatte et al., 2016). These formulations established the basic principle that graph structure can support CNN-like or operator-like transformations without requiring Euclidean regularity.

2. Learnable propagation and representation geometry

A major line of work treats the graph propagation rule itself as the principal object of study. One influential parameterization defines the family of symmetrized spatial convolutions

$\{a_j,u_j\}_{j=1}^N$2

with $\{a_j,u_j\}_{j=1}^N$3, together with the row-normalized family

$\{a_j,u_j\}_{j=1}^N$4

The central distinction is that $\{a_j,u_j\}_{j=1}^N$5 is a weighted sum over neighbors, whereas the row-normalized version is a convex combination of neighbor features. The reported geometric consequence is that row-normalized outputs lie in the convex hull of the neighborhood and act as local averaging, while symmetrized operators preserve magnitude differences and topological signals. The same analysis identifies a phase transition at $\{a_j,u_j\}_{j=1}^N$6: for $\{a_j,u_j\}_{j=1}^N$7, high-degree nodes can spread farther from the origin, whereas for $\{a_j,u_j\}_{j=1}^N$8, they are pulled toward the origin; $\{a_j,u_j\}_{j=1}^N$9 acts as additional mass in the degree term and interpolates between a standard GNN regime and an MLP-like regime when $u_j=F(a_j)$0 is large (Jeong et al., 2022).

This operator-level distinction is not only theoretical. On eight node-classification datasets, changing $u_j=F(a_j)$1 and $u_j=F(a_j)$2 while leaving the rest of the architecture fixed materially changed performance. Symmetric operators showed strong dependence on $u_j=F(a_j)$3, whereas row-normalized operators were much more robust. On homophilic graphs, symmetric operators could retain more useful structure; on low-homophily graphs such as Cornell and Wisconsin, their gains were weaker and row-normalized operators could be more robust (Jeong et al., 2022). The broader implication is that aggregation determines what the network preserves and how latent space is organized.

A second route to operator learning dispenses with fixed Laplacian propagation and learns the spatial operator from graph paths. In pathGCN, a path $u_j=F(a_j)$4 beginning at node $u_j=F(a_j)$5 is filtered by a learned kernel $u_j=F(a_j)$6 through

$u_j=F(a_j)$7

and the final operator averages over $u_j=F(a_j)$8 sampled random paths. Its deterministic expectation becomes

$u_j=F(a_j)$9

This yields a polynomial filter in the transition matrix $$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$0, but with coefficients learned from path sampling rather than prescribed by a Laplacian-based template. The method explicitly argues that learning the spatial operator helps avoid over-smoothing, and its ablations rank the variants as $$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$1 (Eliasof et al., 2022).

Taken together, these results recast graph convolution as a tunable geometric mechanism. Whether the operator is a convex averaging rule, a degree-sensitive symmetrized sum, or a path-based polynomial in a transition matrix materially changes the embedding geometry, the spectral bias, and the relation between topology and features.

3. Spectral, shift, and invariant operator spaces

Another strand of the literature enlarges the admissible graph operator space by parameterizing graph shift operators directly. Centrality Graph Shift Operators replace degree normalization by a diagonal centrality matrix $$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$2, giving the basic operator

$$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$3

Here $$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$4 may be degree, PageRank, $$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$5-core number, or fixed-length walk count. The same work defines a learnable family

$$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$6

and proves that $$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$7 is self-adjoint with respect to a weighted inner product and therefore diagonalizable with real eigenvalues. Empirically, the work reports that no single centrality is best on all datasets, that $$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$8-core and PageRank often outperform walk-count centrality, and that combining local and global CGSOs can improve performance (Abbahaddou et al., 2024).

A more compact unification is provided by the Generalized Linear Graph Representation,

$$v_{t+1}(x)=\sigma\!\left( W v_t(x)+\int_D \kappa_\phi\!\big(x,y,a(x),a(y)\big)\,v_t(y)\,\nu_x(dy) \right),$$9

This bounded two-parameter family contains adjacency, degree, Laplacian, signless Laplacian up to scaling, Nikiforov’s $$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$0, Wang’s $$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$1, and transition-type operators inside a compact rectangular domain. Its quadratic form is decomposed into a balance between local smoothness and global degree-weighted energy, and the paper derives Lipschitz spectral perturbation bounds as well as graph-dependent sufficient conditions for positive semi-definiteness. Adaptive GLGR Convolution then makes $$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$2 trainable layerwise in end-to-end GNNs (Zhao et al., 12 Jun 2026). This replaces isolated hand-designed matrices by an explicit compact operator space.

The operator space can also become noncommutative. In graph-tuple neural networks, several graphs share a common vertex set and each layer is parameterized by matrices of noncommutative polynomials $$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$3. A monomial $$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$4 is evaluated as the ordered composition $$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$5, so the order of modalities matters whenever the corresponding operators do not commute. The same work derives stability bounds via expansion constants $$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$6 and $$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$7, develops a graphon-tuple limit theory, and proves a universal transferability theorem for convergent graph-tuple sequences (Velasco et al., 2024). This extends graph neural operator theory from single propagation matrices to several simultaneous graph operators.

Scalability and invariance introduce further structure. Gauge-Invariant Spectral Transformers approximate Laplacian-eigenmap-like embeddings with FastRP random projections and then force all meaningful attention scores to depend only on inner products of the projected embeddings. The architecture achieves end-to-end complexity

$$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$8

and is designed so that arbitrary sign flips, rotations in degenerate eigenspaces, or solver-dependent spectral choices do not alter the attention kernel in a gauge-dependent way (Rigotti et al., 17 Mar 2026). At a more geometric extreme, Finslerian graph neural networks replace Laplace–Beltrami-type isotropic operators with graph discretizations of the Finsler Laplacian

$$v_{t+1}(x)=\sigma\Bigl(Wv_t(x)+\frac{1}{|N(x)|}\sum_{y\in N(x)}\kappa_\phi(e(x,y))\,v_t(y)\Bigr),$$9

and show that the resulting discrete operator converges to the manifold operator as sampling density increases (Roddenberry et al., 15 Jun 2026). Both cases illustrate that graph-based neural operators increasingly encode invariance and continuum geometry directly into the operator class.

4. Multiscale and PDE-surrogate formulations

Graph-based neural operators are especially prominent in PDE surrogate modeling. A central multiscale construction starts from the graph neural operator and its multipole extension, in which long-range interactions are approximated through coarse-to-fine kernel compositions such as

$e(x,y)=(x,y,a(x),a(y))$0

This yields V-, F-, and W-cycle architectures built from downscale kernels, intra-scale kernels, and upscale kernels. On 2D Darcy flow and 1D viscous Burgers’ equation, these multiscale message-passing graph neural operators reported about 80–90% reduction in training error relative to single-scale baselines and 10–50% improvement in test error for MGNO variants, although the more expressive F- and W-cycle variants also increased optimization difficulty and motivated orthogonal initialization (Migus et al., 2022).

Irregular-domain operator learning motivates hybrid spectral-graph designs. GOLA constructs a graph from irregularly sampled spatial points, uses a learnable Fourier encoder with complex frequencies $e(x,y)=(x,y,a(x),a(y))$1, and then applies an attention-enhanced GNN decoder. The operator is written as

$e(x,y)=(x,y,a(x),a(y))$2

with training by relative $e(x,y)=(x,y,a(x),a(y))$3 error on Darcy Flow, Advection, Eikonal, and Nonlinear Diffusion. With only 30 training samples, the reported test errors were $e(x,y)=(x,y,a(x),a(y))$4 for Darcy Flow, $e(x,y)=(x,y,a(x),a(y))$5 for Advection, $e(x,y)=(x,y,a(x),a(y))$6 for Eikonal, and $e(x,y)=(x,y,a(x),a(y))$7 for Nonlinear Diffusion, each outperforming the listed DeepONet, FNO, GKN, or GCN baselines in that regime (Li et al., 25 May 2025). The same study emphasizes that graph-based operator learning is especially natural when spatial samples are sparse, nonuniform, or lie on unstructured meshes.

Time dependence introduces a different hybridization. Instead of learning a global solution operator directly, Graph Neural Simulators model the instantaneous time derivative on a graph and advance the state by explicit Euler,

$e(x,y)=(x,y,a(x),a(y))$8

On 2D Burgers’ scalar, 2D coupled Burgers’ vector, and 2D Allen–Cahn, the reported outcome was under $e(x,y)=(x,y,a(x),a(y))$9 relative $$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$0 error with only 30 training trajectories across all three PDEs. Averaged across all cases, GNS reduced autoregressive error by $$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$1 relative to FNO-AR and $$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$2 relative to DON-AR (Nayak et al., 7 Sep 2025). The paper explicitly presents this as evidence that graph-based local inductive bias plus explicit time stepping can be more data-efficient than standard neural operators in low-data settings, rather than as a universal replacement for them.

Physics-informed training on arbitrary geometries requires differentiable graph operators. The mollified graph neural operator replaces the hard neighborhood indicator $$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$3 in GNO with a smooth weight $$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$4, enabling automatic differentiation through the operator

$$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$5

This was reported as the first method in that setting to compute exact gradients on arbitrary geometries by autograd. On a regular-grid PDE example, mGNO paired with autograd reduced the $$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$6 relative data error by $$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$7 compared to finite differences, although training was slower; on unstructured point clouds, it yielded errors consistently 2 orders of magnitude lower than Meta-PDE for comparable runtimes and delivered speedups from 1 to 3 orders of magnitude compared to the numerical solver for similar accuracy (Lin et al., 11 Apr 2025). Graph-based neural operators thus function not only as forward surrogates but also as differentiable operators inside physics losses and inverse-design loops.

5. Classification, structural learning, and numerical discretization

Although PDE surrogates dominate the term’s contemporary use, graph-based neural operators also appear in supervised classification and in learned numerical building blocks. In spatial transcriptomics, the tissue is treated as a discretized spatial field over a graph of spots, and the graph neural operator is written as

$$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$8

On a mouse brain atlas with about 34,000 spots, 23,000 genes, and 75 tissue samples, with labels collapsed from 15 brain-region classes to 3 high-level classes, the GraphPDE model achieved $$v_{t+1}(x) = \sigma\Big( Wv_t(x)+b(x)+\int_D \kappa_\phi(x,x')\,v_t(x')\,d\nu(x') \Big),$$9 accuracy and $\mathcal{N}(x)$0 F1, reported as nearly $\mathcal{N}(x)$1 F1 and higher than the other graph and baseline approaches on the primary class-balanced metric (Ahmed et al., 2023). This use of the operator viewpoint is explicitly tied to learning mappings between discretized samples of continuous spatial structure.

Graph kernels provide a different notion of graph-based operator. In Graph Kernel Neural Networks, each node-centered radius-$\mathcal{N}(x)$2 neighborhood subgraph $\mathcal{N}(x)$3 is compared against learned structural masks $\mathcal{N}(x)$4 through

$\mathcal{N}(x)$5

The resulting Graph Kernel Convolution does not aggregate neighbor features by message passing; it compares local graph structure directly to learned graph templates. The masks are optimized by discrete randomized descent, and the architecture can plug in graph kernels such as WL, WL with Optimal Assignment, Graphlet, Propagation, or Pyramid match. On graph classification, the paper reports competitive performance, including $\mathcal{N}(x)$6 on MUTAG and $\mathcal{N}(x)$7 on PROTEINS (Cosmo et al., 2021). Here the operator is fully structural and remains outside Laplacian-based or attention-based propagation.

At a lower numerical level, graph-based operator learning can target discrete differential stencils rather than solution maps. The mesh-free operator framework

$\mathcal{N}(x)$8

is parameterized by a GNN that maps normalized relative positions $\mathcal{N}(x)$9 to stencil weights $$\mathbf{f}_{\boldsymbol{\theta}} *_G \mathbf{x} = \sum_{i=0}^{k} \theta_i \mathbf{L}^i \mathbf{x},$$0, and trained self-supervised by polynomial moment constraints derived from truncated Taylor expansions. The resulting operators are reported as resolution-agnostic and reusable across particle configurations and governing equations. In convergence tests, they improved over SPH and offered a favorable accuracy-cost trade-off relative to a representative high-order mesh-free method, while also being inserted into a weakly compressible Navier–Stokes solver (Starepravo et al., 25 Mar 2026). This broadens the meaning of graph-based neural operator from end-to-end surrogate map to learned numerical primitive.

6. Benchmarking, trade-offs, and recurrent technical themes

Large-scale engineering benchmarks make clear that graph-based neural operators are useful but not uniformly dominant. A comparative study on six 3D industry-scale engineering design datasets treats GNO and EA-GNO as representative graph-based operators because they can operate directly on FEM meshes without remeshing or new graph construction. The same benchmark also reports that graph-based models are generally not the top performers, are computationally expensive because of graph construction, neighbor search, and repeated message passing, and can be much slower than branch-trunk, grid-based, or point-based alternatives. Representative per-epoch training times reported for GNO were 45.0 s on Heat sink, 230 s on Bracket, 372 s on Bracket-time, 1200 s on JEB, and 3600 s on DrivAerNet; EA-GNO was slower still on several datasets (Zhong et al., 7 Oct 2025). Mesh-native deployment is therefore a practical advantage, but not a guarantee of superior accuracy or scalability.

A recurring misconception is that graph-based neural operators form a single architectural class. The literature instead spans continuum integral operators discretized by message passing, learned graph shifts and propagation families, graph-kernel structural operators, graph-tuple filters in several noncommuting modalities, differentiable mollified kernels for physics-informed learning, and geometry-constrained discretizations of nonlinear differential operators (Li et al., 2020, Velasco et al., 2024, Lin et al., 11 Apr 2025, Roddenberry et al., 15 Jun 2026). Another misconception is that graph locality by itself determines performance. The operator choice can privilege local averaging, topology-sensitive weighted sums, global centrality, random-path filtering, long-range spectral attention, or multiscale coarse-to-fine transfer, and the reported best choice is consistently task-dependent (Jeong et al., 2022, Abbahaddou et al., 2024, Eliasof et al., 2022, Rigotti et al., 17 Mar 2026).

The aggregate picture is therefore one of operator design rather than architecture branding. Symmetrized versus row-normalized aggregation changes latent geometry and degree sensitivity; compact operator families such as GLGR turn propagation into a bounded search space; gauge-invariant spectral constructions are introduced precisely because naive spectral coordinates can break discretization transfer; mollification is introduced because hard neighborhood indicators obstruct automatic differentiation; and low-data PDE results favor graph-based local inductive bias when causal local structure matters (Zhao et al., 12 Jun 2026, Rigotti et al., 17 Mar 2026, Lin et al., 11 Apr 2025, Nayak et al., 7 Sep 2025). This suggests that graph-based neural operators are best understood as a research program centered on how graph-supported operators should be parameterized, discretized, regularized, and matched to the geometry of the underlying problem.