[table]capposition=top

# Parametric Encoding with Attention and Convolution Mitigate Spectral Bias of
Neural Partial Differential Equation Solvers

Mehdi Shishehbor These authors contributed equally to this work. Department of
Mechanical and Aerospace Engineering, University of California, Irvine  Shirin
Hosseinmardi11footnotemark: 1 Department of Mechanical and Aerospace
Engineering, University of California, Irvine  Ramin Bostanabad  Corresponding
Author<EMAIL_ADDRESS>
GitHub Repository: https://github.com/Bostanabad-Research-Group/pgcan
Department of Mechanical and Aerospace Engineering, University of California,
Irvine

Abstract
Deep neural networks (DNNs) are increasingly used to solve partial
differential equations (PDEs) that naturally arise while modeling a wide range
of systems and physical phenomena. However, the accuracy of such DNNs
decreases as the PDE complexity increases and they also suffer from spectral
bias as they tend to learn the low-frequency solution characteristics. To
address these issues, we introduce Parametric Grid Convolutional Attention
Networks (PGCANs) that can solve PDE systems without leveraging any labeled
data in the domain. The main idea of PGCAN is to parameterize the input space
with a grid-based encoder whose parameters are connected to the output via a
DNN decoder that leverages attention to prioritize feature training. Our
encoder provides a localized learning ability and uses convolution layers to
avoid overfitting and improve information propagation rate from the boundaries
to the interior of the domain. We test the performance of PGCAN on a wide
range of PDE systems and show that it effectively addresses spectral bias and
provides more accurate solutions compared to competing methods.

Keywords: Partial differential equations, Parametric encoding, Convolution,
Spectral bias, Attention mechanism.

##  1 Introduction

Partial differential equations (PDEs) are ubiquitously used to describe
various systems and physical phenomena such as the climate (Washington and
Parkinson, 2005) and turbulence (Reynolds, 1976). Traditional numerical
approaches such as the finite element method (FEM) have long been used to
solve PDEs (Meng et al., 2022) but these methods have some limitations (e.g.,
high implementation costs or dependency of the obtained solution to the
discretization) (Cuomo et al., 2022; Wang et al., 2022) that hamper their
usage for compute-intensive applications such as uncertainty quantification.
To address this challenge, surrogates such as Gaussian processes (GPs)
(Yousefpour et al., 2023) and deep neural networks (DNNs) are increasingly
used (Blechschmidt and Ernst, 2021; Mora et al., 2024).

A major issue in building accurate surrogates is that their accuracy heavily
relies on the size of their training dataset. Physics-informed neural networks
(PINNs) have shown great promise in tackling this issue by infusing domain
knowledge or constraints into the learning process. As explained in Section 2,
the loss function of a PINN differs from conventional DNNs in that it includes
additional terms to ensure that the model’s prediction satisfies the PDE as
well as the initial and boundary conditions (Raissi et al., 2019) (IC and
BCs). Despite their success, vanilla PINNs (vPINN) with multi-layer perceptron
(MLP) or fully-connected feed-forward (FCFF) architectures suffer from some
major challenges such as spectral bias (Cuomo et al., 2022; Kang et al., 2023;
Wang et al., 2021; Rahaman et al., 2019) which refers to the tendency of these
models to preferentially learn lower-frequency components of a function.

Spectral bias in FCFF networks stems from the fact that such networks have
shared parameters whose values affect the model’s output everywhere in the
domain. This issue is exacerbated when the input space is low-dimensional
which is also the case in vPINN since their inputs typically represent time
and/or space (Kang et al., 2023). Among the various types of PDEs tackled by
vPINN the so-called stiff PDEs pose a unique challenge since their solutions
can change rapidly over a small spatial domain or temporal window due to,
e.g., interplay of multiple scales, high sensitivity to IC/BCs, or conflicting
BCs that give rise to discontinuities (Sharma et al., 2023).

The training process of vPINN, which must satisfy the PDE residuals and IC/BCs
simultaneously, is similar to multi-task learning where the network’s loss
function is comprised of a (weighted) sum of multiple terms where each term
measures the ability of the network in learning a particular task. Such multi-
term loss functions also challenge training accurate vPINN since the network
minimize one term better than the others; resulting in an overall sub-optimal
solution (Daw et al., 2022). This challenge is further exacerbated by spectral
bias since the scale and gradient of the IC/BC loss terms is typically quite
different than those of the PDE residuals. While approaches for automatically
satisfying IC/BCs have been developed (Lagari et al., 2020; Lagaris et al.,
1998; Dong and Ni, 2021; McFall and Mahan, 2009; Berg and Nyström,
2018)111These approaches dispense with the IC/BC terms in the loss function.,
they typically complicate the training dynamics and can produce models whose
prediction errors considerably increase as the distance between the query
point and the IC/BCs increases (Karnakov et al., 2024).

As detailed in Section 2, recent works have found some success in addressing
the above challenges by developing novel neural architectures and training
mechanisms. In this paper, we build on these works and introduce Parametric
Grid Convolutional Attention Networks (PGCANs) that increase the ability of
vPINN in learning high-frequency solution features without resulting in over-
fitted models. As demonstrated with a host of experiments in Section 5, PGCANs
particularly improve the performance of existing technologies in approximating
PDE solutions with high gradients. Our models achieve these attractive
features using three primary ingredients. Motivated by encoding practices
adopted in neural rendering (Tewari et al., 2022), the first ingredient of
PGCAN is an encoding scheme that projects the spatiotemporal inputs into a
higher-dimensional feature space. For encoding, we use a parametric multi-cell
grid (with trainable features located at its vertices) as it provides a
structured representation of the spatiotemporal scales. Such a representation
enables capturing localized features in the PDE solution and reduces spectral
bias since not all the model parameters contribute to the predictions made for
an arbitrary query point.

The second ingredient of PGCAN are convolution layers which locally combine
the feature values of adjacent vertices of the grid. In sharp contrast to
neural rendering works where the models are trained on labeled data (Müller et
al., 2022; Dong et al., 2022), we argue that in the context of solving PDE
systems encoders with convolved features are superior to their counterparts
since convolution layers $(1)$ help to propagate the information from IC/BCs
into the domain, $(2)$ reduce the number of required grid cells (and hence
trainable parameters), and $(3)$ avoid overfitting which typically manifests
as a salt-and-pepper noise in predictions of DNNs that leverage parametric
grids as encoders.

The third ingredient of our approach is a self-attention mechanism inspired by
(Wang et al., 2021; McClenny and Braga-Neto, 2020) which enables PGCAN to more
effectively learn complex solution characteristics (e.g., high-frequency or
large gradients) by assigning more attention (i.e., weight) to the features of
the grid cell where such characteristics appear.

We summarize our contributions as follows: (1) We propose a novel architecture
that uniquely merges parametric grid encoding with self-attention mechanisms,
a fusion designed to significantly enhance the precision of PDE solutions and
mitigate spectral bias. (2) Develop a versatile model capable of effectively
solving a wide spectrum of PDEs, ranging from simple to highly challenging
scenarios. (3) Proposing a framework for assessing spectral bias in models by
comparing the Power Spectral Density (PSD) of their errors.

The rest of our paper is organized as follows. We provide some technical
background and review the relevant works in Section 2. We introduce our
approach in Section 3 and then in Section 4 develop the directional PSD-based
metric for assessing the spectral bias of DNNs when they are used to solve
PDEs. We compare the performance of PGCAN against competing methods on a host
of PDE systems in Section 5 where we also study the spectral bias of these
methods. We provide concluding remarks and future research directions in
Section 6.

###  1.1 Nomenclature

Unless otherwise stated, throughout the paper we denote scalars, vectors, and
matrices or tensors of dimensions greater than two with regular, bold lower-
case, and bold upper-case letters, respectively (e.g., $x,\boldsymbol{x},$ and
$\boldsymbol{X}$). Vectors are by default column vectors. Scalars sampled from
matrices or higher-dimension tensors are indicated by the parentheses
containing the corresponding indices. (e.g., $X(i,j)$ is a scalar sampled from
$\boldsymbol{X}\in\mathbb{R}^{n_{i}\times n_{j}}$).

##  2 Background and Related Works

Figure 1: Vanilla physics-informed neural networks (vPINNs) for solving $1D$
Burgers’ equation: The model parameters, collectively denoted by
$\boldsymbol{\theta}$, are optimized by minimizing the three-component loss
function that encourages the network to satisfy the PDE inside the domain
while reproducing the IC/BCs. These loss components are obtained by querying
the network on a set of test points that are distributed inside the domain and
on its boundaries.

