Graph Fourier Neural Kernels (G-FuNK): Learning Solutions of Nonlinear Diffusive Parametric PDEs on Multiple Domains

## Abstract

Abstract Predicting the time-dependent dynamics of complex systems governed by non-linear partial differential equations (PDEs), with varying parameters and domains, is a difficult problem that is motivated by applications in many fields. We introduce a novel family of neural operators based on a G raph F o u rier N eural K ernel (G-FuNK), for learning solution generators of nonlinear PDEs with varying coefficients, across multiple domains, for which the highest-order term in the PDE is diffusive.
G-FuNKs are constructed by combining components that are parameter- and domain-adapted, with others that are not. The latter components are learned from training data, using a variation of Fourier Neural Operators, and are transferred directly across parameters and domains. The former, parameter- and domain-adapted components are constructed as soon as a parameter and a domain on which the PDE needs to be solved are given. They are obtained by constructing a weighted graph on the (discretized) domain, with weights chosen so that the Laplacian on that weighted graph approximates the highest order, diffusive term in the generator of the PDE, which is parameter- and domain-specific, and satisfies the boundary conditions. This approach proves to be a natural way to embed geometric and directionally-dependent information about the domains, allowing for improved generalization to new test domains without need for retraining. Finally, we equip G-FuNK with an integrated ordinary differential equation (ODE) solver to enable the temporal evolution of the system’s state.
Our experiments demonstrate G-FuNK’s ability to accurately approximate heat, reaction diffusion, and cardiac electrophysiology equations on multiple geometries and varying anisotropic diffusivity fields.
We achieve low relative errors on unseen domains and fiber fields, significantly speeding up prediction capabilities compared to traditional finite-element solvers.

## 1 Introduction

##### Neural Operators for PDEs

In scientific machine learning, data-driven deep learning methods aim at predicting the solutions of partial differential equations (PDEs), avoiding the computationally expensive numerical integration methods needed for large-scale simulations. This is especially beneficial for applications like domain optimization and precision medicine, where multiple PDEs need to be solved for varying parameters or domains, requiring significant computational resources. Neural operators learn mappings between high-dimensional function spaces, allowing them to generalize across a family of PDEs without retraining for varying parameters or conditions Kovachki et al. (). Mathematically, a neural operator $\mathcal{N}_{\theta}$ is defined as a mapping from an input domain function space $\mathcal{A}(\Omega_{\alpha};\mathbb{R}^{d_{a}})$ to an output target function space $\mathcal{U}(\Omega_{\alpha};\mathbb{R}^{d_{u}})$, represented as

$$ $\mathcal{N}_{\theta}:\mathcal{A}(\Omega_{\alpha};\mathbb{R}^{d_{a}}) \rightarrow\mathcal{U}(\Omega_{\alpha};\mathbb{R}^{d_{u}}),$ (1) $$

where $\theta\in\Theta$ denotes the neural operator’s parameters, and $\Omega_{\alpha}\subset\mathbb{R}^{d}$ (or a $d$-dimensional manifold) represents the spatial domain on which the functions are defined and $\alpha\in\mathscr{A}$ denotes the shape of the domain. The dimensions $d_{a}$ and $d_{u}$ denote the respective sizes of the input and output function spaces, often subspaces of Sobolev spaces. Neural operators can be applied to various problems such as those described by continuous functions $C(\Omega;\mathbb{R}^{d_{a}})$ or Sobolev spaces $H^{s}(\Omega;\mathbb{R}^{d_{a}})$ for some $s\geq 0$. The neural operator $\mathcal{N}_{\theta}$ is an approximation of a true target operator $\mathcal{N}$ (e.g., the solution operator of a PDE) obtained by training on input-output function pairs $\{a_{i},u_{i}\}_{i=1}^{m}$, where $a_{i}\in\mathcal{A}$ and $u_{i}=\mathcal{N}(a_{i})\in\mathcal{U}$. These pairs could be simulation data representing a known, high-fidelity numerical approximation of the PDE. In cardiac electrophysiology, for example, $\mathcal{A}$ might represent the space of initial electrical activation patterns across cardiac tissue, whereas $\mathcal{U}$ could correspond to the resultant electrical potential fields over time.

##### Problem Setup

In this context, we focus on the family of second-order nonlinear differential operators, governing the evolution of a function $u:(\mathbf{x},t)\rightarrow\mathbb{R}$, $\mathbf{x}\in\Omega_{\alpha}$, $t\geq 0$, in the form,