As schematically demonstrated in Figure 1 for the $1D$ viscous Burgers’
equation, the essential idea of vPINN is to parameterize the solution of a PDE
system with a DNN whose weights and biases (collectively denoted by
$\boldsymbol{\theta}$) are optimized such that the DNN satisfies the
differential equations as well as the IC/BCs. To demonstrate, we consider the
following generic PDE system:

$\displaystyle\mathcal{N}_{\scriptscriptstyle{\boldsymbol{x},t}}[u(\boldsymbol{x},t)]=f(\boldsymbol{x},t),\quad\boldsymbol{x}\in\Omega,\quad
t\in[0,T],$ (1a) $\displaystyle
u(\boldsymbol{x},t)=g(\boldsymbol{x},t),\quad\boldsymbol{x}\in\partial\Omega,\quad
t\in[0,T],$ (1b) $\displaystyle
u(\boldsymbol{x},0)=h(\boldsymbol{x}),\quad\boldsymbol{x}\in\Omega,$ (1c)

where $\boldsymbol{x}=[x,y,z]$ and $t$ denote, respectively, spatial
coordinates and time which we combine as
$\boldsymbol{\zeta}=[\boldsymbol{x},t]^{\top}$ for brevity. $\Omega$ and
$\partial\Omega$ denote the spatial domain and its boundary, and
$u(\boldsymbol{x},t)$ is the PDE solution and $\mathcal{N}_{\boldsymbol{x},t}$
denotes the differential operator acting on it, $f(\boldsymbol{x},t)$ is a
known function, and the prescribed IC and BC are characterized via
$h(\boldsymbol{x})$ and $g(\boldsymbol{x},t)$, respectively. The loss function
of vPINNs is expressed as:

$\displaystyle\mathcal{L(\boldsymbol{\theta})}=\mathcal{L}_{\scriptscriptstyle
PDE}(\boldsymbol{\theta})+\lambda_{\scriptscriptstyle
BC}\mathcal{L}_{\scriptscriptstyle
BC}(\boldsymbol{\theta})+\lambda_{\scriptscriptstyle
IC}\mathcal{L}_{\scriptscriptstyle IC}(\boldsymbol{\theta}),$ (2a)
$\displaystyle\mathcal{L}_{\scalebox{0.5}{$PDE$}}(\boldsymbol{\theta})=\frac{1}{n_{\scalebox{0.5}{$PDE$}}}\sum_{i=1}^{n_{\scalebox{0.5}{$PDE$}}}\left|\mathcal{N}_{\scriptscriptstyle{\boldsymbol{x},t}}[\hat{u}(\boldsymbol{\zeta}_{i};\boldsymbol{\theta})]-f(\boldsymbol{\zeta}_{i})\right|^{2},$
(2b) $\displaystyle\mathcal{L}_{\scriptscriptstyle
BC}(\boldsymbol{\theta})=\frac{1}{n_{\scriptscriptstyle
BC}}\sum_{i=1}^{n_{\scriptscriptstyle
BC}}\left|\hat{u}(\boldsymbol{\zeta}_{i};\boldsymbol{\theta})-g(\boldsymbol{\zeta}_{i})\right|^{2},$
(2c) $\displaystyle\mathcal{L}_{\scriptscriptstyle
IC}(\boldsymbol{\theta})=\frac{1}{n_{\scriptscriptstyle
IC}}\sum_{i=1}^{n_{\scriptscriptstyle
IC}}\left|\hat{u}(\boldsymbol{x}_{i},0;\boldsymbol{\theta})-h(\boldsymbol{x}_{i})\right|^{2},$
(2d)

where $\hat{u}$ represents the predicted solution from a vPINN and the terms
$\mathcal{L}_{PDE}$, $\mathcal{L}_{BC}$, and $\mathcal{L}_{IC}$ represent the
PDE residuals and the network’s errors in reproducing the BCs and IC,
respectively. $\lambda_{BC}$ and $\lambda_{IC}$ are the weights that scale the
BC/IC loss terms to ensure that $\mathcal{L(\boldsymbol{\theta})}$ is not
dominated by $\mathcal{L}_{PDE}$.

vPINNs typically leverage FCFF networks and a loss function similar to that in
Equation 2. These choices result in spectral bias (Rahaman et al., 2019) which
actually offers an intuitive explanation for the good generalization power of
over-parameterized neural networks (Xu et al., 2019) that inherently adhere to
a class of models with lower complexity during early training stages (i.e.,
early stopping prevents overfitting). For instance, Rahaman et al. (2019)
empirically demonstrated that FCFF networks with ReLu activation functions
first learn lower frequency features that exhibit high robustness to random
perturbations in the data.

While spectral bias can be leveraged in a constructive manner to prevent
overfitting, it challenges the training process of vPINNs since high-frequency
and localized features (e.g., solution discontinuities or high gradients)
often indicate important characteristics of PDE systems. To address this
issue, a range of innovative strategies have been developed that, as we have
broadly classified below, either alter the training process or design new
architectures.

###  2.1 Training Process Adaptation

The majority of strategies in this category rely on modifying the loss
function components (or their contributions to the overall loss). Wang et al.
(2022) recently used the neural tangent kernel (NTK) theory (Jacot et al.,
2018) to demonstrate that vPINNs trained via stochastic gradient descent (SGD)
suffer from spectral bias where the network tends to approximate low-frequency
solution features that prioritize minimizing $\mathcal{L}_{pde}$ over the
other terms in Equation 2a. Following this observation, algorithms based on
minimax weighting (Liu and Wang, 2021), nonadaptive weighting (Wight and Zhao,
2020), learning rate annealing (Wang et al., 2021), self-adaptive mask
functions (McClenny and Braga-Neto, 2020), data augmentation of high-frequency
components (Fridovich-Keil et al., 2022), and NTK theory (Wang et al., 2022)
are proposed to help vPINNs in effectively minimizing all the loss components
on the right-hand-side of Equation 2a.

In our approach, we leverage the method developed by Wang et al. (2021) which
scales the BC and IC loss terms based on the PDE residuals such that the
gradients of these loss terms with respect to the network’s parameter have
comparable magnitudes. Specifically, first the dynamic weight for the $i^{th}$
loss term is calculated as:

$\hat{\lambda}_{i}=\frac{max_{\boldsymbol{\theta}_{n}}|\nabla_{\boldsymbol{\theta}}\mathcal{L}_{PDE}(\boldsymbol{\theta}_{n})|}{\overline{|\nabla_{\boldsymbol{\theta}}\mathcal{L}_{i}(\boldsymbol{\theta}_{n})|}},$
(3)

where $\boldsymbol{\theta}_{n}$ denotes the parameter values at the $n^{th}$
epoch and the bar indicates averaging. Then, a moving average (over the
training epochs) version of Equation 3 is used as the coefficient of IC and BC
losses in Equation 2a.

###  2.2 Architectural Designs

In addition to improving the training process as described in Section 2.1,
major architectural changes are needed for accurately solving complex
nonlinear PDE systems via DNNs. In general, these changes aim to improve the
backpropagated loss gradients such that the model’s learning capacity is more
effectively leveraged in approximating localized features.

One approach is to use layers that perform wavelet transforms such as Wavelet
Weighted Product (Wang et al., 2021) and Mscale networks (Liu et al., 2020)
that are frequency-aware. Another approach is based on Fourier features
(Tancik et al., 2020; Wang et al., 2021; Liu et al., 2020; Mildenhall et al.,
2021) where high-frequency modes are explicitly introduced in the network such
that it simultaneously approximates low- and high-frequency features. These
features encode the network inputs into a high-dimensional space by passing
them through sinusoidal functions of various frequencies which can be either
trained or selected from a defined spectral range. While Fourier features are
useful, choosing the appropriate frequencies for encoding is not
straightforward since the approach loses its power if the selected frequencies
do not align well with the high-frequency components of the target function.

In the following two subsections, we review in more details two recent
approaches that are also used in our model. These approaches modify the
architecture of the networks and are also used as baselines in Section 5. We
note two points in passing here. Firstly, besides input encoding, activation
functions (Hong et al., 2022; Liu et al., 2020) can also be modified to tackle
spectral bias. For instance, it has been shown that piecewise linear B-spline
or the Hat function can remove spectral bias (Hong et al., 2022) but their
effectiveness strongly depends on the PDE characteristics and network
architecture. Secondly, the combination of these approaches can also be used
to address spectral bias. For instance, Aldirany et al. (2024) iteratively
builds a series of sequentially connected FCFF networks with Fourier features
to capture the high-frequency solution features in boundary value problems.

####  2.2.1 Transformer-based Architectures: M4

Motivated by the attention mechanism, the M4 architecture (Wang et al., 2021)
integrates two transformer networks and leverages the learning rate annealing
procedure described in Section 2.1. In this model the inputs are first
projected into two high-dimensional feature spaces, $\boldsymbol{\phi}_{1}$
and $\boldsymbol{\phi}_{2}$, which are then passed through a series of hidden
layers as follows:

$\displaystyle\boldsymbol{\phi}_{1}=\sigma(\boldsymbol{W}^{1}\boldsymbol{\zeta}+\mathbf{b}^{1}),$
(4a)
$\displaystyle\boldsymbol{\phi}_{2}=\sigma(\boldsymbol{W}^{2}\boldsymbol{\zeta}+\mathbf{b}^{2}),$
(4b)
$\displaystyle\mathbf{h}^{1}=\sigma(\boldsymbol{W}^{z,1}\boldsymbol{\zeta}+\mathbf{b}^{z,1}),$
(4c)
$\displaystyle\mathbf{z}^{k}=\sigma(\boldsymbol{W}^{z,k}\mathbf{h}^{k}+\mathbf{b}^{z,k}),\quad
k=1,\dots,L,$ (4d)
$\displaystyle\mathbf{h}^{k+1}=(1-\mathbf{z}^{k})\odot\boldsymbol{\phi}_{1}+\mathbf{z}^{k}\odot\boldsymbol{\phi}_{2},\quad
k=1,\dots,L,$ (4e)
$\displaystyle\hat{u}(\boldsymbol{\zeta})=\boldsymbol{W}\mathbf{h}^{L+1}+\mathbf{b}.$
(4f)

where $\odot$ indicates element-wise multiplication, $\sigma$ is a nonlinear
activation function, and
$\boldsymbol{\theta}=\\{\boldsymbol{W}^{1},\mathbf{b}^{1},\boldsymbol{W}^{2},\mathbf{b}^{2},(\boldsymbol{W}^{z,k},\mathbf{b}^{z,k})_{k=1}^{L},\boldsymbol{W},\mathbf{b}\\}$
are the model parameters. In M4, the projections in Equation 4a and Equation
4b evolve into more complex representations before connecting with the network
outputs. These representations especially benefit from the attention mechanism
introduced via Equation 4e which blends them via the initially learnt features
$\boldsymbol{\phi}_{1}$ and $\boldsymbol{\phi}_{2}$.

####  2.2.2 Parametric Grid-based Encoding: PIXEL

Motivated by recent works in neural rendering (Mildenhall et al., 2021; Dong
et al., 2022; Lindell et al., ; Müller et al., 2022; Tancik et al., ; Xiangli
et al., ; Tewari et al., 2022; JAMBON et al., 2023), PIXEL (Kang et al., 2023)
leverages parametric grid-based encoding to map the spatiotemporal inputs into
a high-dimensional learnable feature space which is then connected to the
output (i.e., the PDE solution) via a relatively shallow FCFF network. The
rationale behind adopting such an encoding in neural rendering and PIXEL is
twofold. Firstly, the domain is structured by grids of $N_{r}$ resolutions
where the grid corresponding to resolution $l\in\\{1,..,N_{r}\\}$ has
$N_{v}^{l,i}$ vertices along dimension $i\in\\{x,y,z,t\\}$. The neighboring
vertices identify a cell wherein the vertex features focus on capturing the
local patterns. Secondly, it decreases computational costs (while increasing
the memory footprint) since a backpropagated gradient associated with
$\boldsymbol{\zeta}$ only updates the vertex features of the cell containing
$\boldsymbol{\zeta}$ (as opposed to an FCFF network whose parameters are all
updated by any backpropagated gradient).

The encoder of PIXEL largely relies on another work (Müller et al., 2022)
where multi-resolution hash encoding is used for fast and accurate image
reconstruction and neural rendering. Hence, we first review the ideas
developed by Müller et al. (2022) and then explain the working principles of
PIXEL. Multi-resolution hash encoding maps the inputs to a multiscale grid-
based feature space (using hash tables and interpolations) and then passes
these features through a shallow FCFF decoder to obtain the outputs. This
process is illustrated in Figure 2 (the use of hash table is not shown) where
the cell boundaries are marked with solid black and dashed red lines at the
first (coarse) and second (fine) resolutions, respectively. The vertices in
the first and second resolutions in Figure 2 are color-coded by black and red
markers, respectively. For the arbitrary query point $\boldsymbol{\zeta}$,
Müller et al. (2022) first locate the cell at each resolution that contains it
and then linearly interpolate the vertex features to obtain the feature vector
of $\boldsymbol{\zeta}$. Afterwards, the feature vectors obtained from the
different resolutions are concatenated and fed into the shallow FCFF network
for decoding.

Figure 2: Parametric grid encoding with $N_{r}=2$ in $2D$: The grid encoding
has $N_{r}=2$ levels which have $3\times 3$ and $9\times 9$ cells at the
coarse and fine resolutions. Each vertex at any resolution is endowed with
some learnable features which are used to obtain the features of the query
point $\boldsymbol{\zeta}$ via interpolation. $\bar{x}$ and $\bar{y}$ denote
the local coordinates of $\boldsymbol{\zeta}$ in the cell that contains it.

The architecture of PIXEL (Kang et al., 2023) resembles the one in Figure 2
with two major modifications. Firstly, PIXEL leverages a high-order
interpolation function rather than a linear one to ensure gradients with
respect to the inputs are properly backpropagated and achieve $C^{1}$
continuity at the vertices (see Section 3.2 for more details). Secondly, PIXEL
uses a single relatively fine resolution to solve PDE systems. Such a
resolution provides the ability to learn high-frequency and localized features
where the PDE solution characteristics can substantially vary across the grid
cells. However, it might result in overfitting and hence Kang et al. (2023)
address it via multi-grid techniques where $n_{rep}$ grids of the same
resolution are used to find $n_{rep}$ feature vectors for any query point.
These repeated grids are obtained by diagonally shifting the original one and
the final feature vector of a query point is obtained by summing the $n_{rep}$
feature vectors corresponding to the $n_{rep}$ grids. While these
modifications enable PIXEL to solve a wide range of PDE systems, they render
it computationally very expensive and still prone to overfitting. This is
especially the case when solving complex PDE systems such as the Navier-Stokes
equations.

##  3 Our Approach: PGCAN

We argue that the parametric grid-based encoding of Müller et al. (2022) acts
as a double-edged sword when it is used as part of a model that aims to solve
PDE systems. On the one hand, it provides a localized learning ability because
during training the features associated with a cell are only updated by
backpropagated gradients of query points that lie in that cell. This ability
indicates that adjacent cells can have quite different features which, in
turn, enables the model to capture high-frequency and sharp changes across
different cells. Using a multi-resolution encoder (see Figure 2 for an
example) particularly helps in this regard as coarse and fine resolutions can
focus on learning low- and high-frequency behaviors.

On the other hand, a grid-based encoder heavily relies on the decoder (which
is typically a shallow FCFF network in neural rendering applications, see for
example Tewari et al. (2022)) to relate the feature values to the output. Such
an encoder-decoder setup, however, dramatically loses its accuracy in solving
PDEs when there are no labeled data points in the domain. This issue is a
direct result of the localization ability of the encoder which largely
prevents information propagation from the IC/BCs to the rest of the domain.
While multi-resolution encoders have the potential to address this issue, they
fail in practice as they quickly converge to local optima during training. In
our studies, we observed a similar issue when the grids of a multi-resolution
encoder (with shared or independent decoders) are trained sequentially from
coarse to fine resolutions or vice versa. We also note that increasing the
depth of the decoder not only negates the entire point of having an encoder,
but also increases the computational costs.

Figure 3: Parametric Grid Convolutional Attention Networks (PGCANs) for
solving the Navier-Stokes equations: encoder-decoder setup (a) The input space
is mapped to a structured high-dimensional space parameterized via features
$\boldsymbol{F}_{0}^{l}\in\mathbb{R}^{n_{rep}\times N_{f}\times
N_{v}^{l,x}\times N_{v}^{l,y}}$. Upon convolution and interpolation for a
query point, these features are passed to the ensuing NN that implements the
projections formulated in Equation 4e. Feature convolution block (b):
Trainable features are arranged in grids covering the domain. These features
are convolved by a 3$\times$3 kernel and followed by a $\tanh$ activation
function. Feature interpolation block (c): The query point is placed in the
convolved feature maps. Each of these grid-like maps is diagonally shifted to
prevent overfitting. Cosine interpolation is then performed based on the local
coordinates of the point in the corresponding unit cells.

The multi-grid technique of PIXEL (Kang et al., 2023) alleviates the above
issues but is very expensive and still prone to overfitting since it relies on
$(1)$ high-resolution grids which have large feature vectors associated with
their vertices, and $(2)$ a large number of grid repetitions (e.g.,
$n_{rep}=96$ or $n_{rep}=192$). As demonstrated in Section 5, this overfitting
issue is pronounced when solving complex PDE systems such as the Navier-Stokes
equations.