$$ $\left\{\begin{aligned} \partial_{t}u(\mathbf{x},t)&=\nabla_{\mathbf{x}}\cdot( \mathbf{K}(\mathbf{x})\nabla_{\mathbf{x}}u(\mathbf{x},t))+\mathcal{S}(\,u( \mathbf{x},t),\mathbf{x},\nabla_{\mathbf{x}}u(\mathbf{x},t)\,)\,=\mathcal{N}( \,u(\mathbf{x},t),\mathbf{K}(\mathbf{x}),\mathbf{x},t\,),\\ \partial_{\mathbf{n}(\mathbf{x})}u(\mathbf{x})&=0\>\,,\quad\forall\,\mathbf{x} \>\in\partial\Omega_{\alpha}\,\end{aligned}\right.$ (2) $$

where $\mathbf{K}$ is a diffusion tensor field on the domain satisfying uniform ellipticity, $\mathcal{N}:\mathbb{R}^{3d}\times\mathbb{R}^{d_{\mathbf{K}}}\rightarrow\mathbb
{R}$ is a vector function that can encompass source terms, sinks, or other linear and nonlinear interactions within the term $\mathcal{S}$, in this work we use no-flux or no-diffusive-flux Neumann boundary conditions, where $\mathbf{n}(\mathbf{x})$ is the normal to $\partial\Omega_{\alpha}$ at $\mathbf{x}$, and $\partial\Omega_{\alpha}$ is the boundary of the spatial domain $\Omega_{\alpha}$ (in $\mathbb{R}^{d}$ or a $d$-dimensional manifold $\mathcal{M}^{d}$) within a family of domains $\{\Omega_{\alpha}\}_{\alpha\in\mathscr{A}}$. The dependence of the solution to (2) on the domain $\Omega_{\alpha}$ is often complex.
Note from 2 that we are focusing here on the semilinear case, i.e. $\mathcal{N}$ is linear in its highest-order term (the diffusive component), and in general nonlinear in the lower-order terms.

The training data consists of trajectories $\{u^{(m)}(x_{i},t_{\ell})\}_{m=1,i=1,\ell=1}^{M,n_{\alpha^{(m)}},L}$, where $x_{i}$ is sampled on a graph $\mathcal{G}_{\alpha^{(m)}}$ (a mesh discretizing $\Omega_{\alpha^{(m)}}$), and time points $0=t_{0}<t_{1}<\dots<t_{L}=T$ are sampled at a fine timescale, together with the parameter $\mathbf{K}^{(m)}$ and the graph $\mathcal{G}_{\alpha^{(m)}}$.
In particular, note that the observation at $t=0$ says that we are given the initial condition, and that different trajectories in the training set may correspond to different (given) values of both the parameter $\mathbf{K}$ and the domain $\mathcal{G}_{\alpha}$.

The objective (in contrast to equation (1)) is to learn from the training data a neural operator $\mathcal{N}_{\theta}$, with $\theta$ denoting the neural network parameters, that approximates the generator of the solutions of (2) by learning a mapping which takes any given state $u(x,t_{\ell})$, parameters $\mathbf{K}$ and domain $\mathcal{G}_{\alpha}$, and outputs an approximation of the left-hand side of (2). Solutions of the PDE may be approximated by integrating

$$ $\partial_{t}u(x,t)\approx\mathcal{N}_{\theta}(u(x,t);\mathcal{G}_{\alpha}, \mathbf{K})\,,$ (3) $$

with $x\in\mathcal{G}_{\alpha}$, on the whole interval $[0,T]$, given any initial condition $u(\cdot,0)$ even with changes in parameters and spatial domain. This solver has the chance of being more efficient than one for equation 2, for example by employing larger time steps, by efficiently incorporating the shape of $\Omega_{\alpha}$ and the boundary conditions, or by reducing the dimensionality of the problem when $\mathbf{x}$ is high-dimensional.

##### Reduced Models and Homogenization

In fact, we are very much interested in situations where the “true” system is driven by equations more general than equation 2, with many more variables, possibly evolving at multiple time scales; this will be the case in examples 2 and 3, where besides the spatial variables and the unknown action potential $u$, there are many other quantities $\mathbf{v}$ evolving in space and time (e.g. modeling ionic concentrations that affect the evolution of $u$, and these are driven by a large system of ODEs).
In these situations, given observations from such systems, our estimated equations equation 3 can be viewed as performing dimension reduction on the variables (e.g. by keeping only the spatial variables) and homogenization, obtaining an effective equation for $u$, without tracking the evolution of the other quantities $\mathbf{v}$.

##### Related Works

Several studies have demonstrated the ability of different types of neural networks to approximate nonlinear operators with strong theoretical guarantees (e.g., PDE solution operator) Chen & Chen (); Lu et al. (); Kovachki et al. (). Fourier Neural Operators (FNOs) and Deep Operator Networks (DeepONets) and their variants are among the most popular computational frameworks under the umbrella of neural operators Lu et al. (); Li et al. (). Vanilla DeepONets, however, are restricted to a fixed grid and resolution, and thus do not generalize well on unseen domains. Yin et al. () propose a work-around spatial transformations of multiple domains to a universal domain where DeepONets can be used to learn a latent operator, the evaluation of which can be mapped back to the target domains via the inverse for evaluation. FNOs, while robust to variations in grid resolution, require regularly spaced square grids to perform FFT Li et al. (). As an improvement, Li et al. () propose a Geometry-Aware Fourier Neural Operator (Geo-FNO), which learns an additional deformation step to transform irregular grids, achieving greater accuracy than interpolating the solution on a uniform grid. However, they explore only a relatively small subset of possible transformations and do not investigate generalization performance on differently shaped grids. For large-scale 3D PDEs, Li et al. () also propose a more general way of adapting to arbitrary domains and irregular grids with Geometry Informed Neural Operators (GINO) which incorporates information about the geometry into the learning process, however, they mention that the trained model is limited to certain family of shapes and to generalize well, abundant, geometrically varying training samples are necessary, a luxury that is not available in real world settings such as computational medicine. Several works have also shown promise with regards to predicting full solution trajectories over a time interval, but typically do not carry the same generalizability over multiple domains as the aforementioned works Chen et al. (); Zhang et al. (); He et al. (); Regazzoni et al. ().

Graph Neural Networks (GNNs) are a versatile class of neural network architectures designed to work directly with graph-structured data. These networks can make predictions at the level of nodes, edges, or entire graphs, typically utilizing a method known as message passing Gilmer et al. (). In this method, neighboring nodes exchange messages, and the aggregated information from these messages is used to update the states of subsequent nodes. GNNs are particularly well-suited for irregular grid data and can be effectively applied to training data generated by finite-element method simulations (considered ground truth). Li et al. () proposed a message-passing operator (using edge-conditioned graph convolutions Simonovsky & Komodakis ()) to learn solution operators for PDEs. Iakolev et al. also utilized a message-passing scheme to learn isotropic PDEs Iakovlev et al. (). Other studies have considered learning time-dependent PDEs using graph-based approaches Li et al. (); Behmanesh et al. (); Pilva & Zareei (), but none of these works combine learning with both varying domains and varying parameters in the PDE. Furthermore, the message-passing frameworks they employed in these works involve edge aggregation schemes that limit the learning of direction-dependent information, which is crucial in the case, for example, of anisotropic diffusions. Finally, several studies have shown that spectral GNNs possess superior expressiveness and interpretability, capturing global information more effectively than spatial GNNs Bo et al. (); Defferrard et al. (); Wang & Zhang (); Yang et al. ().

##### Contributions

We propose a new family of neural operators consisting of novel Graph Fourier Neural Kernel (G-FuNK) layers for learning generators of solutions to time-dependent PDEs with varying coefficients across multiple domains as in equation (2). This method combines ideas from both GNNs and FNOs: it leverages the spectral domain of graphs, unlike GNNs that rely on localized message passing, and is suitable on general graphs rather than relying on grids, as FNOs.
G-FuNK employs the Graph Fourier Transform (GFT) to convert functions on graphs into functions on the graph spectral domain. This facilitates global information integration across the entire graph, allowing G-FuNK to capture both local and global dependencies effectively, as FNOs do in the classical Euclidean grid setting.

G-FuNK integrates components that are both parameter- and domain-adapted with others that are not. The non-adapted components are learned from training data, and can be directly transferred across different parameters and domains, independent of the mesh resolution or the time step size.
The adapted components are constructed for specific parameters and domains by forming a parameter-dependent, tailored weighted graph on the discretized domain, such that the corresponding graph Laplacian approximates the highest-order diffusive term in the PDE with those specific values of the parameter, on that domain.

We remark that here we do wish to predict entire trajectories of temporal dynamics on $[0,T]$, which is difficult as inaccuracies can lead to dramatic accumulation of errors, depending on the PDE. Many existing methods Li et al. (; ; ) are either restricted by design on learning the map from initial conditions (or several observation points) to a solution at a fixed time $t$, or while they could in principle generate solutions at all times, experiments in this setting are not reported in the corresponding papers.
With this goal in mind, our framework is equipped with an integrated ODE solver to enable generating the temporal evolution of the dynamics of the system. This allows for efficient, mesh-independent prediction of the system’s solutions from new initial conditions, with new values of the parameter $\mathbf{K}$, on a new domain $\Omega_{\alpha}$. To our knowledge, there is no current neural operator framework which naturally handles directionally dependent information of anisotropic domains while predicting time-evolving trajectories of the solution operator on multiple geometries.

We demonstrate G-FuNK’s capabilities by first considering the basic example of the anisotropic heat equation on the square, where the tensor field $\mathbf{K}(\mathbf{x})$ reflects the directional dependence of diffusivity. We then move to nonlinear, more complex models of reaction-diffusion equations for cardiac electrophysiology (EP), first on families of rectangular domains, then on 3D geometries representing the outer surfaces of real human left atrial (LA) chambers from computed tomography (CT) data. The PDE operators on the patient-specific LA chambers involve effects from high-dimensional coupled ODEs.
Our method accurately captures the complex dynamics of the transmembrane potential on varying geometries and anisotropic diffusivity (fiber) fields.
Our experiments show that G-FuNK can significantly speed up patient-specific cardiac electrophysiology simulations, which are used for making real-time quantitatively informed clinical decisions Boyle et al. (); Trayanova et al. (), achieving low relative errors on unseen domains.

## 2 G-FuNK: G raph F o u rier N eural K ernel

### 2.1 Graphs, Laplacians, and the Graph Fourier Transform

Let $\mathcal{G}_{\alpha}(\mathcal{V},\mathcal{E})$ be an undirected graph where $\mathcal{V}=\{x_{i}\}_{i=1}^{n_{\alpha}}$ is the set of $n_{\alpha}$ nodes (or vertices) constructed from a finite element mesh discretization or down-sampled point cloud representation of $\Omega_{\alpha}$ and $\mathcal{E}\subseteq\{(x_{i},x_{j})|x_{i},x_{j}\in\mathcal{V}\}$ is the set of edges that connects pairs of nodes which
are defined as the k-nearest neighbors of each $x_{i}$.
Let $W\in\mathbb{R}^{n_{\alpha}\times n_{\alpha}}$ be the weighted adjacency matrix, with $W_{ij}=W_{ji}\geq 0$, with strict inequality if and only if $(x_{i},x_{j})\in\mathcal{E}$, and $\mathrm{deg}$ is the $n_{\alpha}\times n_{\alpha}$ diagonal degree matrix of the graph, defined by $\mathrm{deg}_{ii}:=\sum_{j}W_{ij}$.
The undirected graph Laplacian is a positive semi-definite matrix with a spectral decomposition given by

$$ $\displaystyle\mathscr{L}_{\mathcal{G}_{\alpha}}:=\mathrm{deg}-W=\Psi\Lambda \Psi^{T},$ (4) $$

where $\Psi=[\psi_{1},\ldots,\psi_{n}]\in\mathbb{R}^{{n_{\alpha}}\times{n_{\alpha}}}$ is the matrix of orthonormal eigenvectors of $\mathscr{L}_{\mathcal{G}_{\alpha}}$, and $\Lambda=\text{diag}(\lambda_{1},\lambda_{2},\ldots,\lambda_{n_{\alpha}})$ is the diagonal matrix of the corresponding non-negative eigenvalues.
Note that $\{\psi_{i}\}_{i=1}^{n_{\alpha}}$ is an orthonormal basis for functions on the graph, since $\mathscr{L}_{\mathcal{G}_{\alpha}}$ is symmetric.
When the graph is constructed on points randomly sampled on a domain or a manifold, this discrete Laplacian is an approximation (in a suitable sense) to a continuous Laplacian in the limit as the number of samples goes to infinity, and its eigenvectors are approximations to the eigenfunctions of the continuous Laplacian, which are the natural generalization of Fourier modes from a torus or rectangular box to general domains and manifolds. Furthermore, the geometry of the domain (graph, in the discrete case, Riemannian manifold in the continuous case) is completely determined by the Laplacian, and therefore by its eigenvalues and eigenvectors.

We define the Graph Fourier Transform (GFT) as follows: for a function $a\in\mathbb{R}^{n_{\alpha}}$ on the vertices of $\mathcal{G}_{\alpha}$, the Graph Fourier Transform (GFT), denoted by $\mathcal{F}$, rewrites $a$ in the basis of eigenvectors of the graph Laplacian $\mathscr{L}_{\mathcal{G}_{\alpha}}$, yielding a function $\hat{a}$ of the eigenvalues $\lambda_{1},\dots,\lambda_{n_{\alpha}}$, defined as

$$ $\hat{a}:=\mathcal{F}a:=(\langle\psi_{i},a\rangle)_{i=1,\dots,{n_{\alpha}}}= \Psi^{T}a\,,$ (5) $$

In the reverse process, the inverse Graph Fourier Transform (IGFT), represented by $\mathcal{F}^{-1}$, reconstructs the function $a$ from its spectral representation $\hat{a}$, through the equation

$$ $a=\mathcal{F}^{-1}\hat{a}:=\sum_{i=1}^{n_{\alpha}}\hat{a}(\lambda_{i})\psi_{i} =\Psi\hat{a}\,.$ (6) $$

It should be noted that computing the $k_{\text{max}}$ lowest eigenvalues for an undirected graph is generally inexpensive with a complexity of $\mathcal{O}(k_{\text{max}}^{2}n_{\alpha}j)$ for a graph with $n_{\alpha}$ vertices, $j$ neighbors per vertex. Eigenpairs are of course computed once per domain/diffusivity parameter, rather than being recomputed at every instance during network training.

Graph diffusion is analogous to the heat equation and describes the diffusive movement of a function across a graph.
The solution to the graph diffusion equation can be expressed in the basis of the spectral decomposition of the graph Laplacian $\mathscr{L}_{\mathcal{G}_{\alpha}}$ as

$$ $\frac{\partial a(t)}{\partial t}=c\mathscr{L}_{\mathcal{G}_{\alpha}}a(t),\quad a (t)=\sum_{i=1}^{n_{\alpha}}\mathcal{F}(a(0))_{i}e^{-c\lambda_{i}t}\psi_{i}$ (7) $$

where $c$ is the diffusivity constant, and $a(t)\in\mathbb{R}^{n_{\alpha}}$ is the function at time $t$, and the initial condition is the function $a(0)=(a_{i}(0))_{i=1}^{n_{\alpha}}$. The solution can be truncated to the $k_{\text{max}}$ lowest modes, where $k_{\text{max}}\leq n_{\alpha}$, for an approximated solution. It should be noted that the form of Eq. (7) is not dependent upon a spatial coordinate but only on the spectral decomposition of $\mathscr{L}_{\mathcal{G}_{\alpha}}$.

By establishing these preliminaries, we set the stage for the development of G-FuNK, which leverages the spectral properties (precomputed lowest eigenpairs) of the graph Laplacian to efficiently learn and generalize dynamic processes across different domains.

### 2.2 The G-FuNK layers and network

In this section, we describe the neural operator framework, a variation of the FNO of Li et al. (), equipped with our novel Graph Fourier Neural Kernel (G-FuNK) layers to learn solution generators to time-dependent PDEs across multiple domains and multiple anisotropic diffusion tensor fields.

Recall that the domain – or, rather, a discretization thereof – is given, and so is the diffusion field $\mathbf{K}$.
We exploit the knowledge of the former to construct a graph $\mathcal{G}_{\alpha}$ which discretizes $\Omega_{\alpha}$ into $n_{\alpha}$ nodes $\mathcal{V}=\{x_{i}\}_{i=1}^{n_{\alpha}}$, and knowledge of the latter to weight the edges in order to approximate the second-order term $\nabla\cdot(\mathbf{K}(\mathbf{x})\nabla)$ in the PDE. By choosing, on each edge $\mathcal{E}_{ij}$,

$$ $w_{ij}^{-1}:=\frac{1}{2}(x_{j}-x_{i})^{T}\cdot(\mathbf{K}(x_{i})^{-1}+\mathbf{ K}(x_{j})^{-1})\cdot(x_{j}-x_{i})\,,$ (8) $$

we weight the edges at $x_{i}$ proportionally to the direction of the diffusion coefficient $\mathbf{K}(x_{i})$ at $x_{i}$, and to the squared inverse of the distance. The average of $\mathbf{K}(x_{i})^{-1}$ and $\mathbf{K}(x_{j})^{-1}$ is used to obtain a symmetric graph.

The neural operator $\mathcal{N}_{\theta}$, equipped with several G-FuNK layers, approximates the solution generator of the PDE by applying a series of transformations, which can be formally expressed as the composition of different functions representing the network layers. The network begins with the following function transformation:

$$ $\displaystyle k_{0}(x)$ $\displaystyle=\mathcal{P}(a(x);\theta_{P}),\quad$ (9) $$

where $\mathcal{P}:\mathbb{R}^{n_{\alpha}\times d_{a}}\rightarrow\mathbb{R}^{n_{
\alpha}\times d_{p}}$ is the initial lifting operation to a higher dimensional feature space.

For each G-FuNK layer $n=1,\ldots,N$, we have, following the structure diagram in Figure 1 from left to right:

$$ $\displaystyle\hat{k}_{n}(\lambda)$ $\displaystyle=\mathcal{F}(k_{n}(x);\Psi^{T})(\lambda),\quad$ (10) $\displaystyle\hat{\ell}_{n}(\lambda)$ $\displaystyle=\mathcal{L}_{n}(\hat{k}_{n}(\lambda);B)(\lambda)\,,$ $\displaystyle\hat{r}_{n}(\lambda)$ $\displaystyle=\mathcal{R}_{n}(\hat{\ell}_{n}(\lambda);\theta_{\mathcal{R}_{n}} )(\lambda)\,,$ $\displaystyle\hat{f}_{n}(x)$ $\displaystyle=\mathcal{F}^{-1}(\hat{r}_{n}(\lambda);\Psi)(x)\,,\quad$ $\displaystyle\omega_{n}(x)$ $\displaystyle=\mathcal{W}_{n}(k_{n}(x);\theta_{\mathcal{W}_{n}})(x)\,,$ $\displaystyle z_{n}(x)$ $\displaystyle=\sigma(\hat{f}_{n}(x)+\omega_{n}(x))(x)=k_{n+1}(x),$ $$

where $\mathcal{F}:\mathbb{R}^{n_{\alpha}\times d_{p}}\rightarrow\mathbb{R}^{k_{\text
{max}}\times d_{p}}$ represents the Graph Fourier Transform that projects the function into the $k_{\text{max}}$ lowest frequencies in the spectral domain,
$\mathcal{L}_{n}:\mathbb{R}^{k_{\text{max}}\times d_{p}}\rightarrow\mathbb{R}^{
k_{\text{max}}\times d_{p}\times p}$ for $n=1,\ldots,N$ is a linear transformation of $\hat{k}_{n}$ together with a set of the eigenvalues raised to $p$ powers (i.e. $B$ is a matrix $\mathbb{R}^{k_{\text{max}}\times p}$ with columns of $(\lambda^{0}_{1},\dots,\lambda^{0}_{k_{\text{max}}}),(\lambda^{1}_{1},\dots,
\lambda^{1}_{k_{\text{max}}}),\dots,(\lambda^{p-1}_{1},\dots,\lambda^{p-1}_{k_
{\text{max}}})$), $\mathcal{R}_{n}:\mathbb{R}^{k_{\text{max}}\times d_{p}\times p}\rightarrow
\mathbb{R}^{k_{\text{max}}\times d_{p^{\prime}}}$ for $n=1,\ldots,N$ are the parameterized linear transformations in the spectral domain corresponding to the $n$-th G-FuNK layer which can be diagional, tri-diagional, full matrices, etc., $\mathcal{F}^{-1}:\mathbb{R}^{k_{\text{max}}\times d_{p^{\prime}}}\rightarrow
\mathbb{R}^{n_{\alpha}\times d_{p^{\prime}}}$ is the inverse Graph Fourier Transform that maps the transformed spectral function back to the graph domain. $\mathcal{W}_{n}:\mathbb{R}^{n_{\alpha}\times d_{p}}\rightarrow\mathbb{R}^{n_{
\alpha}\times d_{p^{\prime}}}$ for $n=1,\ldots,N$ are additional linear mappings applied after transforming back to the graph domain in each G-FuNK layer, $\sigma:\mathbb{R}^{n_{\alpha}\times d_{p^{\prime}}}\rightarrow\mathbb{R}^{n_{
\alpha}\times d_{p^{\prime}}}$ an activation function applied in each G-FuNK layer.

After processing through all $N$ G-FuNK layers, the final output is obtained by the projection layer:

$$ $\displaystyle y(x)=\mathcal{Q}(z_{N}(x);\theta_{\mathcal{Q}})$ (11) $$

where $\mathcal{Q}:\mathbb{R}^{n_{\alpha}\times d_{p^{\prime}}}\rightarrow\mathbb{R}^
{n_{\alpha}\times d_{u}}$ projects from the last layer’s output onto the target space.
All together, the neural operator with $N$ G-FuNK layers can be expressed as the following composition of the afforementioned transformations to approximate $\partial_{t}u(\mathbf{x},t)$:

$$ $\mathcal{N}_{\theta}(a(x)):=\mathcal{Q}\circ\underbrace{\sigma(\mathcal{W}_{N} +\mathcal{F}^{-1}\circ\mathcal{R}_{N}\circ\mathcal{L}_{N}\circ\mathcal{F})}_{ \text{G-FuNK Layer $N$}}\circ\cdots\circ\underbrace{\sigma(\mathcal{W}_{1}+ \mathcal{F}^{-1}\circ\mathcal{R}_{1}\circ\mathcal{L}_{1}\circ\mathcal{F})}_{ \text{G-FuNK Layer $1$}}\circ\mathcal{P}(a(x);\theta),$ (12) $$

We denote $\theta=\{\theta_{\mathcal{P}_{n}},\theta_{\mathcal{R}_{n}},\theta_{\mathcal{W}
_{n}},\theta_{\mathcal{Q}_{n}}\}_{n=1}^{N}$ as a collection of all learnable parameters in the network that are optimized during the training phase to minimize the discrepancy between the neural operator’s output and the known numerically computed trajectories. Figure 1 illustrates a summary of the network architecture.

Figure: Figure 1: Top: Network structure of the Neural Operator. The input is passed as $a$. (1) It is lifted to a higher dimensional space via a lifting operator $\mathcal{P}$. (2) Multiple G-FuNK layers are then applied. (3) The network projects to the output dimensional space via $\mathcal{Q}$. Bottom: Single G-FuNK layer structure. Starting from input $k$, (1) $\mathcal{F}$ denotes the GFT of the input which is then reduced to the $k_{\text{max}}$-lowest modes. (2) The function is expanded onto a set of the eigenvalues raised to different powers in $\mathcal{L}$ (3) A linear transform is applied with learnable parameters in $\mathcal{R}$. (4) $\mathcal{F}^{-1}$ represents the inverse GFT, mapping back to the original domain. (5) The input function $k$ undergoes a linear mapping in $\mathcal{W}$. (6) The outputs of the top and bottom branches are combined and then passed through an activation function $\sigma$. The output of the G-FuNK layer is $z$ and is operated on by the next G-FuNK layer.
Refer to caption: x1.png

Given discrete observations $\{u({x_{i}},t_{\ell})\}_{\ell=1}^{L}$ of the state variable at a location ${x_{i}\in\mathcal{G}_{\alpha}}$ at different time steps $0=t_{1}<\dots<t_{L}$, we can integrate the approximation over time using an ODE solver, $\mathcal{I}$:

$$ $\displaystyle u(x_{i},t_{\ell}+\Delta t)=\mathcal{I}\left(\mathcal{N}_{\theta} (u(x_{i},t_{\ell},\Psi,\Lambda)\right),$ (13) $$

where $\Delta t$ is the time step size, and $\mathcal{I}$ is a neural ODE numerical algorithm (we use the one proposed by Chen et al. ()). This integrator uses the adjoint method during the gradient computation in backpropagation (since the loss function involves predicting the solution at multiple future time points), which makes training computationally demanding.

In the G-FuNK framework, the eigenvectors of the Laplacian perform the same action as the Fourier transform in the FNO framework, with the important difference that these eigenvectors are adapted to the diffusion coefficient, capturing the parametric dependence of the highest order differential term in the PDE. We suggest using the true eigenfuctions of the shape which embed global structural and physical properties as corroborated by numerical and statistical methods which represent PDE solutions in terms of known basis functions that contain information about the solution structure Bhattacharya et al. (); Nagy (); Almroth et al. (). Furthermore, in FNOs, the eigenvalues are not needed as there is no domain change; here by incorporating the eigenvalues we capture both the dependency on $\mathbf{K}$ and on the domain $\mathcal{G}_{\alpha}$.

Finally, we comment that $\mathcal{L}_{n}$ and $\mathcal{R}_{n}$ allow G-FuNK to approximate spectral multipliers, while $\mathcal{W}_{n}$ allows for the approximation of spatial pointwise multipliers. At a high level, this is particularly well-suited for reaction-diffusion equations, where the diffusion component is expected to be easy to learn through spectral multiplier, and the nonlinear reaction terms act pointwise and can be expected to be captured by the spatial multipliers. Of course, such PDEs have complex interactions between the spatial and spectral domain, so this intuition is very much qualitative.

##### Mesh-Independence of G-FuNK

Although the proposed framework involves mesh-dependent steps, such as forming a weighted graph on the discretized domain, the method maintains a high degree of mesh-independence. More precisely, as the number of points increases, the Laplacian we construct approaches its limit, allowing our framework to effectively handle finer discretizations. Moreover, there are established methods for interpolating eigenfunctions to an underlying continuous domain Coifman et al. (), which further supports the transition from discrete to continuous representations. This characteristic is crucial, as it sheds light on the ability of our method to transition from discrete approximations to continuous domains as the mesh size approaches zero. Our framework’s capability to interpolate and generalize across varying mesh sizes ensures that it remains robust and accurate in capturing the underlying dynamics of PDEs in the limit of fine discretizations.

## 3 Numerical Experiments

In this section, we present several numerical experiments demonstrating the ability of the G-FuNK operator learning framework to accurately predict the solution generator for PDEs of the form in (2), towards the goal of accelerating precision medicine in cardiac electrophysiology. Results for all examples and comparisons to baseline models are summarized in Table 1.

Heat Equation.
The introduction of anisotropic diffusion is critical in simulating phenomena such as the propagation of heat in materials with fibrous structures, where the diffusive properties vary along different directions.
We consider the classical heat equation with very strong anisotropic diffusion on a 2D unit square domain with no-flux Neumann boundary conditions, with $\mathbf{K}$ as a diffusion tensor representing the anisotropy:

$$ $\partial_{t}u=\nabla\cdot(\mathbf{K}(\textbf{x})\nabla u)\,,\quad\mathbf{x}\in \Omega_{\alpha},t>0,\quad\partial_{\mathbf{n}(x)}u(\mathbf{x})=0\>\,\forall\, \mathbf{x}\,\in\,\partial\Omega_{\alpha}\,$ (14) $$

In our experiments, we consider a diffusion tensor with an anisotropy ratio of 9:1, given by:

$$ $\mathbf{K}(\textbf{x})=\left[\begin{array}[]{cc}\hat{F}^{(1)}(\textbf{x})&\hat {F}^{(2)}(\textbf{x})\\ \end{array}\right]\begin{bmatrix}9&0\\ 0&1\end{bmatrix}\left[\begin{array}[]{cc}\hat{F}^{(1)}(\textbf{x})&\hat{F}^{(2 )}(\textbf{x})\\ \end{array}\right]^{T},$ (15) $$

where $\hat{F}^{(1)}(\textbf{x})$ and $\hat{F}^{(2)}(\textbf{x})$ are orthonormal vectors fields that define the longitudinal and transverse directions of diffusion and are described in Appendix section A.1.1, equation (17). The longitudinal direction of the anisotropic diffusion were defined as a superposition of geometric and linear functions, with random parameters sampled independently for each trajectory in both training and test data.

Figure: Figure 2: Heat Equation: G-FuNK vs Target. A comparison between the target and the G-FuNK predicted heat equation. On the left is a stream plot example of the primary direction of the diffusion, described as $\hat{F}^{(1)}$ in (17), which was unique for the test trajectory shown on the right. The prediction of G-FuNK is on the bottom and the target is above for three different time points after the initial condition at $t_{0}$.
Refer to caption: x2.png

For this case, the eigenvalues in $B$ of (10) were raised to only the $0^{\text{th}}$ power, and fiber information was not used as an input to the network. Therefore, learning the influence of anisotropic diffusion was done using only the influence of the eigenvectors. For this case, $\theta_{R}$ was only a diagonal matrix of learnable parameters.
Figure 2 contains a depiction of a diffusion field $\hat{F}^{(1)}$. Additional information about the network can be found in A.1.

##### 2D Nonlinear Reaction Diffusion.

We consider a reaction diffusion system in a 2D domain with no-diffusive-flux Neumann boundary conditions. This system of equations is very stiff and involve multiple non-linear equations to describe the propagation of action potentials in the form of a bioelectrical wave Courtemanche et al. (). The solutions of these equations typically have the range between $-85$ and $20\>\text{mV}$, with the wavefront having a rate of change $>100\>\frac{\text{mV}}{\text{ms}}$, making them steep and close to discontinuous.
For these systems, an external stimulus current is applied to start the wave propagation from a random point within the domain for each trajectory.

In terms of the transmembrane potential $u$, the system is of the form

$$ $\partial_{t}u=\nabla\cdot(\mathbf{K}(\textbf{x})\nabla u)+\sum_{s}J_{s}(u, \mathbf{v}),\quad\frac{d\mathbf{v}}{dt}=\Upsilon(u,\mathbf{v}),\quad\mathbf{K} (\mathbf{x})\nabla u\cdot\mathbf{n}(\mathbf{x})=0\>\forall\>\mathbf{x}\>\in\> \partial\>\Omega_{\alpha}\,,$ (16) $$

where $\mathbf{x}\in\Omega_{\alpha},t>0$, $\sum_{s}J_{s}$ is the sum of 12 ionic and external stimulus currents described by non-linear functions Courtemanche et al. (). In this case, no diffusive flux Neumann boundary conditions were used.

In this example, each trajectory had a unique rectangular domain, $\Omega_{\alpha^{(m)}}$, with edge length drawn independently from 2 uniform distributions on $[15,30]$. The diffusive tensor $\mathbf{K}(\textbf{x})$ was set at 5 to 1 in the x direction across each domain. Results are shown in Figure 3. Additional information about the network can be found in A.2.

Figure: Figure 3: Reaction Diffusion on Random Rectangle: G-FuNK vs Target. A and B are two comparisons of G-FuNK for 2 different test geometries. Both panels use the colorbar on the right which is min-max scaled from -85 to 20 mV and $\Delta t$ is 10 milliseconds.
Refer to caption: x3.png

Cardiac Electrophysiology.
In the field of personalized cardiac electrophysiology (EP), the ability to describe the electrical propagation in the Left Atria (LA) is of the utmost importance to provide better healthcare to patients inflicted by electrical abnormalities. Using cardiac computed tomography (CT) scans from 25 atrial fibrillation (AF) patients, triangulated meshes of the LA surface were constructed from image segmentations. The mitral valve and four pulmonary veins were removed via visual inspection resulting in a topology with five holes. The fiber orientation for each mesh was mapped using an atlas-based procedure Roney et al. (); Ali et al. (); Roney et al. () and produced fibers similar to the examples shown in Figure 4. The process of manually annotating the geometry to define the fiber fields adds inherent noise from one geometry to another.

Finite element simulations were computed using the openCARP software Plank et al. () in which the Courtemanche electrophysiological ionic model Courtemanche et al. () was solved with no-flux Neumann boundary conditions to describe the electrical wave propagation in each domain. The Courtemanche equations follow the same form as (16). The diffusive ratio was set as $5$ to $1$ in the fiber direction with a base diffusive value of $0.625\>\frac{\text{mS}}{\text{mm}}$. Each simulation incorporated a stimulus current, represented as a delta spike, chosen at a random location within the domain and the initial condition was chosen at 10 ms after this stimulus was applied. For this work, the aim was to start from an initial condition $u_{0}$ and compute the next time step repeatedly. This was done using 24 domains, G-FuNK was trained to infer the solution from $0\text{ms}$ to $90\text{ms}$ which focused on the steep wavefront, and an additional, unseen geometry was used as an out-of-training test set as shown in
Figure 5. Information about the network can be found in A.3.

Figure: Figure 4: Example of a patient-specific fiber field providing anisotropic information for the PDE on the 3D left atrial chamber surface.
Refer to caption: x4.png

Figure: Figure 5: Comparing target and G-FuNK predictions on an out-of-training test geometry. A Left atrial posterior wall view at $2\Delta t$ (20 ms) after initial condition. B Left atrial anterior wall view at $4\Delta t$ (40 ms) after initial condition.
Refer to caption: x5.png

##### Results

The performance of G-FuNK and other models is summarized in Table 1 below.

**Table 1: Performance of G-FuNK, FNO (or Geo-FNO) and a GNN equipped with a neural ODE on the different problems presented in 3. The performance is on a test set with new initial conditions and, where applicable, new domains and diffusive fields not seen during training.**
| Problem | Geometries | G-FuNK<br>Parameters<br>Rel $\ell_{2}$ | FNO<br>Parameters<br>Rel $\ell_{2}$ | GNN<br>Parameters<br>Rel $\ell_{2}$ |
| --- | --- | --- | --- | --- |
| Heat<br>Equation | 2D Unit Square<br>Varying Diffusion<br>Training Trajectories: $M=500$<br>(GNN: $M=1000$) | 197465<br>0.0302 | FNO<br>2017157<br>0.0134 | With Edge Weights:<br>12226<br>0.1875<br>Without Edge Weights:<br>35026<br>0.2428 |
| Reaction<br>Diffusion | 2D Random Rectangles<br>No Varying Diffusion<br>Training Trajectories: $M=1000$<br>Test Trajectories Rotated $90^{\circ}$ | 134965<br>0.1189<br>Invariant<br>0.1189 | Geo-FNO<br>2515307<br>0.1162<br>Not Invariant<br>0.5681 | With Edge Weights:<br>20751<br>0.3345<br>Not Invariant<br>0.3663 |
| Cardiac EP | Multiple 3D Patient Geometries<br>Varying Diffusion<br>Training Trajectories: $M=1000$ | 283040<br>0.1642 | Not<br>Applicable | With Edge Weights:<br>248851<br>0.4034 |

##### Discussion

Altogether, G-FuNK was adaptable to all the examples discussed above with low relative error. The results for the various tasks are listed in Table 1.
G-FuNK performed well on all examples with a fairly small amount of learnable parameters compared to FNOs and GeoFNOs Li et al. (; ).

Message-passing GNNs aggregate information at the nodes, and as expected are unable to accurately handle directionally dependent information (varying diffusive fields). We notice that when we incorporate edge weights, as defined in equation (8), into the GNN model for the Heat equation example, the results improve and the model requires much fewer parameters. Therefore, for all comparisons following, we incorporated the edge weights as an input into the different networks.

Prior works such as FNO and GeoFNO, with the code provided, do not present results or compatibility for learning full trajectories of time-evolving processes, where errors can accumulate due to spatial misalignment of predictions Li et al. (; ). To compare with FNO and GeoFNO, we modified the approaches to be compatible with the same neural ODE solver we use for time-dependent predictions with G-FuNK. For the heat example, it should be noted that G-FuNK was able to learn the effects of anisotropic diffusion from only the eigenvectors. The other networks were given the primary diffusive vector as inputs at each sample point which has less meaning outside of planar geometries. In the 2D reaction-diffusion problem with random rectangles, both GeoFNO and G-FuNK perform similarly on the test set, as expected. In regular rectangular domains, the graph Fourier and fast Fourier bases converge, implying that FNO and GeoFNO are subsets of G-FuNK equipped neural operators, corroborated by the similar performance. For training in this example, we still considered anisotropy, but with all the fiber fields strictly pointing in the horizontal axis. We show that G-FuNK and GeoFNO perform quite similarly on test domains with the same directionality presented in the fiber fields as the training samples. However, when the test set domains and fiber fields were rotated by 90 degrees, GeoFNO’s performance significantly degraded, with a relative $\ell_{2}$ error of 0.5681, as shown in Table 1 and Figure 6. In contrast, G-FuNK, trained using only the initial condition and graph eigenpairs, remained rotation-invariant to the unseen fiber orientation.

Figure: Figure 6: Geo-FNO Performance: Target vs Prediction for the random rectangle case. A Geo-FNO performance on the original test set. B Geo-FNO performance on the same original test set rotated $90$ degrees. For A and B, the colorbar is min-max scaled from $-85$ to $20$ mV and $\Delta t$ is $10$ milliseconds.
Refer to caption: x6.png

For the cardiac EP example, we do not show results with GeoFNO since the topologically complex 3D geometries of the human left atrium with five holes cannot be mapped diffeomorphically to a cube or torus where the Fourier transform is available. For this example, our G-FuNK equipped operator learning framework is able to provide respectable predictions, where a majority of the relative $\ell_{2}$ error is due to a temporal misalignment of the prediction and the ground truth which can be specifically attributed to a small lag of about $1.62$ ms in the wavefront. This shift resulted in a higher $\ell_{2}$ error but the cross-correlation is $0.941$ implying that the prediction is correctly shaped with a small lag.

Our method predicts entire trajectories in under 1 second, significantly outperforming traditional numerical methods. For example, cardiac EP simulations typically take at least 15 minutes on 12 CPU cores for one given set of initial conditions. To make quantitatively informed clinical predictions, one must perform comprehensive parameter sweeps to identify optimal treatment strategies for a given patient, which could take hours to days with the finite element approach. These trajectories are predicted by the G-FuNK equipped operator learning framework in less than 1 second. The inclusion of more geometries in the 3D cardiac EP examples is of course expected to lead to lower error. Additionally, small changes in the eigenvalues across domains can lead to mismatches in the order of the eigenvalues between geometries which could be a source in the reported error: this can be avoided with an eigenvector-matching procedure, which will be the subject of future investigations.

## 4 Conclusion

Data-driven methods for solving PDEs rely heavily on the quality and quantity of the data provided. In this work, we introduce a new family of neural operators based on our novel Graph Fourier Neural Kernel (G-FuNK) layers, which integrate GNNs with a variant of FNOs. We show that our framework achieves high accuracy on pedagogical examples, including the anisotropic 2D heat and reaction-diffusion equations. Even in the relatively data-sparse context of cardiac electrophysiology explored here, G-FuNK demonstrates respectable accuracy. Overall, we highlight the method’s expressiveness and versatility by learning PDE solution generators in the spectral domain of graphs, providing a natural way to embed geometric and directionally dependent information about anisotropic domains. Our G-FuNK-equipped operator learning framework effectively predicts the temporal dynamics of complex systems from a single initial condition, handling unstructured data, varying anisotropic diffusion tensor fields, highly nonlinear reaction terms, and different domains, including topologically complex 3D left atrial geometries.

## Acknowledgements:

Z.A. was supported by grant n. 24PRE1196125- American Heart Association (AHA) predoctoral fellowship and grant n. T32 HL007024 Cardiovascular Epidemiology Training Grant from the National Institute of Health (NIH) National Heart, Lung, and Blood Institute. N.T. was supported by National Institute of Health (NIH) grant n. R01HL166759.

## Author Contribution:

SEL: Conceptualization, Methodology, Software Implementation, Simulation, Formal Analysis, Data Curation, Visualization, Writing (original draft)
ZA: Methodology, Software Implementation, Simulation, Formal analysis, Visualization, Writing (original draft)
SYA: Data Curation (Fiber Generation)
CY: Data Curation (Left Atrial Segmentation)
DP: Conceptualization, Methodology, Software Implementation
AY: Data Curation (Left Atrial Mesh Generation, Universal Atrial Coordinates Generation)
YL: Data Curation (Left Atrial Mesh Generation, Universal Atrial Coordinates Generation)
NT: Conceptualization, Formal Analysis, Data Acquisition, Resources, Supervision, Project Administration, Funding Acquisition
MM: Conceptualization, Methodology, Formal Analysis, Supervision, Project Administration, Writing (review)

## Appendix A Appendix

### A.1 Heat Equation

#### A.1.1 Heat Data Generation

For this section of work, the heat equation was solved using the NDSolve function from the Mathematica software on the 2D unit square where the initial condition, $u_{0}$, was unique for each trajectory and chosen as a superposition of three binomial distributions where each center, deviation, and covariance were picked from uniform distributions with ranges of $[0,1]$, $[0.1,0.2]$ and $[-0.1,0.1]$, respectively. The primary direction of diffusion was defined as $F^{1}$ in the equation below.

$$ $F^{(1)}(x,y)=\begin{bmatrix}\sin(a_{1}(\frac{x+a_{2}}{2\pi}))+\cos(a_{3}(\frac {y+a_{4}}{2\pi}))+a_{5}x+a_{6}+a_{7}y+a_{8}\\ \sin(b_{1}(\frac{x+b_{2}}{2\pi}))+\cos(b_{3}(\frac{y+b_{4}}{2\pi}))+b_{5}x+b_{ 6}+b_{7}y+b_{8}\end{bmatrix}$ (17) $$

The parameters $a_{1}$ to $a_{8}$ and $b_{1}$ to $b_{8}$ were derived from the following uniform distributions and were unique for each trajectory in the training and test data sets:

| Uniform Distribution | Parameters |
| --- | --- |
| [0,2] | $a_{1}$, $b_{1}$, $a_{3}$, $b_{3}$ |
| [-1,1] | $a_{2}$, $b_{2}$, $a_{4}$, $b_{4}$, $a_{6}$, $b_{6}$, $a_{8}$, $b_{8}$ |
| [-2,2] | $a_{5}$, $b_{5}$, $a_{7}$, $b_{7}$ |

This allowed for vector fields like the ones shown in the Figure 7 below.

Figure: Figure 7: Three different and unique examples of the primary diffusion vector fields used for the Heat equation example in 3. Each example represents a diffusive field that was randomly determined for each trajectory
Refer to caption: x7.png

After all simulations were completed, a perturbed Delaunay triangulation grid was constructed and linear interpolation was used to generate the trajectories used for training at discretized points within the domain with a max cell measure of $0.0002$ and from 0 0 to $20$ms at every $1$ ms intervals. In addition, all other data, like diffusion and the vector denoting the primary direction of diffusion orientation, were calculated at each discretized point for future uses. Training and test data sets consisted of $1000$ and $100$ trajectories, respectively.

#### A.1.2 Network Description for the Anisotropic Heat Example

For each trajectory, an undirected Laplacian matrix was created such that each point was connected to its $30$ nearest neighbors and weighted using the heat kernel method described in the main text. The graph Laplacian was constructed using the NetworkX software and spectral decomposition was done using the Scipy software’s eigsh function.

The network was trained and optimized using the PyTorch software with an $\ell_{2}$ loss. The network is described as so, a. Linear mapping from three inputs $(u,x,y)$ to $200$ width. b. three G-FuNK layers each with a width of $200$ and using $50$ modes and using a GELU activation function. c. Two Linear layers going from a width of $200$ to $32$, with a GELU activation, and then $32$ to $1$. The network was integrated using the NeuralODE software with a forward Euler integrator with a time step of $0.25$ ms. G-FuNK was optimized with an $\ell_{2}$ loss function and an Adam optimizer with an initial learning rate of $5\cdot 10^{-4}$ with a step decayed learning rate which halved every $100$ epochs. G-FuNK utilized early stopping and was stopped when the validation loss did not improve for over $50$ epochs.

For the comparison models, FNO was trained via the code found on their [GitRepo](https://github.com/neuraloperator). The 2D FNO network is described as so a. A linear layer lifting the solution, primary diffusion vector and positional coordinates from a width of $5$ to $50$, b. next, $4$ FNO layers using $10$ modes for both directions with a width of $50$ and all used a GeLU activation function and finally c. a projection layer consisting of $2$ linear layers with the first expanding to a width of $128$ and the second to the output size of $1$. For FNO, solutions were linearly interpolated onto a 2D unit grid discretized by $32$ by $32$ points. For the GNN, a message-passing neural network was used similar to Iakovlev et al. () which consisted of $2$ networks. For the weighted edge case, the network is described as so a. A message network with an input width of 6 where solution at node $i$, solution difference between node $i$ and neighbors $j$, relative position between node $i$ and neighbors $j$ and difference in primary diffusion vectors at node $i$ and neighbors $j$ were used as inputs with a hyperbolic tangent activation function. There were no hidden layers used and next was an output linear layer of width $75$. b. An aggregation network consisting of an input linear layer with a width of $78$ to include the message, the solution at node $i$, and the primary diffusion vector at node $i$, and a hyperbolic tangent function. This network has no hidden layers and next was an output linear layer which projected to a width of $1$. The GNN used the $30$ nearest neighbors as for connecting nodes with edges and the weights of these edges were calculated using the heat kernel described in 8.
The GNN without edge weights is described as so a. A message network with an input width of 6 where solution at node $i$, solution difference between node $i$ and neighbors $j$, relative position between node $i$ and neighbors $j$ and difference in primary diffusion vectors at node $i$ and neighbors $j$ were used as inputs. There were $3$ hidden layers of width $75$ with hyperbolic tangent activation functions and an output linear layer of width $75$. b. An aggregation network consisting of an input width of $78$ to include the message, the solution at node $i$ and the primary diffusion vector at node $i$. This network has 3 hidden layers of width $75$ with hyperbolic tangent activation functions and an output linear layer of width $1$. Both the FNO and GNN networks were integrated using the NeuralODE software with a forward Euler integrator with a time step of $1$ ms. Both used the same training protocol as G-FuNK.

### A.2 2D Reaction Diffusion

#### A.2.1 Data Generation - Random Rectangle Example

For this section of work, the 2D Reaction Diffusion equations were solved using the openCARP software which used a finite element solver on a random rectangle to solve the Courtemanche equations. The solver used a Crank-Nicolson method with a time step of $0.2$ ms and an average edge length of $0.2$ mm. Each edge of the rectangle was drawn from an independent uniform distribution of $[15,30]$ cm. The primary direction of diffusion was defined in the $x$-direction with a $5$ to $1$ anisotropic ratio and a base diffusion value of 0.0625 $\nicefrac{{mS}}{{mm}}$.

After all simulations were completed, a down-sampled point cloud was calculated from the finite element mesh and nearest neighbor interpolation was used to generate the trajectories used for training at discretized points within the domain with a resolution of one point per $0.25$ mm voxel (via open3d voxel_down_sample_and_trace), and from 0 0 to $100$ ms at every $10$ ms intervals. In addition, all other data like diffusion, vector orientation, etc. was calculated at each discretized point for future uses. Training and test data sets consisted of $1000$ and $100$ trajectories, respectively.

#### A.2.2 Network Description for Random Rectangle Example

For each trajectory, an undirected graph Laplacian matrix was created such that each point was connected to its $30$ nearest neighbors and weighted using the heat kernel method described in the main text. The graph Laplacian was constructed using the NetworkX software and spectral decomposition was done using the Scipy software’s eigsh function.

The network was trained and optimized using the PyTorch software with a loss function with the form $||\hat{u}-u||_{2}+5\cdot||\nabla\hat{u}-\nabla u||_{2}$. This loss was used to minimize the error at the wavefront where the difference in the gradient is typically the largest. G-FuNK was optimized with the Adam optimizer with a step decay learning rate where the learning rate halved every 100 epochs starting at $5\cdot 10^{-4}$. For this case, the eigenvalues used in $B$ in (10) were raised to powers 0 0, $1$, and $2$ for these cases. For both cases, $\theta_{R}$ was a full matrix of learnable parameters.

For the random rectangle example, the network is described as so, a. Linear mapping from $1$ input ($u$) to 200 width. b. Three G-FuNK layers each with a width of $200$ and using $50$ modes and each using a GELU activation function. c. Two Linear layers going from a width of $200$ to $32$, with a GELU activation, and then $32$ to $1$. The network was integrated using the NeuralODE software with a forward Euler integrator with a time step of $1$ ms.

For the comparison models, Geo-FNO was trained via the code found on their [GitRepo](https://github.com/neuraloperator). The 2D Geo-FNO network is described as so a. A linear layer lifting the solution, fiber vector and positional coordinates from a width of $5$ to $50$, b. next, $4$ FNO layers using $10$ modes for both directions with a width of $50$ and all used a GeLU activation function and finally c. a projection layer consisting of two linear layers with the first expanding to a width of $128$ and the second to the output size of $1$. The latent space was defined as a unit square discretized in $32$ by $32$ points. For the GNN, a message-passing neural network was used similar to Iakovlev et al. () which consisted of two networks as so: a. A message network with an input width of $6$ where solution at node $i$, solution difference between node $i$ and neighbors $j$, relative position between node $i$ and neighbors $j$ and difference in primary diffusion vectors at node $i$ and neighbors $j$ were used as inputs. There were $3$ hidden layers of width $40$ with hyperbolic tangent activation functions and an output linear layer of width $40$. b. An aggregation network consisting of in input width of $43$ to include the message, the solution at node $i$ and the primary diffusion vector at node $i$. This network has $3$ hidden layers of width $40$ with hyperbolic tangent activation functions and an output linear layer of width 1. The GNN used the $30$ nearest neighbors as neighbors for each node and edges weights were calculated using the heat kernel described in 8. Both the Geo-FNO and GNN networks were integrated using the NeuralODE software with a forward Euler integrator with a time step of $1$ ms. Geo-FNO was trained using a $\ell_{2}$ loss and the GNN was trained using a relative $H_{1}$ loss function. Both used the same training protocol as G-FuNK.

### A.3 Cardiac Electrophysiology

#### A.3.1 Data Generation for the Cardiac Electrophysiology Examples

From cardiac CT scans of $25$ patients, the left atria were segmented and a surface mesh was extracted using the Seg3D software. Manual corrections were made using MeshMixer and average edge length was set to $0.25$ mm. Surfaces were then annotated for calculating universal atrial coordinates for mapping atrial fibers with the atlas based method. The Courtemanche equations Courtemanche et al. () were solved using the openCARP software which used a finite element solver which implemented a Crank-Nicolson scheme with a time step of $0.2$ ms.

The diffusion tensor $\mathbf{K}(\mathbf{x})$ was calculated using the following equations.

$$ $\mathbf{\hat{F}}(\mathbf{x})=\left[\begin{array}[]{ccc}\vrule&\vrule&\vrule\\ \hat{F}^{(1)}(\mathbf{x})&\hat{F}^{(2)}(\mathbf{x})&\hat{F}^{(3)}(\mathbf{x}) \\ \vrule&\vrule&\vrule\end{array}\right]$ (18) $$

$$ $\mathbf{K}(\mathbf{x})=\mathbf{\hat{F}}(\mathbf{x})\left[\begin{array}[]{ccc} \sigma^{(1)}&0&0\\ 0&\sigma^{(2)}&0\\ 0&0&\sigma^{(3)}\end{array}\right]\mathbf{\hat{F}}^{T}(\mathbf{x})$ (19) $$

Where $\hat{F}^{(1)}$ is the mapped fiber field from the atlas-based method described in Roney et al. (); Ali et al. (); Roney et al. (), $\hat{F}^{(1)}$, $\hat{F}^{(2)}$, and $\hat{F}^{(3)}$ are orthonormal and $\sigma^{(2)}$ and $\sigma^{(3)}$ are the base diffusive coefficient of $0.625$ $\nicefrac{{mS}}{{mm}}$ with $\sigma^{(1)}$ set as five times greater.

After all simulations were completed, a down-sampled point cloud was calculated from the finite element mesh and nearest neighbor interpolation was used to generate the trajectories used for training at discretized points within the domain with a resolution of one point per $0.4$ mm voxel (via open3d’s function voxel_down_sample_and_trace), and the trajectory was downsampled from every 1 ms to every 10 ms in the range 0 0 to $300$ ms. In addition, all other data like diffusion and the fibers vector orientation were calculated at each discretized point for future uses. For the multiple geometries case, training and test data sets consisted of $2400$ and $100$ trajectories, respectively where each geometry had $100$ trajectories and the test set consisted of all 100 trajectories of a hold-out geometry.

#### A.3.2 Network Description for Cardiac Electrophysiology Examples

For each trajectory, an undirected graph Laplacian matrix was created such that each point was connected to its $6$ nearest neighbors and weighted using the heat kernel method described in the main text. The graph Laplacian was constructed using the NetworkX software and spectral decomposition was done using the Scipy software’s eigsh function.

For this case, G-FuNK was trained and optimized using the PyTorch software with a relative $H_{1}$ loss function with the form $\frac{||\hat{u}-u||_{2}}{||u||_{2}}+\frac{||\nabla\hat{u}-\nabla u||_{2}}{||
\nabla u||_{2}}$. This loss was used to minimize the error at the wavefront where the difference in the gradient is typically the largest.
An Adam optimizer equipped with a step decay learning rate with an initial rate of $5\cdot 10^{-4}$ that halves every 100 epochs was used and was early stopped when the validation loss did not decrease after 50 epochs. The eigenvalues used in $B$ in (10) were raised to the powers $0,1,2$ in these cases. $\theta_{R}$ was a full matrix of learnable parameters.

For this case, the G-FuNK network is described as so, a. Linear mapping from $1$ input ($u$) to $300$ width. b. Three G-FuNK layers each with a width of 300 and using 25 modes and each using a GELU activation function. c. Two Linear layers going from a width of $300$ to $32$, with a GELU activation, and then $32$ to $1$. The network was integrated using the NeuralODE software with a forward Euler integrator with a time step of $1$ ms.

For the comparision GNN model, a message-passing neural network was used similar to Iakovlev et al. () which consisted of two networks as so a. A message network that operated on the edges with an input width of $8$. The input features consisted of the solution at node $i$, solution difference between node $i$ and neighbor $j$, relative position between node $i$ and neighbor $j$, and difference in primary diffusion vectors at node $i$ and neighbors $j$ were used as inputs. There were three hidden layers of width $175$ with hyperbolic tangent activation functions and an output linear layer of width $175$. b. An aggregation network consisting of an input width of $179$ to include the message, the solution at node $i$ and the primary diffusion vector at node $i$. This network has three hidden layers of width $175$ with hyperbolic tangent activation functions and an output linear layer of width $1$. The GNN used the $6$ nearest neighbors as neighbors for each node and edges weights were calculated using the heat kernel described in 8. The network was integrated using the NeuralODE software with a forward Euler integrator with a time step of $1$ ms. A relative $H_{1}$ loss function was used to optimize the model and it implemented the same training protocol as G-FuNK.

### A.4 Other Examples

#### A.4.1 2D Reaction Diffusion with Random Fiber Fields and Anisotropic Diffusion Ratio

##### Data Generation - Random Fiber Field and Ratio Example

For this section of work, the 2D Reaction Diffusion equations were solved using the openCARP software which used a finite element solver on a $20$ cm by $20$ cm square to solve the Courtemanche equations. The solver used a Crank-Nicolson method with a time step of $0.2$ ms and an average edge length of $0.2$ mm. The primary direction of diffusion was defined as $F^{1}$ in the equation below.

$$ $F^{(1)}(x,y)=\begin{bmatrix}\frac{a_{1}x}{L}+a_{2}+\cos(\frac{a_{3}\pi y}{L}+a _{4})\\ \frac{b_{1}x}{L}+b_{2}+\sin(\frac{b_{3}\pi y}{L}+b_{4})\end{bmatrix}$ (20) $$

The parameter $L$ was set to the edge length of the square, 20 cm. The parameters $a_{1}$ to $a_{4}$ and $b_{1}$ to $b_{4}$ were derived from the following uniform distributions and were unique for each trajectory in the training and test data sets:

| Uniform Distribution | Parameters |
| --- | --- |
| [.35,.65] | $a_{1}$, $b_{1}$, $a_{2}$, $b_{2}$, $a_{3}$, $b_{3}$, $a_{4}$, $b_{4}$ |

For this case, $\mathbf{K}$ is a diffusion tensor representing the anisotropy with the same form as (15) and (19) except the anisotropic ratio is a random parameter $\gamma$ drawn from a uniform distribution on $[4,6]$ with a base diffusion of 0.0625 $\nicefrac{{mS}}{{mm}}$. The diffusion was calculated using the equation 21 below.

$$ $\mathbf{K}(\textbf{x})=\left[\begin{array}[]{cc}\hat{F}^{(1)}(\textbf{x})&\hat {F}^{(2)}(\textbf{x})\\ \end{array}\right]K_{0}\begin{bmatrix}\gamma&0\\ 0&1\end{bmatrix}\left[\begin{array}[]{cc}\hat{F}^{(1)}(\textbf{x})&\hat{F}^{(2 )}(\textbf{x})\\ \end{array}\right]^{T},$ (21) $$

Where $K_{0}$ is the base diffusive coefficient of 0.0625 $\nicefrac{{mS}}{{mm}}$ and $\gamma$ is a random parameter drawn from a uniform distribution of [4,6].

After all simulations were completed, a down-sampled point cloud was calculated from the finite element mesh and nearest neighbor interpolation was used to generate the trajectories used for training at discretized points within the domain with a resolution of 250 and from 0 to 100ms at every 10 ms intervals. In addition, all other data like diffusion and vector orientation were calculated at each discretized point for future uses. Training and test data sets consisted of 1000 and 100 trajectories, respectively.

For the single geometry case, the network is described as so, 1. Linear mapping from 1 input (u) to 200 width. 2. Three G-FuNK layers each with a width of 200 and using 100 modes and each using a GELU activation function. 3. Two Linear layers going from a width of 200 to 32, with a GELU activation, and then 32 to 1. The network was integrated using the NeuralODE software with a forward Euler integrator with a time step of 1 ms.

The prediction vs. the ground truth for this example is shown in figure 8 below.

Figure: Figure 8: Target (top) vs Prediction (bottom) at different times after the initial condition $t_{o}$ for one domain in the random fiber and random anisotropic diffusion ratio case. The majority of the error is associated with the wavefront. $\Delta t$ is $10$ milliseconds and the colorbar is min-max scaled from -85 to 20 mV.
Refer to caption: x8.png

For this case, the network used is described as so, a. Linear mapping from three inputs $(u,x,y)$ to $50$ width. b. Four G-FuNK layers each with a width of $50$ and using $200$ modes and each using a GELU activation function. $c.$ Two Linear layers going from a width of $50$ to $32$, with a GELU activation, and then $32$ to $1$. The network was integrated using the NeuralODE software with a forward Euler integrator with a time step of $1$ ms.

#### A.4.2 Single Atrial Example

G-FuNK was trained on a single left atria using the same setup as the multiple EP example. The network was trained and tested on the same points and fiber fields with random initial conditions. The data consisted of training and test trajectories from time $[0,90]$ ms with $1000$ and $100$ trajectories, respectively.

Figure: Figure 9: A: Vector representation of the fiber orientation across a left atrium. The direction of the white arrows is the primary direction of diffusion which influences the wave propagation. B: Posterior view of the LA at $3\Delta t$ ($30$ ms) after $t_{0}$ for the single geometry case. C: Anterior view of an LA at $5\Delta t$ (50 ms) after $t_{0}$ for the single geometry case.
Refer to caption: x9.png

For this case, the network is described as so, a. Linear mapping from 1 input ($u$) to $100$ width. b. Three G-FuNK layers each with a width of 100 and using $25$ modes and each using a GELU activation function. c. Two linear layers going from a width of $200$ to $32$, with a GELU activation, and then $32$ to $1$. The network was integrated using the NeuralODE software with a forward Euler integrator with a time step of $1$ ms.

The network was trained and optimized using the PyTorch software with a relative $H_{1}$ loss function with the form $\frac{||\hat{u}-u||_{2}}{||u||_{2}}+\frac{||\nabla\hat{u}-\nabla u||_{2}}{||
\nabla u||_{2}}$. This loss was used to minimize the error at the wavefront where the difference in the gradient is typically the largest. For this case, G-FuNK was optimized with a relative-$H_{1}$ loss function $\frac{||\hat{u}-u||_{2}}{||u||_{2}}+\frac{||\nabla\hat{u}-\nabla u||_{2}}{||
\nabla u||_{2}}$. An Adam optimizer equipped with a step decay learning rate with an initial rate of $5\cdot 10^{-4}$ that halves every 100 epochs was used and was early stopped when the validation loss did not decrease after 50 epochs. The eigenvalues used in $B$ in (10) were raised to the powers $0,1,2$ in these cases. $\theta_{R}$ was a full matrix of learnable parameters.

#### A.4.3 Multiple Atrial Example

In the same set up as Cardiac Electrophysiology in 3, the model was allowed to predict from $t=[0,190]$ ms which focused on the entire action potential cycle.

For the multiple geometries case from time $[0,190]$ ms, the network is described as so, a. Linear mapping from one input ($u$) to $100$ width. b. Three G-FuNK layers each with a width of $100$ and using $25$ modes and each using a GELU activation function. c. Two Linear layers going from a width of $200$ to $32$, with a GELU activation, and then $32$ to $1$. The network was integrated using the NeuralODE software with a forward Euler integrator with a time step of $1$ ms.

The network was trained and optimized using the PyTorch software with a relative $H_{1}$ loss function with the form $\frac{||\hat{u}-u||_{2}}{||u||_{2}}+\frac{||\nabla\hat{u}-\nabla u||_{2}}{||
\nabla u||_{2}}$. This loss was used to minimize the error at the wavefront where the difference in the gradient is typically the largest. For this case, G-FuNK was optimized with a relative-$H_{1}$ loss function $\frac{||\hat{u}-u||_{2}}{||u||_{2}}+\frac{||\nabla\hat{u}-\nabla u||_{2}}{||
\nabla u||_{2}}$. An Adam optimizer equipped with a step decay learning rate with an initial rate of $5\cdot 10^{-4}$ that halves every 100 epochs was used and was early stopped when the validation loss did not decrease after 50 epochs. The eigenvalues used in $B$ in (10) were raised to the powers $0,1,2$ in these cases. $\theta_{R}$ was a full matrix of learnable parameters.

#### A.4.4 Other Example Results

| Problem | Geometries | Parameters | Training<br>Trajectories | Metrics (Test Set) | Rel $\ell_{2}$ |
| --- | --- | --- | --- | --- | --- |
| 2D Reaction Diffusion | square | 12865 | $M$=1000 | Rel. $\ell_{2}$ | 0.0796 |
| Cardiac EP, $t=[0,90]\text{ms}$ | single | 327665 | $M$=1000 | Rel. $\ell_{2}$ | 0.0717 |
| Cardiac EP, $t=[0,190]\text{ms}$ | multiple | 35640 | $M$=2400 | Rel. $\ell_{2}$ | 0.0823 |

### A.5 Computational Resources

All data generation was completed using less than $12$ CPU cores with less than $64$ GB of RAM in no more than $25$ minutes per trajectory (Computing time varies based on size of the mesh). Training of neural networks was completed using either a NVIDIA RTX A6000 with $48$ GB of virtual ram or a NVIDIA RTX A4500 with $24$ GB of virtual RAM. The Cardiac Electrophysiological example was the most extensive and took no more than five days to train.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 3
presentation: 3
questions: **Major:**

- Is there an error estimator in the prediction? How can one trust the outcome of fitted G-FUNK model for an unseen problem?
- Is there a way to enforce conservation laws (or some notion of structure preservation) in prediction if the underlying PDE admits such constraint?
- Does the solution/data need to be smooth? Can you try out viscous burger's equation with emerging discontinuity as viscosity goes to zero?
- Is the outcome ODE stable? Is there a guarantee on stability of ODE?
- Can you compare the Data-Generation/Training/Prediction time of proposed method versus the ones from a standard finite element/volume solver in the presented test cases? I believe comparing only the time/complexity of prediction against standard solver is very much misleading.
- If we know the underlying PDE, wouldn't it make sense to incorporate that in the loss, similar to what PINN does?
- Can the author show a case of extrapolation? To me, similar to PINNs, the proposed method can only be used as an efficient interpolator within the space of training data. This makes me doubt the claims on "parameter and domain-adaptation". For example, if you train your model to estimate heat equation in 1d problem inside the domain [0,1], can you test it on domain [-10,10] with different boundary conditions?

**Minor:**

- Abstract: “… for which the highest-order term in the PDE is diffusive...” What does it mean? Do you mean the highest-order term is even or, it has to be 2? I’m guessing second-order PDEs, which needs to be clarified in Abstract.
- Abstract: "without the need for retraining" and not "without need for retraining”.
soundness: 3
strengths: The paper proposes merging graphs within the standard surrogate models, which allows estimating the solution of 2nd order PDEs with varying coefficients on complicated geometries.
summary: The authors proposes an interesting surrogate model that combines Graph as the discretization method for Fourier neural operators. I believe the manuscript may be accepted after a major revision.
weaknesses: I believe there are several important details missing in the manuscript. Below, I list them.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: NA
soundness: 2
strengths: The numerical demonstration is the cardiac EP example  is interesting since I understood the least. I will rely on an expert in this area to give a meaningful comment on this example.
summary: This paper proposes a neural network model to learn a solution operator of time-dependent second-order semi-linear PDEs that takes diffusion tensor and random sample points of a family of domains as inputs.
weaknesses: This paper proposes a neural network model to learn a solution operator of time-dependent second-order semi-linear PDEs that takes diffusion tensor and random sample points of a family of domains as inputs.
While the first input has been considered by many authors, the difficulty here is to allow the domain to be chosen from a family of domains, denoted by $\\{\Omega_\alpha\\}_{\alpha\in \mathcal{A}}$. My first thought is whether the setup makes sense since there is no discussion on the class of domains that are imposed. I believe this won't work on arbitrary classes of domains. For example, e.g., for Riemannian manifolds, I would believe that if any pair of manifolds in this class have Riemannian metrics that are diffeomorphic (or even stronger, such as conformally equivalent), then the learning problem makes sense. One would need a notion of continuity between any pair of domains in the class, otherwise, it is not feasible to interpolate (to have a map that can interpolate between the training domains).  In the numerical examples shown in the paper,  there is an affine transformation between any pair of arbitrary side lengths in the 2D Nonlinear Reaction-Diffusion. For the Cardiac Electrophysiology, although the measured data come from 25 patients, the PDE is solved on processed domains (finally 24 of these), and I suspect that these domains are diffeomorphic.

The only interesting numerical demonstration is the cardiac EP example since I understood the least. I will rely on an expert in this area to give a meaningful comment on this example. In terms of methodology, I cannot understand why the Graph Laplacian structure is helpful, unless when the derivatives in Eq. (2) are defined with respect to the Riemannian metric of the embedded manifolds. It is also not obvious to me why the construction of the G-FUNK layer should be a way to go since I cannot reason it from any basic principle. 

While some numerical results look interesting, I don't really understand why the approach should work in general due to the lack of theoretical justification. I am also not sure how the approach behaves if one increases the number of layers or parameters in each G-FUNK layer. Finally, the three numerical examples are low-dimensional problems (2D or 3D), I would naively believe the standard PDE solvers should be able to solve the problem accurately in a reasonable time. Solutions to these (FEM in the cardiac EP example) are being used to train the G-FUNK model.  Based on these concerns, I believe this paper is technically (or mathematically) not interesting and not suitable for publication in ICLR. I would urge the authors to consider submitting this work to a domain science journal that is relevant to personalized cardiac EP.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. Why do the authors use the neural ODE  model? Any gain from the specific model?
soundness: 2
strengths: 1. G-FuNK leverages parameter- and domain-adapted spectral method such that it is well-suited for problems involving anisotropic diffusion and complex geometries.

2. The application on Cardiac Electrophysiology is very interesting.

3. The paper presents a detailed comparison with FNO and GNN methods, showing that G-FuNK can outperform these methods
summary: The paper proposed Graph Fourier Neural Kernels (G-FuNK), which aim to solve time-dependent, nonlinear partial differential equations (PDEs) with varying parameters and domains. G-FuNK leverages parameter- and domain-adapted spectral method. These operators are particularly well-suited for problems involving anisotropic diffusion and complex geometries. The paper demonstrates G-FuNK's effectiveness on several applications, such as heat equation simulations, reaction diffusion equations, and cardiac electrophysiology, showing promising results in terms of accuracy and computational speed.
weaknesses: 1. The limitations should be addressed.

2. The so-called graph fourier transform is actually spectral method, which need the eigenvectors pre-calculated first. This procedure could make the method not useful for large scale problem. Computational efficiency and scalability should be reported, including offline computational for the eigenvectors and online computation time. Comparison with FNO and GNN in terms of efficiency is also absent.

3. The novelty of graph fourier transform is limited. This is a widely studyed area in graph neural network community.

---

confidence: 3
contribution: 3
initial_rating: 3
presentation: 3
questions: Given that the paper focuses on PDEs on graphs, it seems that the "Multipole Graph Neural Operator for Parametric Partial Differential Equations" (MGKN), which claims to be mesh-invariant, would be a more suitable baseline for comparison. While the original MGKN framework does not explicitly tackle changing geometry, it seems it can still be applied in this setting, since only the graph is used as an input. There are also many similarities between these two approaches: the Fourier transform and inverse Fourier transform in G-FuNK play a role similar to kernel convolutions in MGKN, which are computed using the multipole algorithm. Additionally, both methods employ some form of truncation to make computation more tractable (e.g., limiting modes or long-range interactions). I would appreciate further discussion or comparison of these two approaches.
soundness: 3
strengths: The writing is clear, and the combination of graph neural networks (GNNs) and FNO is well motivated and novel to the best of my knowledge. The numerical examples show a gradual increase in complexity, leading to the final example on cardiac electrophysiology, which is complex and demonstrates the practical utility of the proposed framework.
summary: This paper proposes  Graph Fourier Neural Kernels (G-FuNK) for learning solution generators of time-dependent partial differential equations (PDEs) on graphs. G-FuNK aims to be "geometry-invariant" by leveraging the spectral domain of graphs through the Graph Fourier Transform (GFT), similar to how the Fourier Neural Operator (FNO) achieve "discretization-invariance" in regular domain.
weaknesses: In Table 1, for the three methods (G-FuNK, FNO, GNN), the number of parameters can differ by orders of magnitude, making it challenging to evaluate the improvement in performance. 

The statement "Our method predicts entire trajectories in under 1 second, significantly outperforming traditional numerical methods. For example, cardiac EP simulations typically take at least 15 minutes on 12 CPU cores for one given set of initial conditions." could be more precise, particularly regarding the numerical solution being compared. It seems the cardiac EP simulations from numerical methods serve as the high-fidelity solutions that are used to generate the training data and evaluate the error. However, G-FuNK learns on reduced modes (k_max), which may lead to limited accuracy. It would be more informative if the numerical solutions are computed on coarser mesh that achieve similar of accuracy as G-FuNK, or using the reduced modes.

Meta-review: The paper introduces a novel family of neural operators, termed Graph Fourier Neural Kernels (G-FuNK), designed to learn the temporal dynamics of diffusive PDEs across multiple anisotropic domains with varying parameters. The method embeds geometric and directional information about the domains by combining graph Laplacian-based constructions for domain-specific components with non-adapted components learned from training data using FNO. Additionally, an integrated ODE solver is used to capture the system’s time evolution. 

The paper presents an interesting approach and some compelling results, although a number of weaknesses were raised .First, there is a lack of theoretical justification for why the proposed method should work in general, particularly regarding the interplay between graph-based Laplacian constructions and the learned Fourier components. A clearer convergence analysis or performance comparison with standard methods (eg. FEM, etc..) would be necessary to substantiate some of the claims made in the paper. Second, reviewers noted that the paper lacks a quantitative understanding of the limitations of the Graph Laplacian FFT, which is central to the method. This raises concerns about scalability, stability, and general applicability beyond the presented examples. Finally, there is limited analysis regarding the impact of architectural choices (e.g., increasing the number of layers or parameters) on the method’s performance.

The panel ultimately recommends rejection. While the method demonstrates potential for domain-specific applications, the lack of theoretical justifications, convergence analysis, and scalability experiments limits its impact for the broader ICLR audience. That said, the contributions could be valuable in application-focused venues where solving PDEs on anisotropic domains is of high practical importance. The authors are encouraged to further develop the theoretical underpinnings of the approach, provide comprehensive comparisons to standard methods, and explore more comprehensively the scalability of the proposed method. With these improvements, this work could make a meaningful contribution to both the scientific computing and machine learning communities.

Decision: Reject