As schematically illustrated in Figure 3 and detailed in Sections 3.1, 3.2 and
3.3, we address the above issues with two primary innovations. Firstly, we
introduce convolutional layers in the encoder to propagate boundary
information (from IC/BCs) into the domain more effectively. The presence of
these layers provides PGCAN with a natural multi-scale representation of the
domain similar to how convolutional layers provide a cascade of
representations in vision models (Szegedy et al., ; Chollet, ; Howard et al.,
2017). It also prevents overfitting to the extent that our encoder has a
single mother resolution222As opposed to the encoder of Müller et al. (2022)
that uses multiple distinct resolutions which has much fewer cells and grid
repetitions compared to PIXEL.

Our second innovation is to use a transformer-type network as the decoder
instead of a shallow FCFF one. This choice is justified as it enables the
decoder to have access to the encoder features at multiple depths to improve
gradient flows. We note that this decoder enables our approach to naturally
extend to higher dimensions and irregular domains. We demonstrate these
attractive features in Section 5.

###  3.1 Parametric Grid with Convolution

PGCAN starts with encoding the input space with the set of features
$\boldsymbol{F}_{0}^{l}$ that partition it into a collection of cells whose
vertices are endowed with trainable features. The encoder characteristics
(e.g., resolution and size of vertex features) are chosen to strike a balance
between computational costs and the ability of the model to learn high-
frequency features (in Section 5 we use the same settings across all examples
to demonstrate the robustness of our approach to these choices). Following the
grid definition, we employ a convolutional layer to $(1)$ bidirectionally
enhance the flow of information between the boundaries and the interior
regions, and $(2)$ capture spatial correlations between features and avoid
overfitting. Our studies show that this convolution layer greatly improves the
gradient flows and avoids trivial solutions that only satisfy the PDE
residuals without respecting the IC/BC.

In both PGCAN and PIXEL $N_{r}=1$ so we omit the superscripts $l$ hereafter.
With a single resolution parametric grid in a $2D$ spatial domain333In our
approach, we treat the temporal dimension identically to spatial ones., we
parameterize the domain via $\boldsymbol{F}_{0}\in\mathbb{R}^{n_{rep}\times
N_{f}\times N_{v}^{x}\times N_{v}^{y}}$ features and index the vertices via
$\boldsymbol{i}=[i_{x},i_{y}]^{\top}$ where $i_{x}\in\\{0,1,..,N_{v}^{x}-1\\}$
and $i_{y}\in\\{0,1,..,N_{v}^{y}-1\\}$. The feature map
$\boldsymbol{F}_{0}^{p,q}\in\mathbb{R}^{N_{v}^{x}\times N_{v}^{y}}$ refers to
the spatial parametric grid isolated at repetition $p\in\\{1,..,n_{rep}\\}$
and feature index $q\in\\{1,..,N_{f}\\}$. A feature-wise $2D$ convolution
followed by a $\tanh$ activation function transforms the feature located at
vertex $\boldsymbol{i}$ as:

$\begin{split}F_{c}^{p,q,i_{x},i_{y}}&=\tanh(\sum_{b=0}^{2}\sum_{a=0}^{2}\\\
&F_{0}^{p,q,i_{x}+a-1,i_{y}+b-1}\cdot W_{q}(a,b))\end{split}$ (5)

where $\boldsymbol{W}_{q}\in\mathbb{R}^{3\times 3}$ is the kernel weight
matrix associated with the $q^{th}$ feature. The convolved feature maps
$\boldsymbol{F}_{c}^{p,q}$ are finally scattered in the feature map
$\boldsymbol{F}_{c}\in\mathbb{R}^{n_{rep}\times N_{f}\times N_{v}^{x}\times
N_{v}^{y}}$ which is used for feature interpolation in our approach (as
opposed to PIXEL which uses $\boldsymbol{F}_{0}$). We note that the above
operations can use arbitrary kernel sizes (we use $3\times 3$ in our studies)
and naturally extend to higher input dimensions (including time) and multi-
resolution grids.

###  3.2 Differentiable Feature Interpolation

To obtain the feature vector of a query point, we first identify the cell that
it lands in (e.g., the black square in Figure 2) and then obtain that point’s
local coordinates in the identified cell. These local coordinates are denoted
by $\bar{\boldsymbol{x}}=(\bar{x},\bar{y})$ and fall in the $[0,1]^{2}$ range
since the vertices of the cell are mapped to
$\boldsymbol{c}\in\\{(0,0),(0,1),(1,0),(1,1)\\}$. Following Kang et al.
(2023), these local coordinates are further transformed via Equation 6 to
ensure the model’s output is differentiable with respect to its inputs:

$\tilde{\boldsymbol{x}}=\frac{1}{2}(1-\cos(\pi\bar{\boldsymbol{x}}))\in[0,1]^{2}.$
(6)

Once $\tilde{\boldsymbol{x}}$ corresponding to the query point
$\boldsymbol{x}$ is obtained, its feature vector
$\mathbf{f}^{p}(\boldsymbol{x})\in\mathbb{R}^{N_{f}}$ is computed for each
grid repetition $p$ by linearly interpolating the trainable feature vectors
$\mathbf{f}_{\boldsymbol{c}}^{p}\in\mathbb{R}^{N_{f}}$ of the vertices of the
cell that contains $\tilde{\boldsymbol{x}}$, that is:

$\begin{split}\mathbf{f}^{p}(\tilde{\boldsymbol{x}})=&(1-\tilde{x})(1-\tilde{y})\mathbf{f}_{(0,0)}^{p}+(1-\tilde{x})\tilde{y}\mathbf{f}_{(0,1)}^{p}\\\
&+\tilde{x}(1-\tilde{y})\mathbf{f}_{(1,0)}^{p}+\tilde{x}\tilde{y}\mathbf{f}_{(1,1)}^{p}.\end{split}$
(7)

Finally, the resulting $n_{rep}$ feature vectors are summed up to obtain the
feature vectors of the query point in the global coordinate:

$\mathbf{f}(\boldsymbol{x})=\sum_{p=1}^{n_{rep}}\mathbf{f}^{p}(\tilde{\boldsymbol{x}})\in\mathbb{R}^{N_{f}}.$
(8)

The feature vector $\mathbf{f}(\boldsymbol{x})$ continuously changes in the
domain and is informed by features outside of the cell that the query point
lands in via two mechanisms: the convolution operations and the multi-grid
scheme. $\mathbf{f}(\boldsymbol{x})$ is eventually split into two equal-size
vectors $\mathbf{f}_{1}(\boldsymbol{x})$ and $\mathbf{f}_{2}(\boldsymbol{x})$
which are passed to the decoder.

We highlight that the above descriptions are based on spatial coordinates but
they can be naturally extended to the temporal dimension.

###  3.3 Decoder with Attention Mechanism

We connect the learnable features $\boldsymbol{F}_{c}^{p,q}$ to the outputs
via a decoder that leverages an attention mechanism to dynamically prioritize
the features that are more predictive of the PDE solution. The decoder in
Figure 3 resembles the one used in the M4 model discussed in Section 2.2.1
with the major difference that unlike M4 (which generates
$\boldsymbol{\phi}_{1}$ and $\boldsymbol{\phi}_{2}$ directly from the inputs
$\boldsymbol{\zeta}$ via an FCFF layer) we use the feature vectors
$\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ that are generated by our encoder. As
demonstrated in Section 5, this choice enables PGCAN to learn spatiotemporal
features that are localized and yet informed about the IC/BCs due to the
convolution operation and the multi-grid scheme embedded in the encoder. We
highlight that a byproduct of using such a decoder (as opposed to a shallow
FCFF one) is that the learnt feature maps may or may not represent the PDE
solution. This behavior is in sharp contrast to the results reported in Müller
et al. (2022) and Kang et al. (2023) where the learnt feature maps resemble
the model’s output.

##  4 Spectral Bias Quantification

In this section we develop quantitative metrics to evaluate our assertion that
PGCAN mitigates the spectral bias of neural PDE solvers. To this end, we
leverage the Fourier transform and (directional) power spectral density (PSD)
which are ubiquitously used in signal processing for identifying patterns and
structures in the frequency domain.

Figure 4: Directional PSD for analytic signals: Unlike Gaussian and uniform
noise, the directional PSD curves of signals with spatially varying
frequencies are not flat.

The Fourier Transform decomposes a signal such as an image into its sine and
cosine components. For the square image $I(x,y)$ of size $N\times N$, the
zero-centered $2D$ discrete Fourier Transform $\mathcal{F}(u,v)$ is calculated
as:

$\mathcal{F}(u,v)=\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}I(x,y)\cdot e^{-2\pi
i\left(\frac{xu}{N}+\frac{yv}{N}\right)}$ (9)

where $u$ and $v$ represent the spatial frequency indices for the $x$ and $y$
axes, respectively. If $N$ is odd,
$u,v\in\\{-\frac{N-1}{2},...,-1,0,1,...,\frac{N-1}{2}\\}$ and otherwise
$u,v\in\\{-\frac{N}{2},...,-1,0,1,...,\frac{N}{2}-1\\}$. Given
$\mathcal{F}(u,v)$, we can quantify the strength of its various frequencies
via the power spectrum $P(u,v)$ which is the squared magnitude of
$\mathcal{F}(u,v)$:

$P(u,v)=|\mathcal{F}(u,v)|^{2}.$ (10)

To obtain the PSD along the $x$ axis, we first define the frequency bins
$b_{s}=[s-0.5,s+0.5)$ where $s\in\\{1,2,...,\lfloor\frac{N}{2}\rfloor\\}$. In
this case, the $x$ frequency components are selected as:

$K_{x}=\\{u|u\in b_{s}\\},$ (11)

so we compute the PSD along the $x$ axis by averaging the power spectrum over
the $y$ axis’ frequencies for each $k_{x}\in K_{x}$. That is:

$\bar{P}(k_{x})=\frac{1}{N}\sum_{v}P(k_{x},v).$ (12)

n our analysis, our focus is on the frequency fluctuations in the directional
PSD of an image that quantifies the prediction errors of a model in solving a
PDE system. Hence, we exclude the effect of total error (i.e., the total
power) as described below to be able to compare error maps associated with
different models.

Assuming $\Gamma$ directional PSDs $\bar{P}(k_{x})^{\gamma}$ for
$\gamma\in\\{1,2,..,\Gamma\\}$ are obtained from $\Gamma$ different images, we
first find $\bar{P}_{1,max}=\max\\{\bar{P}(1)^{\gamma}\\}$ and then shift all
$\bar{P}(k_{x})^{\gamma}$’s to start from the same point. That is:

$PSD_{x}^{\gamma}(k_{x})=\bar{P}(k_{x})^{\gamma}+(\bar{P}_{1,max}-\bar{P}(1)^{\gamma})$
(13)

Following a similar procedure we obtain $PSD_{y}$. In Figure 4 the directional
PSD curves of a few signals are obtained following the above procedure. As it
can be observed, the curves corresponding to the noisy signals are flat;
indicating that there is no spectral bias in the underlying signal. However,
this is not the case if there are directions with dominant frequencies in a
signal.

##  5 Results and Discussions

Problem | Parameter | PDE | IC/BC | PGCAN | M4 | PIXEL | vPINN
---|---|---|---|---|---|---|---
Burgers |  | $\nu=\frac{1}{\pi}$
---
$\nu=\frac{0.01}{\pi}$
$u_{t}+uu_{x}-\nu u_{xx}=0$ |  | $u(0,x)=-\sin(\pi x)$
---
$u(-1,t)=u(1,t)$
$u_{x}(-1,t)=u_{x}(1,t)$
| 2.90E-03
---
1.45E-02
| 2.91E-03
---
1.38E-02
| 6.46E-02
---
8.26E-02
| 7.55E-03
---
1.49E-02
Convection |  | $\beta=5$
---
$\beta=30$
$u_{t}+\beta u_{x}=0$ |  | $u(x,0)=\sin x$
---
$u(0,t)=u(2\pi,t)$
| 9.05E-04
---
7.77E-03
| 8.07E-04
---
6.51E-01
| 1.49E-03
---
1.57E-01
| 1.27E-03
---
1.27E-01
Helmholtz |  | $a_{2}=1$
---
$a_{2}=10$
$\nabla^{2}u+a^{2}u=q(x,y)$ | $u(x,y)=0,(x,y)\in\partial[-1,1]^{2}$ |  | 5.94E-04
---
2.70E-03
| 1.23E-03
---
6.25E-02
| 3.44E-02
---
2.37E-01
| 9.07E-02
---
1.73E+00
LDC |  | $A=1$
---
$A=5$
| $\nabla.\boldsymbol{u}=0$
---
$\rho(\boldsymbol{u}.\nabla)\boldsymbol{u}=-\nabla
p+\mu\nabla^{2}\boldsymbol{u}$
| $\boldsymbol{u}(0,y)=\boldsymbol{u}(1,y)=(0,0)$
---
$\boldsymbol{u}(x,0)=(0,0)$
$\boldsymbol{u}(x,1)=(A\sin(\pi x),0)$
$p(0,0)=0$
| 1.22E-03
---
1.42E-02
| 1.19E-02
---
2.79E-01
| 2.56E-01
---
6.25E-01
| 3.84E-02
---
6.61E-01
Table 1: PDE systems and summary of the results: The median of $L_{2}^{re}$
errors across $10$ repetitions are reported after $50k$ training epochs. For
the LDC problem, $L_{2}^{re}$ is reported for the velocity magnitude. Our
approach consistently provides more accuracy as the complexity of the
underlying PDE system increases (the version in the second row of each problem
is more difficult to solve).

We evaluate the performance of PGCAN on four PDE systems against three models:
(1) vPINNs (Raissi et al., 2019), (2) M4 (Wang et al., 2021), and (3) PIXEL
(Kang et al., 2023). In our experiments, we fix the architecture and training
mechanism of each model to assess the effect of PDE systems’ complexities on
their accuracy. Specifically, PGCAN uses $9\times 9$ grid (with $N_{f}=128$
and $n_{rep}=2$) and a three-layer decoder with $64$ neurons, vPINN has $8$
hidden layers (each with $40$ neurons), M4 has $4$ hidden layers (each with
$40$ neurons), and PIXEL uses a $16\times 16$ encoding grid (with $N_{f}=4$
and $n_{rep}=96$) and a single-layer FCFF decoder with $16$ neurons. The
number of trainable parameters of these four models (for a single output PDE)
are roughly $35$, $12$, $7$, and $98$ thousands, respectively (as shown below
with some sample studies, the accuracy of vPINNs and M4 models is
insignificantly affected if their size is increased to match the size of our
model).

In all of our experiments, we use Adam as the optimizer and measure the
accuracy based on relative $L_{2}$:

$L_{2}^{re}=\frac{\|\widehat{\boldsymbol{u}}-\boldsymbol{u}\|_{2}}{\|\boldsymbol{u}\|_{2}},$
(14)

where $\widehat{\boldsymbol{u}}$ and $\boldsymbol{u}$ represent the vectors of
predicted and reference solutions, respectively (see Figure 5 for reference
solution fields). We train each model $10$ times (with random initializations)
for $50k$ epochs and evaluate the $L_{2}^{re}$ every $5k$ epoch (below we
report the median error values across the $10$ repetitions and provide more
details in Detailed Error Statistics). The learning rate is initialized at
$1\times 10^{-3}$ and reduced by $90\%$ every $10k$ epochs. We use $20k$
collocation points for all experiments except the lid-driven cavity (LDC)
problem where $5k$ points are used in all models. These points are resampled
every $100$ epochs to improve accuracy. Both M4 and PGCAN leverage the dynamic
weight balancing algorithm explained in Section 2.1 every $100$ epochs (except
for Helmholtz for which the weights are updated every single epoch).

The four PDE systems studied in this work are summarized in Table 1. For each
system, we consider two cases by changing a parameter that increases the
complexity of the underlying system: reducing $\nu$ in Burgers’ introduces
features that resemble a shock wave in the solution, increasing $\beta$ and
$a_{2}$ in the convection and Helmholtz problems increases the solution
frequency, respectively, and changing $A$ in the LDC problem significantly
affects the velocity gradients and vortices.

###  5.1 Summary of Comparative Studies

Figure 5: Reference solutions and error maps: The first column includes the
reference solutions while the rest of them indicate the absolute errors
associated with each model. Due to similarity, we only show the error maps of
one out of $10$ repetitions.

The results of our studies are summarized in Table 1 and indicate that, except
for the Burgers’ problem, our approach outperforms other methods by almost an
order of magnitude (we provide error statistics beyond the median in Detailed
Error Statistics). In particular, as the solution gradients or frequency
increase, we observe that PGCAN’s close competitor, PIXEL, suffers from
overfitting. The M4 model performs better than PIXEL but also loses accuracy
as PDE complexity increases. These trends are more clearly illustrated via the
error maps in Figure 5 which indicate that M4, PIXEL, and vPINN fail in
learning high-frequency features in the case of Helmholtz and Convection PDEs.
Similarly, they fail to capture the large gradient changes across the primary
vortex in the LDC problem.

We observe in Table 1 that the performance of all models decreases as the
complexity of the PDE system increases (compare the errors associated with a
model in the first and second rows). This trend is expected since the
architecture and training mechanism of all the models is fixed for a given PDE
system and hence the performance drops as $\frac{1}{\nu},\beta,a_{2}$, or $A$
increase.

In the case of Burgers’ equation, M4 marginally outperforms PGCAN as $\nu$
decreases. We attribute this behavior to the fact that PGCAN aims to learn
both local and global features; as $\nu$ decreases the sharp gradients around
the $x=0$ line become more and more localized which forces PGCAN to focus on
learning a smaller and smaller region (see the error maps of PGCAN in Figure 5
for the two values of $\nu$). Such a focus, slightly decreases the solution
accuracy elsewhere which reduces the overall performance.

An interesting behavior is also observed in the case of the convection problem
whereas $\beta$ increases from $5$ to $30$, unlike PGCAN whose performance
drops by only one order of magnitude, the performance of other models
especially M4 drops substantially (a similar trend is noticed in Helmholtz).
We believe these performance drops are due to the high-frequency nature of the
PDE solution (see the reference solutions in Figure 5) which renders the PDE
residuals quite sensitive to the network parameters (this sensitivity can be
quantified based on the eigendecomposition of the Hessian of the PDE residuals
with respect to $\boldsymbol{\theta}$). This sensitivity (which is decreased
in PGCAN due to its encoder) hampers the optimization process and reduces the
model performance.

###  5.2 Model Features

To gain more insights into the training mechanism and performance of these
models, in Figure 6 we visualize the median $L_{2}^{re}$ of each one during
$50k$ training epochs. We observe that in most cases PGCAN provides a higher
accuracy at any stage during the training process. This behavior is especially
noticeable in complex cases such as the LDC problem where having more than
$50k$ epochs seems to have benefited PGCAN more.

In the case of the LDC problem, Figure 6 contains two additional M4 and vPINN
models whose error histories are indicated with dashed lines. The only
difference between these two models compared to the original M4 and vPINN
networks (whose performance is shown via solid lines) is their size.
Specifically, we increase the size of these models such that the number of
their trainable parameters roughly matches PGCAN. As it can be observed,
merely increasing the size of the models fails to improve their performance.

Figure 6: Median of $L_{2}^{re}$ over training epochs: In most cases, PGCAN
achieves smaller errors for a fixed number of training epochs. We also observe
better error reduction rates for PGCAN in complex PDE systems such as the LDC
problem. The dashed lines in the case of LDC problem correspond to M4 and
vPINN models whose number of trainable parameters are approximately the same
as PGCAN. Figure 7: Resolution invariance: The accuracy of each model is
compared against reference solutions obtained via the finite element method
whose discretization is refined as $N$ increases. The errors converge at a
large enough $N$ after which the reference solution does not change.

We highlight that even though PGCAN discretizes the domain with a parametric
grid, it provides resolution-invariant solutions because its encoder is
followed by a continuous decoder. To demonstrate this feature, we obtain the
reference solution of the LDC problem via the finite element method while
refining the discretization by increasing the spatial discretization parameter
$N$. Once these series of reference solutions are obtained, we compare them to
the outputs provided by each model. As shown in Figure 7 for both $A=1$ and
$A=5$, the reported errors for each model converge once the reference solution
is stabilized.

Compared to neural rendering works and PIXEL, PGCAN has a more complex decoder
which, while improving the gradient flows and the network’s overall
performance, results in features maps that do not necessarily resemble the PDE
solution. To demonstrate this behavior, in Figure 8 we visualize the reference
solution for the Burgers’ problem along with the average interpolated feature
values for the $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ tensors shown in Figure
3. As it can be observed, the features continuously vary in the domain and
only resemble the solution in the case of $\mathbf{f}_{2}$. We note that each
spatial feature map in Figure 8 demonstrates the average of $64$ maps and
hence there is a chance one or more features maps in $\mathbf{f}_{1}$ and/or
$\mathbf{f}_{2}$ better resemble the reference solution. Such maps can be
found by solving a simple optimization problem but we opt out of this process
as in our approach we do not require the feature maps to resemble the PDE
solution.

(a) Average feature values in $\mathbf{f}_{1}$

(b) Average feature values in $\mathbf{f}_{2}$

(c) Reference solution

Figure 8: Interpolated features of PGCAN and the reference solution for the
Burgers’ problem: We visualize the average feature values in $\mathbf{f}_{1}$
and $\mathbf{f}_{2}$ in (a) and (b), respectively (we use the feature values
after $50000$ training epochs). The reference solution is shown in (c) which
slightly resembles the plot in (b) but not in (a).

Lastly, we study the scalability of PGCAN to higher input dimensions and
irregular domains as real world applications typically involve such scenarios.
To this end, we consider the torus in Figure 9(a) and solve Poisson’s equation
in it:

$\displaystyle\nabla^{2}u=-e^{x+y+z},\quad\boldsymbol{x}\in\Omega,$ (15a)
$\displaystyle u=\sin(\pi xyz),\quad\boldsymbol{x}\in\partial\Omega.$ (15b)

where $u(x,y,z)$ denotes the potential, $\nabla^{2}$ is Laplacian operator,
and $\partial\Omega$ is given by $(1-\sqrt{x^{2}+y^{2}})^{2}+z^{2}=0.5^{2}$.

As shown in Figure 9(a) our grid-based encoder does not have to conform to the
topology of the domain since the grid is essentially a mapping function that
associates inputs with a high-dimensional feature vector. Compared to an
encoder adapted to the domain shape, our strategy may result in an encoder
with more parameters which are not trained/used if their vertex belongs to a
cell that does not cover any part of the domain. However, we still prefer a
grid-based encoder whose axes are parallel to the input coordinates as it is
much easier to implement while providing a computationally efficient mechanism
for locating the cell that contains an arbitrary query point.

The performance of PGCAN in solving Equation 15 is compared to the reference
solution on selected $2D$ cross-sections in Figure 9(b) which indicates that
the errors are quite small. To put this performance into perspective, we
obtain $L_{2}^{re}$ for all the models during the training. As seen in Figure
9(c), PGCAN consistently outperforms other models during the entire training
process. We attribute this performance gap to PGCAN’s effective use of its
parametric encoder which accelerates the convergence as well (compare the
errors at $5k$ epochs).

(a) Torus and the overlaid grid.

(b) Prediction errors of PGCAN.

(c) Median $L_{2}^{re}$ across $50k$ training epochs.

Figure 9: Poisson equation solved in a torus: (a) The parametric grid-based
encoder of PGCAN contains the torus. Features whose vertices are located
outside of the domain marginally contribute to solution approximation. (b)
Point-wise error of PGCAN against reference solution, $u-\hat{u}$, on six
cross sections $y\in\\{-1.2,-0.6,0.,0.6,1.2\\}$. (c) Based on $L_{2}^{re}$,
PGCAN outperforms other models at any epoch during training.

###  5.3 Spectral Bias

We have claimed PGCANs can mitigate the spectral bias of vPINN and hence in
this section we test this claim. To this end, we use PSD as detailed in
Section 4 to see if the error signals associated with a model’s predictions
are dominated by any specific frequency ranges.

Specifically, we obtain the PSD of error maps generated by comparing the
predicted and reference solutions at uniformly spaced points in the domain. As
detailed in Section 4, we obtain the PSD curves along two orthogonal input
dimensions and shift their starting points so that the curves corresponding to
all models in each problem have the same initial values (these shifts are done
since the overall accuracy of these models and hence the power of the
resulting error maps are quite different. Since in this section our focus is
on the error frequency characteristics and not error magnitudes, these shifts
improve the interpretability of the results). With this setup, the flatness of
a PSD curve provides a measure for spectral bias where flatter curves indicate
that the underlying error maps are not dominated via high-frequency terms
(i.e., there is no spectral bias) and the underlying error map more closely
resembles a Gaussian noise (see Figure 4).

The results of our studies are summarized in Figure 10 and indicate that the
predictions obtained by PGCAN result in flatter PSD curves overall (across
both directions); indicating that our model is more effective in mitigating
spectral bias (at least in these two directions. Note that for some problems
such as the convection one the dominant frequency direction is not along any
of the axes). While M4 also performs quite well based on this metric, it fails
to reduce the overall power (i.e., error). We believe that these findings are
in line with previous research that establish the spectral bias of vPINN but
our metrics (while needing the reference solutions) are more general and do
not assume that the underlying solution is, e.g., a summation of sinusoidal
terms.

Figure 10: PSD curves of the error maps along two orthogonal input dimensions
(denoted by $x$ and $y$): Flat curves are desired as they indicate that the
underlying error map is not dominated by high- or low-frequencies. In general,
PGCAN provides flatter PSD curves (the insets magnify parts of the plot). For
Burgers’ equation, $y$ denotes time and for LDC we focus on the horizontal
velocity component.

##  6 Conclusions and Future Works

In this paper, we introduced PGCAN for solving PDE systems with DNNs that
leverage parametric grid encoding as well as convolution and attention
mechanisms. We tested the performance of PGCAN against three other methods on
a wide range of problems and demonstrated that our method provides higher
accuracy as the complexity of the underlying problem increases. Using
Poisson’s equation, we also showed that PGCAN naturally extends to higher
dimensions and irregular domains.

While developing PGCAN, we encountered several challenges such as effective
handling of IC/BCs. We hypothesized that the performance of the model would
improve by constraining certain features to closely resemble the IC/BC (or by
adding additional loss components as opposed to constraining). However, this
approach did not meet our expectations as it increased the complexity of
minimizing the loss function. Another significant challenge was overfitting
(or converging to a trivial solution that does not satisfy the IC/BCs) which
we eventually addressed via the convolution layers. Our initial studies relied
on attention mechanism (to select the most important features), feature
engineering, and variants of multi-grid techniques that increase the non-
locality feature of the encoder. Except for the convolution layers, none of
these techniques performed consistently well across different PDE systems and
hence we do not adopt them in our baseline model.

Currently, we believe the major limitation of our approach is the uniform
partitioning of the input space into cell. As demonstrated in Section 5 for
the Burgers’ equation, this feature may prevent the model from reducing the
global accuracy. We believe this shortcoming can be addressed via adaptive
domain decomposition (similar to adaptive meshing) or perhaps using spatially
varying weights for collocation points as suggested in McClenny and Braga-Neto
(2020).

## Appendix

### Detailed Error Statistics

We provide detailed error statistics associated with our comparative studies
in Table 2. These statistics are obtained for the error metric in Equation 14
across $10$ repetitions. We observe that in the majority of the cases, the
mean and median values are very close and the standard deviations are
relatively small.

Table 2: Statistics of $L_{2}^{re}$ associated with PINN, PIXEL, M4, and PGCAN: Each model is trained $10$ times on each problem to obtain the statistics of $L_{2}^{re}$. Model | mean | std | median | best | mean | std | median | best
---|---|---|---|---|---|---|---|---
| Burger
| $\nu=0.01/\pi$ | $\nu=1/\pi$
PGCAN | 1.48E-02 | 8.63E-04 | 1.45E-02 | 1.39E-02 | 2.90E-03 | 6.33E-05 | 2.90E-03 | 2.81E-03
M4 | 1.44E-02 | 1.18E-03 | 1.38E-02 | 1.37E-02 | 3.21E-03 | 7.91E-04 | 2.91E-03 | 2.72E-03
PIXEL | 9.67E-02 | 5.73E-02 | 8.26E-02 | 3.44E-02 | 6.57E-02 | 1.19E-02 | 6.46E-02 | 4.75E-02
vPINN | 1.58E-02 | 2.99E-03 | 1.49E-02 | 1.36E-02 | 9.73E-03 | 5.02E-03 | 7.55E-03 | 4.04E-03
| Convection
| $\beta=5$ | $\beta=30$
PGCAN | 9.01E-04 | 8.87E-05 | 9.05E-04 | 7.12E-04 | 7.56E-03 | 5.31E-04 | 7.77E-03 | 6.47E-03
M4 | 1.50E-03 | 1.50E-03 | 8.07E-04 | 5.51E-04 | 5.58E-01 | 2.34E-01 | 6.51E-01 | 9.31E-02
PIXEL | 1.36E-03 | 6.73E-04 | 1.27E-03 | 7.59E-04 | 1.29E-01 | 2.81E-02 | 1.27E-01 | 9.24E-02
vPINN | 4.79E-03 | 1.69E-03 | 4.39E-03 | 2.50E-03 | 9.46E-02 | 6.25E-02 | 7.75E-02 | 1.92E-02
| Helmholtz
| $a_{2}=1$ | $a_{2}=10$
PGCAN | 6.30E-04 | 2.11E-04 | 5.94E-04 | 3.68E-04 | 2.64E-03 | 7.95E-04 | 2.70E-03 | 1.37E-03
M4 | 1.95E-03 | 1.77E-03 | 1.23E-03 | 2.94E-04 | 1.00E-01 | 1.08E-01 | 6.25E-02 | 2.83E-02
PIXEL | 3.32E-02 | 3.60E-03 | 3.44E-02 | 2.67E-02 | 2.38E-01 | 1.05E-01 | 2.37E-01 | 1.11E-01
vPINN | 1.15E-02 | 8.57E-03 | 9.07E-03 | 4.68E-03 | 1.88E+00 | 1.11E+00 | 1.73E+00 | 8.41E-01
| LDC
| $A=1$ | $A=5$
u values |  |  |  |  |  |  |  |
PGCAN | 7.75E-03 | 1.32E-02 | 1.32E-03 | 7.54E-04 | 1.44E-02 | 2.94E-03 | 1.34E-02 | 1.20E-02
M4 | 1.17E-02 | 3.55E-03 | 1.11E-02 | 7.64E-03 | 2.75E-01 | 6.72E-02 | 2.62E-01 | 1.94E-01
PIXEL | 2.38E-01 | 2.85E-02 | 2.44E-01 | 1.81E-01 | 7.10E-01 | 3.68E-02 | 7.12E-01 | 6.27E-01
vPINN | 3.53E-02 | 5.26E-03 | 3.59E-02 | 2.54E-02 | 5.92E-01 | 2.48E-02 | 5.99E-01 | 5.40E-01
v values |  |  |  |  |  |  |  |
PGCAN | 1.24E-02 | 2.07E-02 | 2.55E-03 | 1.25E-03 | 1.76E-02 | 4.39E-03 | 1.62E-02 | 1.42E-02
M4 | 2.04E-02 | 5.97E-03 | 1.96E-02 | 1.37E-02 | 3.63E-01 | 8.50E-02 | 3.48E-01 | 2.60E-01
PIXEL | 3.75E-01 | 5.55E-02 | 3.81E-01 | 2.84E-01 | 1.88E+00 | 4.09E-02 | 9.67E-01 | 8.60E-01
vPINN | 6.68E-02 | 9.35E-03 | 6.84E-02 | 4.88E-02 | 1.04E+00 | 2.56E-02 | 1.03E+00 | 1.01E+00
p values |  |  |  |  |  |  |  |
PGCAN | 2.56E-02 | 2.20E-02 | 1.70E-02 | 1.52E-02 | 2.78E-02 | 5.13E-03 | 2.64E-02 | 2.34E-02
M4 | 4.24E-02 | 3.86E-03 | 4.08E-02 | 3.80E-02 | 4.09E-01 | 7.57E-02 | 3.94E-01 | 3.05E-01
PIXEL | 2.99E-01 | 3.80E-02 | 2.91E-01 | 2.62E-01 | 8.37E-01 | 2.89E-02 | 8.33E-01 | 8.03E-01
vPINN | 1.16E-01 | 9.39E-03 | 1.20E-01 | 9.76E-02 | 9.85E-01 | 7.20E-03 | 9.82E-01 | 9.79E-01

Acknowledgements We appreciate the support from the Office of the Naval
Research (award number N000142312485), NASA’s Space Technology Research Grants
Program (award number 80NSSC21K1809), and National Science Foundation (award
number 2211908).

## Declarations

Conflict of interest The authors declare that they have no conflicts of
interest.

Replication of results Our data, models, and codes are accessible via GitHub
upon publication.

## References

  * Washington and Parkinson (2005) Washington, W.M., Parkinson, C.L.: An Introduction to Three-Dimensional Climate Modeling, 2nd edn. (2005)
  * Reynolds (1976) Reynolds, W.C.: Computation of turbulent flows. Annual Review of Fluid Mechanics 8(1), 183–208 (1976)
  * Meng et al. (2022) Meng, C., Seo, S., Cao, D., Griesemer, S., Liu, Y.: When physics meets machine learning: A survey of physics-informed machine learning. arXiv preprint arXiv:2203.16797 (2022)
  * Cuomo et al. (2022) Cuomo, S., Di Cola, V.S., Giampaolo, F., Rozza, G., Raissi, M., Piccialli, F.: Scientific machine learning through physics–informed neural networks: Where we are and what’s next. Journal of Scientific Computing 92(3), 88 (2022)
  * Wang et al. (2022) Wang, H., Planas, R., Chandramowlishwaran, A., Bostanabad, R.: Mosaic flows: A transferable deep learning framework for solving pdes on unseen domains. Computer Methods in Applied Mechanics and Engineering 389, 114424 (2022) https://doi.org/10.1016/j.cma.2021.114424
  * Yousefpour et al. (2023) Yousefpour, A., Foumani, Z.Z., Shishehbor, M., Mora, C., Bostanabad, R.: Gp+: A python library for kernel-based learning via gaussian processes. arXiv preprint arXiv:2312.07694 (2023)
  * Blechschmidt and Ernst (2021) Blechschmidt, J., Ernst, O.G.: Three ways to solve partial differential equations with neural networks—a review. GAMM-Mitteilungen 44(2), 202100006 (2021)
  * Mora et al. (2024) Mora, C., Yousefpour, A., Hosseinmardi, S., Bostanabad, R.: Neural networks with kernel-weighted corrective residuals for solving partial differential equations. arXiv preprint arXiv:2401.03492 (2024)
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., Karniadakis, G.E.: Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics 378, 686–707 (2019)
  * Kang et al. (2023) Kang, N., Lee, B., Hong, Y., Yun, S.-B., Park, E.: Pixel: Physics-informed cell representations for fast and accurate pde solvers. In: Proceedings of the AAAI Conference on Artificial Intelligence, vol. 37, pp. 8186–8194 (2023)
  * Wang et al. (2021) Wang, S., Wang, H., Perdikaris, P.: On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering 384, 113938 (2021)
  * Rahaman et al. (2019) Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F., Bengio, Y., Courville, A.: On the spectral bias of neural networks. In: International Conference on Machine Learning, pp. 5301–5310 (2019). PMLR
  * Sharma et al. (2023) Sharma, P., Evans, L., Tindall, M., Nithiarasu, P.: Stiff-pdes and physics-informed neural networks. Archives of Computational Methods in Engineering, 1–30 (2023)
  * Daw et al. (2022) Daw, A., Bu, J., Wang, S., Perdikaris, P., Karpatne, A.: Mitigating propagation failures in pinns using evolutionary sampling (2022)
  * Lagari et al. (2020) Lagari, P.L., Tsoukalas, L.H., Safarkhani, S., Lagaris, I.E.: Systematic construction of neural forms for solving partial differential equations inside rectangular domains, subject to initial, boundary and interface conditions. International Journal on Artificial Intelligence Tools 29(05), 2050009 (2020)
  * Lagaris et al. (1998) Lagaris, I.E., Likas, A., Fotiadis, D.I.: Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions on neural networks 9(5), 987–1000 (1998)
  * Dong and Ni (2021) Dong, S., Ni, N.: A method for representing periodic functions and enforcing exactly periodic boundary conditions with deep neural networks. Journal of Computational Physics 435, 110242 (2021)
  * McFall and Mahan (2009) McFall, K.S., Mahan, J.R.: Artificial neural network method for solution of boundary value problems with exact satisfaction of arbitrary boundary conditions. IEEE Transactions on Neural Networks 20(8), 1221–1233 (2009)
  * Berg and Nyström (2018) Berg, J., Nyström, K.: A unified deep artificial neural network approach to partial differential equations in complex geometries. Neurocomputing 317, 28–41 (2018)
  * Karnakov et al. (2024) Karnakov, P., Litvinov, S., Koumoutsakos, P.: Solving inverse problems in physics by optimizing a discrete loss: Fast and accurate learning without neural networks. PNAS Nexus (2024) https://doi.org/10.1093/pnasnexus/pgae005
  * Tewari et al. (2022) Tewari, A., Thies, J., Mildenhall, B., Srinivasan, P., Tretschk, E., Yifan, W., Lassner, C., Sitzmann, V., Martin-Brualla, R., Lombardi, S., Simon, T., Theobalt, C., Nießner, M., Barron, J.T., Wetzstein, G., Zollhöfer, M., Golyanik, V.: Advances in neural rendering. Computer Graphics Forum 41(2), 703–735 (2022) https://doi.org/10.1111/cgf.14507
  * Müller et al. (2022) Müller, T., Evans, A., Schied, C., Keller, A.: Instant neural graphics primitives with a multiresolution hash encoding. arXiv preprint arXiv:2201.05989 (2022)
  * Dong et al. (2022) Dong, W., Lao, Y., Kaess, M., Koltun, V.: Ash: A modern framework for parallel spatial hashing in 3d perception. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1–18 (2022) https://doi.org/10.1109/TPAMI.2022.3214347
  * Wang et al. (2021) Wang, S., Teng, Y., Perdikaris, P.: Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing 43(5), 3055–3081 (2021)
  * McClenny and Braga-Neto (2020) McClenny, L., Braga-Neto, U.: Self-adaptive physics-informed neural networks using a soft attention mechanism. arXiv preprint arXiv:2009.04544 (2020)
  * Xu et al. (2019) Xu, Z.-Q.J., Zhang, Y., Luo, T., Xiao, Y., Ma, Z.: Frequency principle: Fourier analysis sheds light on deep neural networks. arXiv preprint arXiv:1901.06523 (2019)
  * Wang et al. (2022) Wang, S., Yu, X., Perdikaris, P.: When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics 449, 110768 (2022)
  * Jacot et al. (2018) Jacot, A., Gabriel, F., Hongler, C.: Neural tangent kernel: Convergence and generalization in neural networks. Advances in neural information processing systems 31 (2018)
  * Liu and Wang (2021) Liu, D., Wang, Y.: A dual-dimer method for training physics-constrained neural networks with minimax architecture. Neural Networks 136, 112–125 (2021)
  * Wight and Zhao (2020) Wight, C.L., Zhao, J.: Solving allen-cahn and cahn-hilliard equations using the adaptive physics informed neural networks. arXiv preprint arXiv:2007.04542 (2020)
  * Fridovich-Keil et al. (2022) Fridovich-Keil, S., Gontijo Lopes, R., Roelofs, R.: Spectral bias in practice: The role of function frequency in generalization. Advances in Neural Information Processing Systems 35, 7368–7382 (2022)
  * Liu et al. (2020) Liu, Z., Cai, W., Xu, Z.-Q.J.: Multi-scale deep neural network (mscalednn) for solving poisson-boltzmann equation in complex domains. arXiv preprint arXiv:2007.11207 (2020)
  * Tancik et al. (2020) Tancik, M., Srinivasan, P., Mildenhall, B., Fridovich-Keil, S., Raghavan, N., Singhal, U., Ramamoorthi, R., Barron, J., Ng, R.: Fourier features let networks learn high frequency functions in low dimensional domains. Advances in Neural Information Processing Systems 33, 7537–7547 (2020)
  * Mildenhall et al. (2021) Mildenhall, B., Srinivasan, P.P., Tancik, M., Barron, J.T., Ramamoorthi, R., Ng, R.: Nerf: Representing scenes as neural radiance fields for view synthesis. Communications of the ACM 65(1), 99–106 (2021)
  * Hong et al. (2022) Hong, Q., Siegel, J.W., Tan, Q., Xu, J.: On the activation function dependence of the spectral bias of neural networks. arXiv preprint arXiv:2208.04924 (2022)
  * Aldirany et al. (2024) Aldirany, Z., Cottereau, R., Laforest, M., Prudhomme, S.: Multi-level neural networks for accurate solutions of boundary-value problems. Computer Methods in Applied Mechanics and Engineering 419, 116666 (2024) https://doi.org/10.1016/j.cma.2023.116666
  * (37) Lindell, D.B., Van Veen, D., Park, J.J., Wetzstein, G.: Bacon: Band-limited coordinate networks for multiscale scene representation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 16252–16262
  * (38) Tancik, M., Casser, V., Yan, X., Pradhan, S., Mildenhall, B., Srinivasan, P.P., Barron, J.T., Kretzschmar, H.: Block-nerf: Scalable large scene neural view synthesis. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 8248–8258
  * (39) Xiangli, Y., Xu, L., Pan, X., Zhao, N., Rao, A., Theobalt, C., Dai, B., Lin, D.: Bungeenerf: Progressive neural radiance field for extreme multi-scale scene rendering. In: European Conference on Computer Vision, pp. 106–122. Springer
  * JAMBON et al. (2023) JAMBON, C., KERBL, B., KOPANAS, G., DIOLATZIS, S., DRETTAKIS, G.: Nerfshop: Interactive editing of neural radiance fields. Proc. ACM Comput. Graph. Interact. Tech 6(1) (2023)
  * (41) Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J., Wojna, Z.: Rethinking the inception architecture for computer vision. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 2818–2826
  * (42) Chollet, F.: Xception: Deep learning with depthwise separable convolutions. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1251–1258
  * Howard et al. (2017) Howard, A.G., Zhu, M., Chen, B., Kalenichenko, D., Wang, W., Weyand, T., Andreetto, M., Adam, H.: Mobilenets: Efficient convolutional neural networks for mobile vision applications. arXiv preprint arXiv:1704.04861 (2017)