11footnotetext:  Corresponding author.

# Limiting Over-Smoothing and Over-Squashing of Graph Message Passing by Deep
Scattering Transforms

Yuanhong Jiang<EMAIL_ADDRESS>
Shanghai Jiao Tong University Dongmian Zou<EMAIL_ADDRESS>
Duke Kunshan University Xiaoqun Zhang<EMAIL_ADDRESS>
Shanghai Jiao Tong University Yu Guang Wang∗<EMAIL_ADDRESS>
Shanghai Jiao Tong University
& University of New South Wales

###### Abstract

Graph neural networks (GNNs) have become pivotal tools for processing graph-
structured data, leveraging the message passing scheme as their core
mechanism. However, traditional GNNs often grapple with issues such as
instability, over-smoothing, and over-squashing, which can degrade performance
and create a trade-off dilemma. In this paper, we introduce a discriminatively
trained, multi-layer Deep Scattering Message Passing (DSMP) neural network
designed to overcome these challenges. By harnessing spectral transformation,
the DSMP model aggregates neighboring nodes with global information, thereby
enhancing the precision and accuracy of graph signal processing. We provide
theoretical proofs demonstrating the DSMP’s effectiveness in mitigating these
issues under specific conditions. Additionally, we support our claims with
empirical evidence and thorough frequency analysis, showcasing the DSMP’s
superior ability to address instability, over-smoothing, and over-squashing.

Keywords: Message Passing Neural Network, Scattering Transform, Over-
smoothing, Over-squashing

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
    1. 2.1 Graphs, Function Space and Spectrum
    2. 2.2 Graph Framelet System
    3. 2.3 Graph Scattering Transform
    4. 2.4 Graph Neural Networks and Message Passing
  3. 3 Related works
    1. 3.1 Message Passing Neural Network
    2. 3.2 Over-smoothness and Over-squashing in MPNN
    3. 3.3 Spectral Methods on GNNs and Scattering Transform
  4. 4 Deep Scattering Message Passing
    1. 4.1 Framelet Scattering Transform
    2. 4.2 Deep Scattering Message Passing Model
    3. 4.3 Limiting Over-smoothing in DSMP
    4. 4.4 Limited Over-squashing in DSMP
    5. 4.5 Stability of DSMP
  5. 5 Numerical Analysis
    1. 5.1 Datasets and Experimental Protocol
    2. 5.2 Evolution of Dirichlet Energy
    3. 5.3 Experiments and Analysis on Over-squashing Issue
    4. 5.4 Stability of DSMP
    5. 5.5 Graph Classification Experiments
    6. 5.6 Computational Cost
  6. 6 Conclusion and Future Work
  7. A Table of notations

##  1 Introduction

Graphs offer an elegant and comprehensive means to describe data’s intricate
topology or structure. A wide array of relational systems and structured
entities have been aptly represented by graphs, encompassing their internal
connections. Notable examples of such applications include recommendation
systems (Wu et al., 2022), prediction of molecular properties (Bumgardner et
al., 2021), analysis of social networks (Min et al., 2021), and particle
system investigations (Li et al., 2009; Wang et al., 2023). The inherent
topology embedded within these graphs enables the establishment of local
communication channels, facilitating the exchange of neighborhood information
(Bruna et al., 2014; Kipf and Welling, 2017; Veličković et al., 2018). In this
paper, we develop a scattering wavelet-based message passing neural network,
which has proved excellent performance in dealing with the commonly prevalent
issues of over-smoothing, over-squashing, and instability in GNNs.

GNNs typically adhere to the _message passing_ paradigm, which is the
foundation of _message passing neural networks_ (MPNNs) (Gilmer et al., 2017).
Most MPNNs extract graph embeddings from the spatial domain and have
demonstrated their utility in diverse real-world tasks. These networks
efficiently capture the complex relationships and interactions within graphs,
making them useful tools for understanding and analyzing structured data.

Despite the remarkable success of MPNNs in various real-world applications,
they still encounter significant challenges with depth scalability and
information bottleneck. One major issue is the over-smoothing phenomenon,
which occurs when node features gradually converge and lose their uniqueness
as the number of layers increases (Li et al., 2018). This convergence not only
diminishes the distinctiveness of node representation but also directly leads
to a decline in performance (Chen et al., 2020). The severity of over-
smoothing is often quantified by analyzing the Dirichlet energy of the
embedded graph data features (Cai and Wang, 2020; Di Giovanni et al., 2022).

Another notable challenge is the over-squashing problem, originally identified
by Alon and Yahav (2021). This arises when node features are unable to
effectively capture information from distant nodes, resulting in a loss of
crucial contextual details. Di Giovanni et al. (2023) delve deep into the
underlying causes of over-squashing, examining various factors such as the
graph topology and neural network architecture, particularly the depth and
width. The severity of the over-squashing of a graph is closely related to the
graph topology. A method to assess the degree of over-squashing in a GNN
involves analyzing the number of non-zero elements within the convolution
matrix. A further approach to examining over-squashing involves a
comprehensive sensitivity analysis conducted on the Jacobian of node features,
as described in Topping et al. (2022), which reveals that negatively curved
edges contribute significantly to this over-squashing issue.

The stability problem is another prevalent issue in GNNs, where the presence
of noisy or perturbed graph signals can adversely affect the performance of
GNN models. Specifically, summation-based MPNNs fail to demonstrate
satisfactory stability and transferability (Yehudai et al., 2021). The process
of feature aggregation in graph representation learning amplifies the negative
impact of irregular entities on their neighborhoods and accumulates noisy
signals during multi-layer propagation, ultimately leading to misleading
embeddings. In fact, the generalization bound associated with this problem is
directly proportional to the largest eigenvalue of the graph Laplacian (Verma
and Zhang, 2019). Therefore, employing spectral graph filters that are robust
to certain structural perturbations becomes a necessary condition for
achieving transferability in graph representation learning.

In this paper, our objective is to overcome the above obstacles by
establishing a message passing framework operating within the spectral domain.
Utilizing wavelet scattering, we can effectively capture long-range patterns
in graphs, thereby mitigating the over-smoothing issue. Additionally, we draw
parallels to the scattering transform to devise a layered message passing
architecture. This neural network performs message passing on graphs with
denser connectivity, addressing the problem of over-squashing. The
meticulously designed DSMP neural network inherits the stability imparted by
scattering wavelet transforms, which conserve energy. This approach empowers
us to effectively capture and manipulate complex graph signal patterns,
advancing the field of GNNs.

Figure 1: Architecture of the proposed scattering message passing (SMP) for
graph signal embedding. In the initial scattering layer $L_{0}$, only the low
pass coefficient is computed with framelet operator
$\bm{{\mathcal{W}}}_{0,J}$. In the following scattering layer $L_{1}$, the
high-pass coefficients are computed with framelet filter
$\bm{{\mathcal{W}}}_{r,l}$ and cascaded with low pass filter
$\bm{{\mathcal{W}}}_{0,J}$ according to the scheme of scattering
transformation. The low pass coefficient of the second scattering layer is
reserved for higher-level scattering transformation. The coefficients are
further processed with the MPNN which is taken to be the GAT model in this
paper. After the aggregation operation, the multi-layer perceptrons (MLPs)
with pointwise nonlinear activation functions are performed to increase the
expressive power of the SMP model.

##  2 Preliminaries

In this section, we introduce some basic notation and properties on graphs,
framelets, and scattering transform, which are based on the works of Chung
(1997); Dong (2017); Wang and Zhuang (2019); Bruna and Mallat (2013); Zou and
Lerman (2020). In accordance with Table A in the appendix, we give a list of
notations used in this paper.

###  2.1 Graphs, Function Space and Spectrum

We consider an undirected attributed graph
$\mathcal{G}=(\mathcal{V},\mathcal{E},\bm{X})$ with a non-empty finite set
$\mathcal{V}$ of vertices, a set
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ of edges between vertices
in $\mathcal{V}$, and the node attributes $\bm{X}\in\mathbb{R}^{N\times d}$
encoding the feature dimensionality of each node, where $d$ denotes the
dimension of features. $\bm{X}$ also denotes a feature matrix of graph signal.
For an undirected graph $\mathcal{G}$, we denote $N:=|\mathcal{V}|$ and
$M:=|\mathcal{E}|$ the numbers of vertices and edges. An edge
$e\in\mathcal{E}$ with vertices $u,v\in\mathcal{V}$ is an unordered pair
denoted by $(u,v)$ or $(v,u)$. The connectivity between $u$ and $v$ is denoted
as $u\sim v$. We build the adjacency matrix $\mathcal{A}\in\mathbb{R}^{N\times
N}$ of $\mathcal{G}$, and $\mathcal{A}$ is a mapping from $\mathcal{E}$ to
$\mathcal{V}\times\mathcal{V}$ by letting

$\mathcal{A}_{uv}=\left\\{\begin{array}[]{cl}1,&u\sim v\\\
0,&\text{otherwise}.\end{array}\right.$

Note that for an undirected graph, the adjacency matrix $\mathcal{A}$ is
symmetric in the sense that $\mathcal{A}(u,v)=\mathcal{A}(v,u),\forall
u,v\in\mathcal{V}$. We denote the degree of a vertex $v\in\mathcal{V}$ by
$\bm{d}(v):=\sum_{p\in\mathcal{V}}\mathcal{A}(v,p)$. The volume of the
$\operatorname{graph}\operatorname{vol}(\mathcal{G}):=\operatorname{vol}(\mathcal{V})=\sum_{v\in\mathcal{V}}\bm{d}(v)=2M$,
which is the sum of degrees of all vertices of $\mathcal{G}$. Given a subset
$V_{0}$ of $\mathcal{V}$, the volume of $V_{0}$ is the sum of degrees of all
nodes in $V_{0}$, that is,
$\operatorname{vol}\left(\mathcal{V}_{0}\right):=\sum_{v\in\mathcal{V}_{0}}d(v)$.

Let
$l_{2}(\mathcal{G}):=l_{2}\left(\mathcal{G},\langle\cdot,\cdot\rangle_{\mathcal{G}}\right)$
be the Hilbert space of vectors $\bm{f}:\mathcal{V}\rightarrow\mathbb{C}$ on
the graph $\mathcal{G}$ equipped with the inner product

$\langle\bm{f},\bm{g}\rangle_{\mathcal{G}}:=\sum_{v\in\mathcal{V}}\bm{f}(v)\overline{\bm{g}(v)},\quad\bm{f},\bm{g}\in
l_{2}(\mathcal{G}),$

where $\bar{\bm{g}}$ is the complex conjugate to $\bm{g}$. The induced norm
$\|\cdot\|_{\mathcal{G}}$ is then given by $\|\bm{f}\|_{\mathcal{G}}:=$
$\sqrt{\langle\bm{f},\bm{f}\rangle_{\mathcal{G}}}$ for $\bm{f}\in
l_{2}(\mathcal{G})$. For simplicity, we shall drop the subscript
$\mathcal{G}$, and use $\langle\cdot,\cdot\rangle$ and $\|\cdot\|$ if no
confusion arises.

Let $\delta_{\ell,\ell^{\prime}}$ be the Kronecker delta function satisfying
$\delta_{\ell,\ell^{\prime}}=1$ if $\ell=\ell^{\prime}$, and
$\delta_{\ell,\ell^{\prime}}=0$ if $\ell\neq\ell^{\prime}$. A finite subset
$\left\\{\bm{u}_{\ell}\right\\}_{\ell=1}^{N}$ of $l_{2}(\mathcal{G})$ is said
an orthonormal basis on $l_{2}(\mathcal{G})$ if
$\left\langle\bm{u}_{\ell},\bm{u}_{\ell^{\prime}}\right\rangle=\delta_{\ell,\ell^{\prime}},1\leq\ell,\ell^{\prime}\leq
N.$ Let $\left\\{\bm{u}_{\ell}\right\\}_{\ell=1}^{N}$ be an orthonormal basis
on $l_{2}(\mathcal{G})$. For $\ell=1,\ldots,N$, let
$\widehat{\bm{f}}_{\ell}:=\left\langle\bm{f},u_{\ell}\right\rangle$ be the
(generalized) Fourier coefficient of degree $\ell$ for $\bm{f}\in
l_{2}(\mathcal{G})$ with respect to $\bm{u}_{\ell}$.

Let
$\widehat{\bm{f}}_{\ell}:=\left(\widehat{\bm{f}}_{1},\ldots,\widehat{\bm{f}}_{N}\right)\in\mathbb{C}^{N}$
be the sequence of the Fourier coefficients for $\bm{f}$. Then, $\bm{f}=$
$\sum_{\ell=1}^{N}\widehat{\bm{f}}_{\ell}\bm{u}_{\ell},\forall\bm{f}\in
l_{2}(\mathcal{G})$, and Parseval’s identity holds:
$\|\bm{f}\|^{2}=\sum_{\ell=1}^{N}\left|\widehat{\bm{f}}_{\ell}\right|^{2},\forall\bm{f}\in
l_{2}(\mathcal{G})$. We say
$\left\\{\left(\bm{u}_{\ell},\lambda_{\ell}\right)\right\\}_{\ell=1}^{N}$ is
an orthonormal eigen-pair for $l_{2}(\mathcal{G})$ if
$\left\\{\bm{u}_{\ell}\right\\}_{\ell=1}^{N}$ is an orthonormal basis for
$l_{2}(\mathcal{G})$ with $u_{1}\equiv 1/\sqrt{N}$ and
$\left\\{\lambda_{\ell}\right\\}_{\ell=1}^{N}\subseteq\mathbb{R}$ is a
nondecreasing sequence of nonnegative numbers satisfying
$0=\lambda_{1}\leq\cdots\leq\lambda_{N}$. A typical example is the eigen-
pairs, that is, the set of all pairs of the eigenvectors and eigenvalues of
the graph Laplacian on $\mathcal{G}$.

We define the Laplacian for graphs $\mathcal{G}$ without loops and multiple
edges. We first consider the unnormalized Laplacian matrix $L$, also known as
the combinatorial graph Laplacian, which is defined as follows:

$L_{uv}=\left\\{\begin{array}[]{cl}d_{v},&u=v\\\ -1,&u\sim v\\\
0,&\text{otherwise},\end{array}\right.$

and let $\mathcal{D}$ the diagonal matrix with the $(v,v)$-th entry having
value $d_{v}$. The Laplacian of $\mathcal{G}$ is defined to be the matrix

$\mathcal{L}_{uv}=\left\\{\begin{array}[]{cl}1,&u=v,d_{v}\neq 0\\\
-\frac{1}{\sqrt{d_{u}d_{v}}},&u\sim v\\\
0,&\text{otherwise}.\end{array}\right.$

We can write

$\mathcal{L}=\mathcal{D}^{-1/2}L\mathcal{D}^{-1/2},$

with the convention $\mathcal{D}^{-1}(v,v)=0$ for $d_{v}=0$. Since
$\mathcal{L}$ is symmetric, its eigenvalues are all real and non-negative. The
set of the eigenvalues is usually called the spectrum of $\mathcal{L}$. The
spectrum of a graph is closely related to its connectivity, and a larger value
of the eigenvalue often indicates denser connectivity. For a graph
$\mathcal{G}$ with $n$ vertices, there are some basic facts about the spectrum
of a graph (Chung, 1997):

  1. (i)

The sum of the eigenvalues

$\sum_{i=1}^{N}\lambda_{i}\leq N$

with equality holding if and only if $\mathcal{G}$ has no isolated vertices.

  2. (ii)

If $\mathcal{G}$ is connected, then $\lambda_{1}>0$. If $\lambda_{i}=0$ and
$\lambda_{i+1}\neq 0$, then $\mathcal{G}$ has exactly $i+1$ connected
components.

  3. (iii)

For all $i\leq N-1$, we have

$\lambda_{i}\leq 2$

with $\lambda_{N-1}=2$ if and only if a connected component of $\mathcal{G}$
is bipartite and nontrivial.

  4. (iv)

The spectrum of a graph is the union of the spectra of its connected
components.

###  2.2 Graph Framelet System

In this section, we will present a concise introduction to the framelet
system, tailored specifically for graph signals. Furthermore, we shall delve
into the fast approximation method, focusing on the Chebyshev approximation
for framelet decomposition and reconstruction operators. This method holds
paramount importance in attaining computationally efficient outcomes. The
original concept was put forward by Dong (2017).

A framelet comprises two basic components: a _filter bank_
$\bm{\eta}:=\\{a;b^{(1)},\dots,b^{(K)}\\}$ and a set of _scaling functions_
$\Psi=\\{\alpha;\beta^{(1)},\dots,\beta^{(K)}\\}$. In this context, $a$
represents the low-pass filter, while $b^{(r)}$, where $r=1,\dots,K$,
represents the $r$-th high-pass filter. These sets of filters extract both the
approximate and detailed information of the input graph signal within the
transformed domain, commonly referred to as the framelet domain.

Different choices of filter masks lead to distinct tight framelet systems
(Dong, 2017). In general, the selected filter bank and its associated scaling
function must satisfy the following relationship:

$\widehat{\alpha}(2\xi)=\widehat{a}(\xi)\widehat{\alpha}(\xi),\quad\widehat{\beta^{(r)}}(2\xi)=\widehat{b^{(r)}}(\xi)\widehat{\alpha}(\xi),$
(1)

where $r=1,\ldots,K$ and $\xi\in\mathbb{R}$. As a common implementation, we
consider the _Haar-type_ filter with a single high-pass filter. Specifically,
for $x\in\mathbb{R}$, we implement the _Haar-type_ filter with one high-pass,
i.e., $K=1$. For $x\in\mathbb{R}$, it defines

$\widehat{a}(x)=\cos(x/2)\mbox{~{}~{}and~{}~{}}\widehat{b^{(1)}}(x)=\sin(x/2),$

with scaling functions

$\widehat{\alpha}(\xi)=\frac{\sin(\xi/2)}{\xi/2},\quad\widehat{\beta^{(1)}}(\xi)=\sqrt{1-\left(\frac{\sin(\xi/2)}{\xi/2}\right)^{2}}.$

By utilizing the defined filter bank and scaling functions, we can establish
the undecimated framelet basis for a graph signal with the eigenpairs
$\\{(\lambda_{\ell},{\bm{u}}_{\ell})\\}_{\ell=1}^{N}$ of the graph Laplacian
$L$. For each level $l=1,\dots,J$, the undecimated framelets corresponding to
node $p$ are defined as follows:

$\begin{array}[]{l}\bm{\varphi}_{l,p}(v):=\sum_{\ell=1}^{N}\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\overline{{\bm{u}}_{\ell}(p)}{\bm{u}}_{\ell}(v),v\in{\mathcal{V}}\\\\[5.69054pt]
\bm{\psi}_{l,p}^{(r)}(v):=\sum_{\ell=1}^{N}\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\overline{{\bm{u}}_{\ell}(p)}{\bm{u}}_{\ell}(v),v\in{\mathcal{V}},\quad
r=1,\ldots,K.\end{array}$ (2)

Here, $\bm{\varphi}_{l,p}(v)$ and $\bm{\psi}_{l,p}^{(r)}(v)$ represent the
low-pass and the $r$-th high-pass framelet basis, respectively.

#### Framelet Decomposition

Next, we introduce the _framelet decomposition operator_ denoted as
$\bm{{\mathcal{W}}}$ for the tight framelet system, which is comprised of
individual operators ${\bm{{\mathcal{W}}}_{k,l}}$ which we will specify
shortly. In Section 2.3, we will utilize this operator to construct the SMP
for the graph signal.

The framelet decomposition operator $\bm{{\mathcal{W}}}_{k,l}$ consists of
orthonormal bases for
$(k,l)\in\\{(0,J)\\}\cup\\{(1,1),\dots,(1,J),\dots(K,1),\dots,(K,J)\\}$. It
transforms a given graph signal ${\bm{X}}$ into a set of _framelet
coefficients_ at multiple scales and levels. These coefficients represent the
spectral representation of the graph signal ${\bm{X}}$ in the transformed
domain. Specifically, $\bm{{\mathcal{W}}}_{0,J}$ includes
$\bm{\varphi}_{J,p}$, where $p\in{\mathcal{V}}$, and constructs the low-pass
framelet coefficients $\bm{{\mathcal{W}}}_{0,J}{\bm{X}}$, preserving the
approximate information in ${\bm{X}}$. These coefficients provide a smooth
representation that captures the global trend of ${\bm{X}}$. On the other
hand, the high-pass coefficients $\bm{{\mathcal{W}}}_{r,l}{\bm{X}}$, where
$\bm{{\mathcal{W}}}_{r,l}={\bm{\psi}_{l,p}^{(r)},p\in{\mathcal{V}}}$, capture
detailed information at scale $r$ and level $l$. They reveal localized
patterns or noise in the signal. A larger scale $k$ contains more localized
information with smaller energy.

The framelet coefficients can be directly projected with
$\langle\bm{\varphi}_{l,p},{\bm{X}}\rangle$ and
$\langle\bm{\psi}_{l,p}^{(r)},{\bm{X}}\rangle$ for node $p$ at scale level
$l$. For instance, if we perform eigendecomposition on a normalized graph
Laplacian $\mathcal{L}$, resulting in
${\bm{U}}=[{\bm{u}}_{1},\dots,{\bm{u}}_{N}]\in\mathbb{R}^{N\times N}$ as the
eigenvectors and $\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{N})$
as the eigenvalues, we can define the corresponding filtered diagonal matrices
with low-pass and high-pass filters as follows:

$\displaystyle\widehat{\alpha}\left(\frac{\Lambda}{2}\right)$
$\displaystyle=\operatorname{diag}\left(\widehat{\alpha}\left(\frac{\lambda_{1}}{2}\right),\dots,\widehat{\alpha}\left(\frac{\lambda_{n}}{2}\right)\right),\
\widehat{\beta^{(r)}}\left(\frac{\Lambda}{2^{l}}\right)$
$\displaystyle=\operatorname{diag}\left(\widehat{\beta^{(r)}}\left(\frac{\lambda_{1}}{2^{l}}\right),\dots,\widehat{\beta^{(r)}}\left(\frac{\lambda_{n}}{2^{l}}\right)\right).$

The associated framelet coefficients for the low-pass and the $r$-th high-pass
signals are given by:

$\displaystyle\bm{{\mathcal{W}}}_{0,J}{\bm{X}}={\bm{U}}\widehat{\alpha}\left(\frac{\Lambda}{2}\right){\bm{U}}^{\top}{\bm{X}},\
$
$\displaystyle\bm{{\mathcal{W}}}_{r,l}{\bm{X}}={\bm{U}}\widehat{\beta^{(r)}}\left(\frac{\Lambda}{2^{l}}\right){\bm{U}}^{\top}{\bm{X}},\quad\forall
l=1,\dots,J.$ (3)

#### Approximate Framelet Transformation

To achieve an efficient framelet decomposition, we employ two strategies: a
recursive formulation on filter matrices and the utilization of Chebyshev-
approximated eigenvectors.

Assuming a filter bank that satisfies (1), we can implement the decomposition
recursively as follows:

$\displaystyle\bm{{\mathcal{W}}}_{r,1}{\bm{X}}={\bm{U}}\widehat{\beta^{(r)}}\left(2^{-R}\Lambda\right){\bm{U}}^{\top}{\bm{X}},$
(4)

for the first level ($l=1$), and

$\displaystyle\bm{{\mathcal{W}}}_{r,l}{\bm{X}}$
$\displaystyle={\bm{U}}\widehat{\beta^{(r)}}\left(2^{-R+l-1}\Lambda\right)\widehat{\alpha}\left(2^{-R+l-2}\Lambda\right)\dots\widehat{\alpha}\left(2^{-R}\Lambda\right){\bm{U}}^{\top}{\bm{X}}\
$ $\displaystyle\quad\forall l=2,\dots,J,$ (5)

where the real-value dilation scale $R$ satisfies $\lambda_{\max}\leq
2^{R}\pi$. In this definition, it is required for the finest scale
$\frac{1}{2^{K+J}}$, that guarantees
$\frac{\lambda_{\ell}}{2^{K+J-\ell}}\in(0,\pi),\forall\ell=1,2,\dots,J$.

Furthermore, we utilize an $r$-order Chebyshev polynomial approximation to
achieve an efficient framelet decomposition. This approach avoids the
eigendecomposition of the graph Laplacian, which can be time-consuming for
large graphs. Let the filters $\alpha$ and $\beta^{(r)}$ be approximated by
${\mathcal{T}}_{0}$ and ${\mathcal{T}}_{r}$, respectively, using $r$-order
approximation. The framelet decomposition operator $\bm{{\mathcal{W}}}_{r,l}$
is then approximated as follows:

$\bm{{\mathcal{W}}}_{r,l}^{\natural}=\begin{cases}{\mathcal{T}}_{r}\left(2^{-R}\mathcal{L}\right),&l=1,\\\\[2.84526pt]
{\mathcal{T}}_{r}\left(2^{-R+l-1}\mathcal{L}\right){\mathcal{T}}_{0}\left(2^{-R+l-2}\mathcal{L}\right)\dots\mathcal{T}_{0}\left(2^{-R}\mathcal{L}\right),&l=2,\dots,J.\end{cases}$
(6)

The above formulation can be efficiently computed using the Laplacian matrix
$\mathcal{L}$, and the resulting transformation matrices exhibit a close
relationship with the adjacency matrix $\mathcal{A}$.

In a spectral-based graph convolution, graph signals are transformed to a set
of coefficients $\hat{\bm{X}}=\bm{{\mathcal{W}}}\bm{X}^{\text{in }}$ in
frequency channels to learn useful graph representation. With $K$ high-pass
filters,
$\hat{\bm{X}}=\left\\{\bm{{\mathcal{W}}}_{0,J}\bm{X};\bm{{\mathcal{W}}}_{r,l}\bm{X}:r=1,\ldots,K,l=1,\ldots,J\right\\}$
constitutes of a low-pass coefficient matrix $\bm{{\mathcal{W}}}_{0,J}\bm{X}$
and $K$ high-pass coefficient matrices $\bm{{\mathcal{W}}}_{r,l}\bm{X}$ at $J$
scale levels. The framelet coefficient at node $i$ gathers information from
its neighbors in the framelet domain of the same channel. Since the original
(undecimated) framelet transforms involve integrals, we employ Chebyshev
polynomial approximation for $\bm{{\mathcal{W}}}_{r,l}$ (including
$\bm{{\mathcal{W}}}_{0,J}$), and obtain the corresponding approximated
framelet transformation matrices $\bm{{\mathcal{W}}}_{r,l}^{\natural}$ as
products of Chebyshev polynomials of $A$. Let $m$ denote the highest order of
the Chebyshev polynomial used. The aggregation of framelet coefficients for
the neighborhood of node $i$ extends up to the $m$-th multi-hop
$\mathcal{N}^{m}(i)$.

###  2.3 Graph Scattering Transform

The scattering transform is a multi-scale signal analysis method that provides
stable representations using framelets. The graph scattering transform extends
this method to analyze graph signals. A scattering transform can be
constructed by cascading three building blocks: a collection of paths
containing the ordering of the scales of framelets, a pointwise modulus
activation function, and a low-pass framelet operator. Following the same
framework, given a graph $\mathcal{G}$ and input feature matrix ${\bm{X}}\in
l^{2}(\mathcal{G})$, we define an analogous graph scattering transformation
$\Phi_{\mathcal{G}}({\bm{X}})$ by cascading three building blocks: a
collection of paths
$\mathbf{p}=\left\\{(1,1),\dots,(1,J),\dots,(K,J)\right\\}^{m-1}$, each
containing the high-pass components of framelet decomposition operators
$\bm{{\mathcal{W}}}_{r,l},\;1\leq r\leq K,\;1\leq l\leq J,$ where $m$ is the
depth of the transformation; a pointwise absolute value function
$\hat{\sigma}$; and a low-pass component $\bm{{\mathcal{W}}}_{0,J}$. More
specifically, the first layer of scattering transformation can be formulated
as

$\phi_{1}(\mathcal{G},{\bm{X}})=\bm{{\mathcal{W}}}_{0,J}{\bm{X}},$ (7)

followed by the second layer

$\phi_{2}(\mathcal{G},{\bm{X}})=\bm{{\mathcal{W}}}_{0,J}\hat{\sigma}(\bm{{\mathcal{W}}}{\bm{X}})=\left\\{\bm{{\mathcal{W}}}_{0,J}\sigma(\bm{{\mathcal{W}}}_{r,l}{\bm{X}})\right\\}_{{1\leq
r\leq K},{1\leq l\leq J}},$ (8)

and the propagation on the $m$-th layer reads

$\phi_{m}(\mathcal{G},{\bm{X}})=\bm{{\mathcal{W}}}_{0,J}\underbrace{\hat{\sigma}(\bm{{\mathcal{W}}}\cdots\hat{\sigma}(\bm{{\mathcal{W}}}{\bm{X}}))}_{(m-1)*{\mathcal{W}}}=\left\\{\bm{{\mathcal{W}}}_{0,J}(\hat{\sigma}(\bm{{\mathcal{W}}}_{r,l}\cdots\hat{\sigma}(\bm{{\mathcal{W}}}_{r,l}{\bm{X}}))\right\\}_{{1\leq
r\leq K},{1\leq l\leq J}}.$ (9)

The representation obtained by concatenating output features from $m$ layers
of such transformation is thus

$\Phi_{\mathcal{G}}({\bm{X}})=\left\\{\phi_{1}(\mathcal{G},{\bm{X}}),\phi_{2}(\mathcal{G},{\bm{X}}),\ldots,\phi_{m-1}(\mathcal{G},{\bm{X}})\right\\}.$
(10)

###  2.4 Graph Neural Networks and Message Passing

Given a graph with adjacency matrix $\mathcal{A}$ and input feature matrix
$\bm{X}^{\text{in }}$, a graph convolution produces a matrix representation
$\bm{X}^{\text{out }}$. Let$\bm{X}_{i}^{\text{in}}$ denote the $i$-th row of
$\bm{X}^{\text{in}}$, that is, the feature for node $i$. At a specific layer,
the propagation for the $i$-th node reads

$\bm{X}_{i}^{\text{out }}=\Gamma\left(\bm{X}_{i}^{\text{in
}},\square_{j\in\mathcal{N}(i)}\phi\left(\bm{X}_{i}^{\text{in
}},\bm{X}_{j}^{\text{in }},\mathcal{A}_{ij}\right)\right),$ (11)

where $\square(\cdot)$ is a differentiable and permutation invariant
aggregation function, such as summation, average, or maximization. The set
$\mathcal{N}(i)$ includes the $i$-th node and its 1-hop neighbors. Both
$\Gamma(\cdot)$ and $\phi(\cdot)$ are MLPs.

#### Over-smoothing

The over-smoothing problem commonly exists in GNN models with message passing
schemes, which makes it difficult to achieve optimal performance and to
implement in real-world applications. The proposed DSMP framework is
theoretically proved to solve the over-smoothing problem in the sense of
Dirichlet energy. The Dirichlet energy quantifies the amount of variation or
oscillation in graph signal values across nodes.

The Dirichlet energy of graph $\mathcal{G}$ is defined as:

$E_{\mathcal{G}}({\bm{X}})=\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\mathcal{A}_{ij}({\bm{X}}_{i}-{\bm{X}}_{j})^{2},$
(12)

where $\mathcal{A}_{ij}$ represents the connection strength between nodes $i$
and $j$ in the adjacency matrix. The Dirichlet energy measures the squared
differences between connected node signal values, weighted by the adjacency
matrix. A higher Dirichlet energy indicates greater signal variation or less
smoothness, while a lower value indicates smoother signal behavior.

#### Over-squashing

According to the recent theory on over-squashing of MPNN (Di Giovanni et al.,
2023), the graph structure plays a pivotal role in the quality of signal
propagation across nodes. A classic measure of ”bottlenecked-ness” of a graph
structure is the Cheeger constant (Chung, 1997)

$h(\mathcal{G}):=\min_{\varnothing\neq\mathcal{S}\subset\mathcal{V}}\frac{|\partial\mathcal{S}|}{\min(\operatorname{vol}\mathcal{S},\operatorname{vol}\overline{\mathcal{S}})},$

where
$\operatorname{vol}\mathcal{S}=\sum_{v\in\mathcal{S}}d_{v},\overline{\mathcal{S}}=\mathcal{V}\backslash\mathcal{S}$
and
$\partial\mathcal{S}=\\{(u,v):u\in\mathcal{S},v\in\overline{\mathcal{S}},(u,v)\in\mathcal{E}\\}$
is the edge boundary of $\mathcal{S}\subset\mathcal{V}$. This definition
considers the connectivity between two complementary subsets of nodes in terms
of the number of connecting edges and their volume, measured by the sum of the
degrees of the nodes. A graph with two tightly connected communities ( $S$ and
$\mathcal{V}\backslash\mathcal{S}$) and few inter-community edges has a small
Cheeger constant. Intuitively, a large value of $h(\mathcal{G})$ indicates
that the graph has no significant bottlenecks, meaning that each part of the
graph is well-connected to the rest through a significant fraction of its
edges. A higher $h(\mathcal{G})$ implies more connections and less over-
squashing on the graph scale.

In general, computing the exact value of $h(\mathcal{G})$ is a computationally
challenging task. However, it is possible to provide lower and upper bounds
for the Cheeger constant in terms of the spectral gap. As mentioned in Section
2.1, the spectral gap of a graph is defined as the difference
$\lambda_{2}-\lambda_{1}$, where $\lambda_{1}=0$. The discrete Cheeger
inequality (Cheeger, 1970; Alon and Milman, 1984) establishes a relationship
between the spectral gap of $\mathcal{G}$ and its Cheeger constant:

$\frac{\lambda_{2}}{2}\leq h(\mathcal{G})\leq\sqrt{2\lambda_{2}}.$ (13)

The discrete Cheeger inequality provides us with the insight that the spectral
gap is non-zero if and only if $h(\mathcal{G})$ is non-zero. Consequently, a
graph $\mathcal{G}$ is considered connected if and only if $\lambda_{2}>0$. As
the spectral gap can be efficiently computed, we utilize it as a measure of
graph bottlenecking, serving as a proxy for the Cheeger constant. Thus, the
Cheeger constant is closely related to the over-squashing issue in GNNs.

##  3 Related works

###  3.1 Message Passing Neural Network

MPNN (Gilmer et al., 2017) has become a computational paradigm for GNN models
to learn from graph-structured data. In one message-passing step, each node of
the graph updates its state by aggregating messages flowing from its direct
neighbors. Different GNN layers can be derived by varying the specific
aggregation operation. For instance, the graph convolutional network (GCN)
(Kipf and Welling, 2017) adds up the degree-normalized node attributes from
neighbors (including itself). The graph attention network (GAT) (Veličković et
al., 2018) and other variants (Xie et al., 2020; Brody et al., 2022) refine
the aggregation weights by utilizing the attention mechanism. The expressive
power of MPNNs has been verified by the Weisfeiler-Leman graph isomorphism
test (Xu et al., 2018; Morris et al., 2019; Bodnar et al., 2021), based on
which the graph isomorphism network (GIN) has been proposed to become more
powerful in expressivity.

###  3.2 Over-smoothness and Over-squashing in MPNN

The over-smoothing issue was initially identified by (Li et al., 2018), who
observed that the output embeddings of different nodes tend to converge to
similar vectors after going through multiple layers of signal smoothing
operations. This issue has been a focal point in the study of GNNs and poses a
challenge to the performance of deep GNNs. On the other hand, using only a
small number of message-passing operations limits the expressivity of the
model. When layers are deeply stacked, the embeddings of connected nodes tend
to converge, becoming indistinguishable.

Another commonly investigated issue is the over-squashing problem. Such
phenomenon has been explored in relation to graph curvature, based on which
SDRF (Topping et al., 2022) and BORF (Nguyen et al., 2023) proposed graph
rewiring methods to address the over-squashing problem. FoSR (Karhadkar et
al., 2023) proposed a rewiring strategy based on spectral expansion theory
(Banerjee et al., 2022), aiming to rewire the graph and obtain a larger
spectral gap. Different from rewiring the input graph, graph transformer
models such as GraphTrans (Wu et al., 2021), Graphormer (Ying et al., 2021),
GPS (Rampášek et al., 2022), and U2GNN (Nguyen et al., 2022) took a special
approach to solve the over-squashing issue. These models implicitly learn on a
fully connected graph inspired by the Transformer architecture (Vaswani et
al., 2017). Using the partial derivative between pairs of graph vertices, (Di
Giovanni et al., 2023) showed different perspective to analyze the over-
squashing problem from the graph topology, depth of message passing scheme,
and number of neurons, based on which more efficient MPNNs could be designed
to alleviate the over-squashing problem. One straightforward solution involves
increasing the number of message passing layers, yet this strategy can
potentially exacerbate the over-smoothing issue and amplify the risk of
vanishing gradients. Alternatively, widening the networks can alleviate the
over-squashing problem, but this may simultaneously introduce instability
concerns.

###  3.3 Spectral Methods on GNNs and Scattering Transform

In recent years, spectral-based graph convolutions have shown promising
performance in transferring trained graph convolutions across different
graphs, indicating their transferability and generalizability (Levie et al.,
2019; Maskey et al., 2022). Furthermore, the stability of spectral-based
methods in the face of input graph perturbations has been emphasized in
various studies (Gama et al., 2020; Ruiz et al., 2021). One notable work by
Zheng et al. (2021) extracts multi-scale features from graph data and
integrates the framelet coefficients with a learnable GNN model, achieving
good results in various graph-related tasks.

Scattering methods, inspired by the layered architecture of framelet
transformations, have been used to emulate deep neural networks and serve as
powerful feature extractors that capture multi-scale group invariant features
from data (Mallat, 2012). The scattering transform has also been employed to
analyze and interpret convolutional neural networks (Bruna and Mallat, 2013).
Similarly to the Euclidean scattering transform, the graph scattering
transform applies a set of framelet filters to graph signal(Gama et al.,
2019b; Gao et al., 2019).

In recent years, prominent researchers including (Zou and Lerman, 2020, 2019;
Gama et al., 2019a, b) have successfully extended the scattering transform to
graph data. By leveraging graph framelets in both spectral and spatial
domains, these extensions have generated stable and equivariant feature
representations for instances and their intricate relationships. Furthermore,
efforts have been made to integrate geometric information with graph
scattering, leading to enhanced methods for graph analysis (Gao et al., 2019;
Chew et al., 2022; Xiong et al., 2023).

Despite the utilization of predefined framelets in scattering transforms, it
is feasible to construct learnable filters within this framework. Such a new
approach was initially implemented in the Euclidean domain (Oyallon et al.,
2017; Gauthier et al., 2022) and more recently in the graph domain (Min et
al., 2020; Zhang et al., 2022). More recently, Zhang et al. (2022) employed
diffusion wavelet techniques to develop a hierarchical diffusion scattering
(HDS) network, seamlessly integrating it with GCN (HDS-GCN) and GAT (HDS-GAT)
models. These hierarchical scattering features significantly enhance the
performance of GNN models, further broadening their capabilities and scope.

##  4 Deep Scattering Message Passing

In this section, we introduce a novel graph embedding neural network called
the DSMP(Deep Scattering Message Passing) model, designed to enhance the
learning capability of scattering transform. The foundational component is the
SMP module, which enables the scalability of the DSMP model. The main obstacle
to scalability for common GNN models is the over-smoothing problem, which is
theoretically proved to be circumvented by the proposed model. We also prove
that the newly designed DSMP can alleviate the over-squashing of message
passing. We also prove the stability of the DSMP.

###  4.1 Framelet Scattering Transform

We utilize framelets to construct the graph scattering transform, which we
refer to as the _framelet scattering transform_. The detailed introduction and
efficient computation of framelet decomposition operators
$\bm{{\mathcal{W}}}_{r,l}$ are provided in section 2.2.

In Figure 1, we illustrate a two-layer scattering propagation process using
the framelet transformation matrices, where the pointwise activation function
$\sigma$ is moved outside. These framelet transformations provide a weighted
aggregation scheme for graph features, which can be regarded as weighted
adjacency matrices for message passing and also rewired graph topology as the
level of framelet transformation $\ell\geq 2$. The low pass coefficient is
directly computed with framelet operator $\bm{{\mathcal{W}}}_{0,J}$ in the
initial scattering layer. When the scattering layer increases, the low pass
coefficient of the following scattering layer is reserved for higher-level
scattering transformation.

###  4.2 Deep Scattering Message Passing Model

The scattering transform is a multi-scale feature extraction method without
learnable parameters, which entangles its learning ability in a data-driven
paradigm. Using message passing scheme, we propose the SMP framework to enable
scattering transform learning ability. The DSMP can be constructed with multi-
layer SMP modules. For the update of graph node feature $\bm{X}^{(t)}$ at
layer $t$, we define the following diffusion propagation process:

$\bm{X}^{(t)}=\bm{X}^{(t-1)}+\bm{Z}^{(t)}.$ (14)

Below, we define one-step propagation of SMP with a depth of 2:

$\displaystyle\bm{X}_{j}^{h_{r,l}}$
$\displaystyle=\sum_{p\in\mathcal{N}(j)}\bm{{\mathcal{W}}}_{r,l}^{p,j}\bm{X}_{p}^{(t-1)},$
(15) $\displaystyle\bm{Z}_{i}^{(t)}$
$\displaystyle=\Gamma\left(\sum_{j\in\mathcal{N}(i)}\bm{\theta}_{0}\bm{{\mathcal{W}}}_{0,J}^{i,j}\bm{X}_{j}^{(t-1)}+\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{\theta}_{r,l}\bm{{\mathcal{W}}}_{0,J}^{i,j}\bm{X}_{j}^{h_{r,l}}\right),$

where the first equation denotes the high-pass coefficients of node $j$
extracted by transformation matrix $\bm{{\mathcal{W}}}_{r,l}$. The message
passing process specifically aggregates neighboring nodes from
$\mathcal{N}(j)$, where the adjacency is determined by the non-zero terms of
$\bm{{\mathcal{W}}}_{r,l}$. In the case of the framelet transformation matrix
$\bm{{\mathcal{W}}}_{r,l}$, if the term in the $p$-th row and $j$-th column is
non-zero, denoted as $\bm{{\mathcal{W}}}_{r,l}^{p,j}\neq 0$, it will influence
the aggregation process. To adjust for different embeddings, $\bm{\theta}_{0}$
and $\bm{\theta}_{r,l}$ are introduced as learnable weighting parameters.

The propagation scheme described above adheres to the fundamental principle of
feature aggregation in message passing. For simplicity on notation, we give
the matrix formulation as

$\displaystyle\bm{X}^{(t)}=\bm{X}^{(t-1)}+\bm{Z}^{(t)},$ (16)
$\displaystyle\bm{Z}^{(t)}=\Gamma\left(\bm{\theta}_{0}\bm{{\mathcal{W}}}_{0,J}\bm{X}^{(t-1)}+\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{\theta}_{r,l}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}\bm{X}^{(t-1)}\right).$

As shown in (16), the first term corresponds to low pass framelet coefficients
weighted by the learnable parameter $\bm{\theta}_{0}$ from the GAT model. The
other part consists of high-pass framelet coefficients multiplied with low-
pass transformation matrix $\bm{{\mathcal{W}}}_{0,J}$, which is designed
following the architecture of two layers of scattering transform.

Following aggregation and weighting in the spectral domain, the MLPs are
applied to increase the expressive power of the SMP model. For the sake of
notational simplicity, we represent the MLPs, denoted as $\Gamma$, with the
parameter matrix $\mathbf{P}$ and nonlinear activation function ReLU
(Krizhevsky et al., 2017) $\sigma(\cdot)$. Specifically, for an input signal
$\bm{X}$, $\Gamma(\bm{X})=\sigma(\mathbf{P}\bm{X})$. Worth noting is that the
enhancement of the traditional pointwise modulus activation function is used
in scattering transformations. This improvement has been achieved by
incorporating the pointwise ReLU function within the MLPs. Importantly, this
modification maintains and indeed enhances the integrity of the underlying
scattering transformation design. Specifically, the ReLU function fulfills a
role analogous to $\sigma(\cdot)$, contributing to the model’s nonlinear
approximation power. As illustrated in Figure 2, the DSMP model is constructed
by cascading multiple SMP modules to expand the scale of this deep learning
model.

Figure 2: Overall architecture of the DSMP scheme. Following the propagation
process of (16), the DSMP is constructed by cascading multiple SMP modules.

###  4.3 Limiting Over-smoothing in DSMP

We utilize the Dirichlet energy as a measure of the level of over-smoothness
in message passing models. Firstly, we establish the energy conservation
property of the framelet transformation. Subsequently, we provide proof
demonstrating how the DSMP model effectively constrains the graph embedding.

######  Proposition 1 (Energy conservation)

Let ${\bm{X}}$ be a feature matrix from the graph signal $\mathcal{G}$ and
$\left\\{\bm{{\mathcal{W}}}_{0,0};\bm{{\mathcal{W}}}_{r,l}:r=1,\ldots,K,l=1,\ldots,J\right\\}$
be the framelet transformation matrix. Then,

$E_{\mathcal{G}}({\bm{X}})=E_{\mathcal{G}}(\bm{{\mathcal{W}}}_{0,J}{\bm{X}})+\sum_{r=1}^{K}\sum_{l=1}^{J}E_{\mathcal{G}}(\bm{{\mathcal{W}}}_{r,l}{\bm{X}}),$
(17)

where framelet transformation breaks the total energy
$E_{\mathcal{G}}({\bm{X}})$ into low pass energy
$E_{\mathcal{G}}(\bm{{\mathcal{W}}}_{0,J}{\bm{X}})$ and high-pass energy
$E_{\mathcal{G}}(\bm{{\mathcal{W}}}_{r,l}{\bm{X}})$.

Proof  By the definition of Dirichlet energy, we have
$E_{\mathcal{G}}({\bm{X}})=\mathrm{tr}({\bm{X}}^{\top}\mathcal{L}{\bm{X}})$
and $\mathcal{L}={\bm{U}}\Lambda{\bm{U}}^{\top}$. Then, according to (3), we
have

$E_{\mathcal{G}}(\bm{{\mathcal{W}}}_{0,J}{\bm{X}})=\mathrm{tr}\left({\bm{X}}^{\top}\bm{{\mathcal{W}}}_{0,J}^{\top}\mathcal{L}\bm{{\mathcal{W}}}_{0,J}{\bm{X}}\right)=\mathrm{tr}\left({\bm{X}}^{\top}{\bm{U}}\widehat{\alpha}\left(\frac{\Lambda}{2}\right)^{2}\Lambda{\bm{U}}^{\top}{\bm{X}}\right),$

and $\forall l=1,\dots,J,$

$E_{\mathcal{G}}(\bm{{\mathcal{W}}}_{r,l}{\bm{X}})=\mathrm{tr}\left({\bm{X}}^{\top}\bm{{\mathcal{W}}}_{0,J}^{\top}\mathcal{L}\bm{{\mathcal{W}}}_{r,l}{\bm{X}}\right)=\mathrm{tr}\left({\bm{X}}^{\top}{\bm{U}}\widehat{\beta^{(r)}}\left(\frac{\Lambda}{2^{l}}\right)^{2}\Lambda{\bm{U}}^{\top}{\bm{X}}\right).$
(18)

Since $\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{N})$ and the
tightness of Framelet system requires the functions in $\Psi$ satisfy

$\displaystyle
1=\left|\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2^{J}}\right)\right|^{2}+\sum_{r=1}^{K}\left|\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{J}}\right)\right|^{2}$
$\displaystyle\quad\forall\ell=1,\ldots,N,$ (19)
$\displaystyle\left|\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\right|^{2}=\left|\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\right|^{2}+\sum_{r=1}^{K}\left|\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\right|^{2}$
$\displaystyle\quad\forall\begin{array}[]{l}\ell=1,\ldots,N,\\\
l=1,\ldots,J-1,\end{array}$ (22)

we thus have

$\widehat{\alpha}\left(\frac{\Lambda}{2}\right)^{2}+\sum_{r=1}^{K}\sum_{l=1}^{J}\widehat{\beta^{(r)}}\left(\frac{\Lambda}{2^{l}}\right)^{2}={\bm{I}}.$
(23)

Therefore, combining (18) and (23), we obtain that

$\displaystyle
E_{\mathcal{G}}(\bm{{\mathcal{W}}}_{0,J}{\bm{X}})+\sum_{r=1}^{K}\sum_{l=1}^{J}E_{\mathcal{G}}(\bm{{\mathcal{W}}}_{r,l}{\bm{X}})$
$\displaystyle=$
$\displaystyle\mathrm{tr}\left({\bm{X}}^{\top}{\bm{U}}\widehat{\alpha}\left(\frac{\Lambda}{2}\right)^{2}\Lambda{\bm{U}}^{\top}{\bm{X}}\right)+\sum_{r=1}^{K}\sum_{l=1}^{J}\mathrm{tr}\left({\bm{X}}^{\top}{\bm{U}}\widehat{\beta^{(r)}}\left(\frac{\Lambda}{2^{l}}\right)^{2}\Lambda{\bm{U}}^{\top}{\bm{X}}\right)$
$\displaystyle=$
$\displaystyle\mathrm{tr}\left({\bm{X}}^{\top}{\bm{U}}\left(\widehat{\alpha}\left(\frac{\Lambda}{2}\right)^{2}+\sum_{r=1}^{K}\sum_{l=1}^{J}\widehat{\beta^{(r)}}\left(\frac{\Lambda}{2^{l}}\right)^{2}\right)\Lambda{\bm{U}}^{\top}{\bm{X}}\right)$
$\displaystyle=$
$\displaystyle\mathrm{tr}\left({\bm{X}}^{\top}{\bm{U}}\Lambda{\bm{U}}^{\top}{\bm{X}}\right)$
$\displaystyle=$ $\displaystyle E_{\mathcal{G}}({\bm{X}}),$

thus completing the proof.

We need the following property of the positive semi-definite matrix.

###### Lemma 2

(Zhou, 2014) Let $\mathcal{A}_{k}\in\mathbb{R}^{N\times N}$ be positive semi-
definite $(k=1,2,\cdots,n)$. Then,

$\mathrm{tr}\left(\mathcal{A}_{1}\mathcal{A}_{2}\cdots\mathcal{A}_{n}\right)\leq\sum_{i=1}^{d}\mu_{i}\left(\mathcal{A}_{1}\mathcal{A}_{2}\cdots\mathcal{A}_{n}\right)\leq\mathrm{tr}\left(\mathcal{A}_{1}\right)\mathrm{tr}\left(\mathcal{A}_{2}\right)\cdots\mathrm{tr}\left(\mathcal{A}_{n}\right),$

where $\mu_{i}(\mathcal{A}_{k})$ represents the $i$-th largest eigenvalue of
matrix $\mathcal{A}_{k}$.

The following theorem shows that the Dirichlet energy of the feature at the
$t$-th layer $\bm{X}^{(t)}$ is not less than that of the $(t-1)$-th layer
$\bm{X}^{(t-1)}$ in DSMP model. In fact, $\bm{X}^{(t)}$ and $\bm{X}^{(t-1)}$
are equivalent up to some constant, and the constant with respect to the upper
bound depends jointly on the level of the framelet decomposition and the
learnable parameters’ bound.

Although the asymptotic behavior of the output $X^{(L)}$ of the GCN as
$L\rightarrow\infty$ is investigated with over-smoothing property, to
facilitate the analysis for our model MPNN, we also sketch the proof for
completeness here.

###### Lemma 3

(Oono and Suzuki, 2019) For any $\bm{X}\in\mathbb{R}^{n\times d}$ and
activation function ReLU $\sigma(\cdot)$, we have

  * •

$E_{\mathcal{G}}\left(\bm{W}_{l}\bm{X}\right)\leq\mu_{\max}^{2}E_{\mathcal{G}}(\bm{X}),$
where $\mu_{\max}$ denotes the eigenvalue of parameter matrix $\bm{W}_{l}$
that has the greatest magnitude.

  * •

$E_{\mathcal{G}}(\sigma(\bm{X}))\leq E_{\mathcal{G}}(\bm{X})$.

######  Theorem 4 (Dirichlet Energy)

Let us consider the matrix formulation of DSMP as (16). Suppose
$\bm{\theta}_{r,l}$ is symmetric positive semi-definite and let

$\mathbf{C}:=\max\\{\mu_{0}(\bm{\theta}_{0}),\max_{\begin{subarray}{c}r=1,\dots,K\\\
l=1,\dots,J\end{subarray}}\mu_{r,l}(\bm{\theta}_{r,l})\\},$

where $\mu_{0}(\bm{\theta}_{0})$ and $\mu_{r,l}(\bm{\theta}_{r,l})$ represent
the largest eigenvalue of matrix $\bm{\theta}_{0},\bm{\theta}_{r,l}$
respectively. For parameter matrix $\mathbf{P}$ of MLPs $\Gamma$, suppose its
largest eigenvalue is consistently bounded by $\mu_{\max}$, then we can bound
the Dirichlet energy of the graph node feature $\bm{X}^{(t)}$ at layer $t$ by

$E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)\leq
E_{\mathcal{G}}\left(\bm{X}^{(t)}\right)\leq\Bigl{(}1+\frac{\sqrt{5}}{2}\mathbf{C}\mu_{\max}\Bigr{)}^{2}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right).$

Proof  First, we consider the energy of term $\hat{\bm{Z}}^{(t)}$, where
$\bm{Z}^{(t)}=\Gamma(\hat{\bm{Z}}^{(t)})=\sigma(\mathbf{P}\hat{\bm{Z}}^{(t)})$.
By definition of $\hat{\bm{Z}}^{(t)}$, we have

$E_{\mathcal{G}}\left(\hat{\bm{Z}}^{(t)}\right)=E_{\mathcal{G}}\left(\bm{\theta}_{0}\bm{{\mathcal{W}}}_{0,J}\bm{X}^{(t-1)}\right)+E_{\mathcal{G}}\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{\theta}_{r,l}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}\bm{X}^{(t-1)}\right).$

By Lemma 2 and the fact that $\mathcal{L}={\bm{U}}\Lambda{\bm{U}}^{\top}$, we
have

$\displaystyle
E_{\mathcal{G}}\left(\bm{\theta}_{0}\bm{{\mathcal{W}}}_{0,J}\bm{X}^{(t-1)}\right)\leq\mathbf{C}^{2}\mathrm{tr}\left[\left(\bm{X}^{(t-1)}\right)^{\top}{\bm{U}}\widehat{\alpha}\left(\frac{\Lambda}{2}\right)^{2}\Lambda{\bm{U}}^{\top}\bm{X}^{(t-1)}\right],$

and according to (3), we have

$\displaystyle
E_{\mathcal{G}}\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{\theta}_{r,l}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}\bm{X}^{(t-1)}\right)$
$\displaystyle\leq\mathbf{C}^{2}\mathrm{tr}\left[\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{X}^{(t-1)\top}\bm{{\mathcal{W}}}_{r,l}^{\top}\bm{{\mathcal{W}}}_{0,J}^{\top}\right)\left(\bm{{\mathcal{W}}}_{0,J}\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{{\mathcal{W}}}_{r,l}\bm{X}^{(t-1)}\right)\right]$
$\displaystyle=\mathbf{C}^{2}\mathrm{tr}\left[\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{X}^{(t-1)\top}\bm{{\mathcal{W}}}_{r,l}^{\top}\right){\bm{U}}\widehat{\alpha}\left(\frac{\Lambda}{2}\right)^{2}\Lambda{\bm{U}}^{\top}\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{{\mathcal{W}}}_{r,l}\bm{X}^{(t-1)}\right)\right]$
$\displaystyle=\mathbf{C}^{2}\mathrm{tr}\left[\bm{X}^{(t-1)\top}{\bm{U}}\widehat{\alpha}\left(\frac{\Lambda}{2}\right)^{2}\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\widehat{\beta^{(r)}}\left(\frac{\Lambda}{2^{l}}\right)\right)^{2}\Lambda{\bm{U}}^{\top}\bm{X}^{(t-1)}\right].$

Since $\Lambda=\operatorname{diag}(\lambda_{1},\dots,\lambda_{N})$, we know
that

$\left|\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)\right|^{2}+\sum_{r=1}^{K}\sum_{l=1}^{J}\left|\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\right|^{2}=1,\forall\ell=1,\ldots,N.$

Therefore, $\forall\ell=1,\ldots,N,$ we have
$\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)^{2}\leq 1,$ and

$\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)^{2}\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\right)^{2}\leq\left(\frac{\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)^{2}+\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\right)^{2}}{2}\right)^{2}\leq\frac{1}{4}.$
(24)

Then, we have
$E_{\mathcal{G}}\left(\hat{\bm{Z}}^{(t)}\right)\leq\mathbf{C}^{2}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)+\frac{\mathbf{C}^{2}}{4}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)\leq\frac{5}{4}\mathbf{C}^{2}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)$.

According to Lemma 3,

$\displaystyle E_{\mathcal{G}}\left(\bm{Z}^{(t)}\right)$
$\displaystyle=E_{\mathcal{G}}\left(\sigma(\mathbf{P}\hat{\bm{Z}}^{(t)})\right)\leq\mu_{max}^{2}E_{\mathcal{G}}\left(\hat{\bm{Z}}^{(t)}\right)$
$\displaystyle\leq\mu_{max}^{2}\left(\mathbf{C}^{2}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)+\frac{\mathbf{C}^{2}}{4}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)\right)$
$\displaystyle\leq\frac{5}{4}\mathbf{C}^{2}\mu_{max}^{2}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right).$

By the definition of Dirichlet energy and the fact that the trace is a linear
operator, we have

$\displaystyle E_{\mathcal{G}}\left(\bm{X}^{(t)}\right)$
$\displaystyle=\mathrm{tr}\left(\left(\bm{X}^{(t)}\right)^{\top}\mathcal{L}\bm{X}^{(t)}\right)$
$\displaystyle=\mathrm{tr}\left[\left(\bm{X}^{(t-1)}+\bm{Z}^{(t)}\right)^{\top}\mathcal{L}\left(\bm{X}^{(t-1)}+\bm{Z}^{(t)}\right)\right]$
$\displaystyle=E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)+E_{\mathcal{G}}\left(\bm{Z}^{(t)}\right)+2\mathrm{tr}\left(\left(\bm{X}^{(t-1)}\right)^{\top}\mathcal{L}\bm{Z}^{(t)}\right).$

Again, by Lemma 2 and the cyclic property of the trace operator, we have

$\displaystyle\mathrm{tr}\left(\left(\bm{X}^{(t-1)}\right)^{\top}\mathcal{L}\bm{Z}^{(t)}\right)$
$\displaystyle=\mathrm{tr}\left(\mathcal{L}\bm{Z}^{(t)}\left(\bm{X}^{(t-1)}\right)^{\top}\right)$
$\displaystyle=\mathrm{tr}\left(\mathcal{L}\right)\mathrm{tr}\left(\bm{Z}^{(t)}\left(\bm{X}^{(t-1)}\right)^{\top}\right)$
$\displaystyle\leq\mathrm{tr}\left(\mathcal{L}\right)\sqrt{\mathrm{tr}(\left(\bm{Z}^{(t)}\right)^{\top}\bm{Z}^{(t)})}\sqrt{\mathrm{tr}(\left(\bm{X}^{(t-1)}\right)^{\top}\bm{X}^{(t-1)})}$
$\displaystyle\leq\sqrt{\mathrm{tr}\left(\mathcal{L}\right)\mathrm{tr}(\left({\bm{Z}^{(t)}}\right)^{\top}\bm{Z}^{(t)})}\sqrt{\mathrm{tr}\left(\mathcal{L}\right)\mathrm{tr}(\left(\bm{X}^{(t-1)}\right)^{\top}\bm{X}^{(t-1)})}$
$\displaystyle\leq\sqrt{E_{\mathcal{G}}\left(\bm{Z}^{(t)}\right)}\sqrt{E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)}$
$\displaystyle\leq\frac{\sqrt{5}}{2}\mathbf{C}\mu_{max}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right).$

Combining the above inequalities, we have

$E_{\mathcal{G}}\left(\bm{X}^{(t)}\right)\leq(1+\frac{\sqrt{5}}{2}\mathbf{C}\mu_{max})^{2}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right).$

Combining the fact that $E_{\mathcal{G}}(Z^{(t)})\geq 0$, we have

$E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right)\leq
E_{\mathcal{G}}\left(\bm{X}^{(t)}\right)\leq(1+\frac{\sqrt{5}}{2}\mathbf{C}\mu_{max})^{2}E_{\mathcal{G}}\left(\bm{X}^{(t-1)}\right).$

The lower and upper bounds suggest that the Dirichlet energy of the DSMP model
will neither diminish to zero nor escalate to an excessively large value.

###  4.4 Limited Over-squashing in DSMP

In this part, we will prove that the DSMP framework can alleviate the over-
squashing problem by connecting multi-hop neighbors. Let $\mathcal{R}$
represent a rewire operator, changing the topology of a graph. The rewired
version of graph $\mathcal{G}$ is denoted as $\mathcal{R}(\mathcal{G})$. We
will use (6) to help construct the transformation matrices and complete the
proof.

Figure 3: The entries of the multi-scale framelet transform may differ from
the adjacency matrix of the input graph. The transformation matrices are
considered as the edge weights and are visualized using a color bar.

###### Lemma 5

Given an adjacency matrix $\mathcal{A}$ of an undirected attributed graph
$\mathcal{G}$, the resulting matrix multiplication of $\mathcal{A}$ is denser
than $\mathcal{A}$.

Proof  When we perform matrix multiplication $\mathcal{A}\times\mathcal{A}$,
we define any entry $\mathcal{B}_{ij}$ of the resulting matrix $\mathcal{B}$
as:

$\mathcal{B}_{ij}=\sum_{k=1}^{n}\mathcal{A}_{ik}\cdot\mathcal{A}_{kj}.$

For $\mathcal{B}_{ij}$ to be non-zero, there must exist at least one $k$ such
that $\mathcal{A}_{ik}=1$ and $\mathcal{A}_{kj}=1$. This implies that there is
a path of length 2 between vertices $i$ and $j$ through vertex $k$.

In the original matrix $\mathcal{A}$, an edge between $i$ and $j$ (i.e.,
$\mathcal{A}_{ij}=1$) directly represents a path of length 1. When we move to
$\mathcal{B}$, we are considering paths of length 2. It is clear that any path
of length 1 in $\mathcal{A}$ will still represent a path (now of length 2) in
$\mathcal{B}$, and thus, the non-zero entry in $\mathcal{A}$ will still be
present in $\mathcal{B}$.

Moreover, $\mathcal{B}$ can have non-zero entries for pairs of vertices that
did not have a direct edge between them in $\mathcal{A}$ but were connected
via another vertex, thus representing a path of length 2. This cannot decrease
the number of non-zero entries; it can only increase or keep them the same.

######  Theorem 6 (Over-squashing)

Let us consider the matrix formulation of DSMP as (16), where the original
graph $\mathcal{G}$ is considered to be rewired as $\mathcal{R}(\mathcal{G})$
by framelet transformation. Designate the binary adjacency matrices of
$\mathcal{G}$ and $\mathcal{R}(\mathcal{G})$ as $\mathcal{A}$ and
$\mathcal{R}(\mathcal{A})$, respectively. We demonstrate that the rewired
graph exhibits denser connectivity, i.e.

$\mathcal{R}(\mathcal{A})-\mathcal{A}\geq 0,$

thus the DSMP model has a less over-squashing problem.

Proof  The issue of over-squashing primarily stems from the topology of the
graph structure. This implies that the over-squashing phenomenon can be
mitigated by increasing the number of edges. Consequently, we aim to
demonstrate that the original adjacency matrix $\mathcal{A}$, corresponding to
the topology of the original graph $\mathcal{G}$, contains less or equal
number of non-zero elements compared to the graph shift operator
$\bm{{\mathcal{W}}}_{0,J}+\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}$,
corresponding to the topology of a rewired graph $\mathcal{R}(\mathcal{G})$.

In (16), we utilize a GAT network to extract embedding from each rewired
graph. This network has the capability to learn the significance of node
relationships, with topology playing a crucial role in the process. Hence, we
consider the binary version of $\bm{{\mathcal{W}}}^{\natural}_{r,l}$, denoted
as $\mathcal{R}(\bm{{\mathcal{W}}}^{\natural}_{r,l})$, and we need to prove
$\mathcal{R}(\bm{{\mathcal{W}}}^{\natural}_{r,l})-\mathcal{A}\geq 0$.

According to (6), the transformation matrices of the multi-scale framelet have
the same dimensions as the adjacency matrix. As depicted in Figure 3, these
transformation matrices can be viewed as new graph structures with different
edge weights. Since the rewired graphs possess distinct adjacency matrices, we
analyze their over-squashing problem separately. Consequently, we obtain:

$\displaystyle\mathcal{R}(\bm{{\mathcal{W}}}^{\natural}_{r,l})$
$\displaystyle=\overline{\bm{{\mathcal{W}}}^{\natural}_{r,l}}=\overline{\mathcal{T}_{r}\left(2^{-R+l-1}\mathcal{L}\right)\mathcal{T}_{0}\left(2^{-R+l-2}\mathcal{L}\right)\ldots\mathcal{T}_{0}\left(2^{-R}\right)}=\overline{\underbrace{A\cdots
A}_{r}},l=2,\dots,J$ $\displaystyle\overline{\bm{{\mathcal{W}}}}_{i,j}$
$\displaystyle=\begin{cases}1,&\bm{{\mathcal{W}}}_{i,j}\neq 0\\\ 0,&{\rm
otherwise}\end{cases}.$

In the above equations, $\mathcal{A}$ is a binary adjacency matrix. According
to Lemma. 5, the matrix multiplication of $\mathcal{A}$ will not decrease the
non-zero entries. Therefore, the original graph is no more dense than the
rewired graph. When contrasted with message passing models such as GCN, GAT,
GIN, and others that aggregate information from one-hop neighbors on the
original graph, the DSMP model should have less over-squashing problem.

###  4.5 Stability of DSMP

In this section, we examine how DSMP improves the stability of message passing
against perturbation in the input node features.

###### Lemma 7

For framelet transforms on graph $\mathcal{G}$, and graph signal $\bm{X}$ in
$\ell_{2}(\mathcal{G})$,

$\displaystyle\left\|\sum_{r=1}^{K}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}\bm{X}\right\|\leq\left\|\bm{X}\right\|,$
$\displaystyle\left\|\bm{{\mathcal{W}}}_{0,J}\bm{X}\right\|^{2}\leq\|\bm{X}\|^{2}.$

Proof  By exploiting the orthonormality of $\bm{u}_{\ell}$, we can express the
coefficients as follows:

$\left\langle\bm{X},\bm{\varphi}_{l,p}\right\rangle=\sum_{\ell=1}^{N}\overline{\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2^{l}}\right)}\widehat{\bm{X}}_{\ell}\>\bm{u}_{\ell}(p),\quad\left\langle\bm{X},\bm{\psi}_{l,p}^{r}\right\rangle=\sum_{\ell=1}^{N}\overline{\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)}\widehat{\bm{X}}_{\ell}\>\bm{u}_{\ell}(p).$

By referencing (22), we can derive the following expression:

$\displaystyle\sum_{\ell=1}^{N}\sum_{r=1}^{K}\left|\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\right|^{2}\left|\widehat{\bm{X}}_{\ell}\right|^{2}\>\left\|\bm{u}_{\ell}\right\|^{2}$
$\displaystyle=\sum_{\ell=1}^{N}\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)^{2}\left[\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2^{l}}\right)^{2}-\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2^{l}}\right)^{2}\right]\bigl{|}\widehat{\bm{X}}_{\ell}\bigr{|}^{2}\>\|\bm{u}_{\ell}\|^{2}$
$\displaystyle=\sum_{\ell=1}^{N}\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)^{2}\sum_{r=1}^{K}\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)^{2}\bigl{|}\widehat{\bm{X}}_{\ell}\bigr{|}^{2}$
$\displaystyle\leq\|\bm{X}\|^{2},$

where the inequality of the last line follows according to (23).

Next, by utilizing (3), we can establish the following relationship:

$\displaystyle\left\|\sum_{r=1}^{K}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}\bm{X}\right\|^{2}=$
$\displaystyle\left\|\left\langle\bm{X},\bm{\varphi}_{0,\cdot}\bm{\psi}_{l,\cdot}^{r}\right\rangle\right\|^{2}$
$\displaystyle=$
$\displaystyle\sum_{\ell=1}^{N}\sum_{r=1}^{K}\left|\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)\widehat{\beta^{(r)}}\left(\frac{\lambda_{\ell}}{2^{l}}\right)\right|^{2}\left|\widehat{\bm{X}}_{\ell}\right|^{2}\>\left\|\bm{u}_{\ell}\right\|^{2}$
$\displaystyle\leq$ $\displaystyle\|\bm{X}\|^{2}.$

Similarly, we have

$\displaystyle\left\|\bm{{\mathcal{W}}}_{0,J}\bm{X}\right\|^{2}=$
$\displaystyle\left\|\left\langle\bm{X},\bm{\varphi}_{0,\cdot}\right\rangle\right\|^{2}$
$\displaystyle=$
$\displaystyle\sum_{\ell=1}^{N}\widehat{\alpha}\left(\frac{\lambda_{\ell}}{2}\right)^{2}|\widehat{\bm{X}}_{\ell}\bigr{|}^{2}\>\|\bm{u}_{\ell}\|^{2}$
$\displaystyle\leq$ $\displaystyle\|\bm{X}\|^{2}.$

thus completing the proof.

######  Theorem 8 (Stability)

Suppose the DSMP has bounded parameters in all layers:
$\|\bm{\theta}_{0}\|\leq\mathbf{C}_{0}$,
$\|\bm{\theta}_{r,l}\|\leq\mathbf{C}_{r,l}$ with constants $\mathbf{C}_{r,l}$
for $r=1,\dots,K,l=1,\dots,J$. For parameter matrix $\mathbf{P}$ of MLPs
$\Gamma$, suppose $\|\mathbf{P}\|\leq\mathbf{C}_{p}$. Then, the DSMP is
Lipschitz continuous in the sense that

$\left\|{\bm{X}}^{(t)}-\widetilde{{\bm{X}}}^{(t)}\right\|\leq
C^{t}\left\|{\bm{X}}^{(0)}-\widetilde{{\bm{X}}}^{(0)}\right\|,$

where

$C:=\mathbf{C}_{p}\left(\max_{\begin{subarray}{c}r=1,\dots,K\\\
l=1,\dots,J\end{subarray}}\mathbf{C}_{r,l}+\mathbf{C}_{0}\right),$

${\bm{X}}^{(0)}$ and $\widetilde{{\bm{X}}}^{(0)}$ are the initial graph node
features, and ${\bm{X}}^{(t)}$ and $\widetilde{{\bm{X}}}^{(t)}$ are node
features propagated at layer $t$.

Proof  Let

${\bm{X}}^{(t)}={\bm{X}}^{(t-1)}+{\bm{Z}}^{(t)},\quad\widetilde{{\bm{X}}}^{(t)}=\widetilde{{\bm{X}}}^{(t-1)}+\widetilde{{\bm{Z}}}^{(t)},$

with initialization
${\bm{X}}^{(0)},\widetilde{{\bm{X}}}^{(0)}\in\ell_{2}({\mathcal{G}})$. It thus
holds that

$\displaystyle\left\|{\bm{X}}^{(t)}-\widetilde{{\bm{X}}}^{(t)}\right\|$
$\displaystyle\leq\left\|{\bm{X}}^{(t-1)}-\widetilde{{\bm{X}}}^{(t-1)}\right\|+\left\|{\bm{Z}}^{(t)}-\widetilde{{\bm{Z}}}^{(t)}\right\|,$

where

${\bm{Z}}^{(t)}=\Gamma\left(\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}{\bm{X}}^{(t-1)}\bm{\theta}_{r,l}+\bm{{\mathcal{W}}}_{0,J}{\bm{X}}^{(t-1)}\bm{\theta}_{0}\right).$

Let
${\bm{Y}}^{(t)}:=\sum_{r=1}^{K}\sum_{l=1}^{J}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}{\bm{X}}^{(t-1)}\bm{\theta}_{r,l}+\bm{{\mathcal{W}}}_{0,J}{\bm{X}}^{(t-1)}\bm{\theta}_{0},$
then we have

$\displaystyle\left\|{\bm{Y}}^{(t)}-\widetilde{{\bm{Y}}}^{(t)}\right\|\leq$
$\displaystyle\sum_{l=1}^{J}\left\|\sum_{r=1}^{K}\bm{{\mathcal{W}}}_{0,J}\bm{{\mathcal{W}}}_{r,l}\left({\bm{X}}^{(t-1)}-\widetilde{{\bm{X}}}^{(t-1)}\right)\bm{\theta}_{r,l}\right\|+\left\|\bm{{\mathcal{W}}}_{0,J}\left({\bm{X}}^{(t-1)}-\widetilde{{\bm{X}}}^{(t-1)}\right)\bm{\theta}_{0}\right\|$
$\displaystyle\leq$
$\displaystyle\sum_{l=1}^{J}\mathbf{C}_{r,l}\left\|{\bm{X}}^{(t-1)}-\widetilde{{\bm{X}}}^{(t-1)}\right\|+\mathbf{C}_{0}\>\left\|{\bm{X}}^{(t-1)}-\widetilde{{\bm{X}}}^{(t-1)}\right\|$
$\displaystyle\leq$
$\displaystyle\Bigl{(}\max_{\begin{subarray}{c}r=1,\dots,K\\\
l=1,\dots,J\end{subarray}}\mathbf{C}_{r,l}+\mathbf{C}_{0}\Bigr{)}\left\|{\bm{X}}^{(t-1)}-\widetilde{{\bm{X}}}^{(t-1)}\right\|.$

Therefore,

$\displaystyle\left\|{\bm{Z}}^{(t)}-\widetilde{{\bm{Z}}}^{(t)}\right\|$
$\displaystyle=\left\|\Gamma\left({\bm{Y}}^{(t)}-\widetilde{{\bm{Y}}}^{(t)}\right)\right\|\leq\mathbf{C}_{p}\left\|{\bm{Y}}^{(t)}-\widetilde{{\bm{Y}}}^{(t)}\right\|$
$\displaystyle\leq\mathbf{C}_{p}\Bigl{(}\max_{\begin{subarray}{c}r=1,\dots,K\\\
l=1,\dots,J\end{subarray}}\mathbf{C}_{r,l}+\mathbf{C}_{0}\Bigr{)}\left\|{\bm{X}}^{(t-1)}-\widetilde{{\bm{X}}}^{(t-1)}\right\|,$

which leads to

$\left\|{\bm{X}}^{(t)}-\widetilde{{\bm{X}}}^{(t)}\right\|\leq
C\left\|{\bm{X}}^{(t-1)}-\widetilde{{\bm{X}}}^{(t-1)}\right\|\leq
C^{t}\left\|{\bm{X}}^{(0)}-\widetilde{{\bm{X}}}^{(0)}\right\|,$

thus completing the proof.

##  5 Numerical Analysis

In this section, we empirically verify the effectiveness and properties of the
DSMP model on graph-level prediction tasks. We assess the performance of the
proposed model by conducting experiments on several popular benchmark datasets
with varying sample sizes. These datasets are commonly used in the field and
serve as reliable indicators of the model’s capabilities in predicting graph-
level properties. We implemented the proposed model with PyTorch-Geometric
(Fey and Lenssen, 2019a) (version 2.4.0) and PyTorch (version 1.12.1) and
tested on NVIDIA Tesla A100 GPU with 6,912 CUDA cores and 80GB HBM2 mounted on
an HPC cluster. The codes for the experiments are available at
https://github.com/jyh6681/DSMP.

###  5.1 Datasets and Experimental Protocol

We conducted numerical experiments on a diverse set of widely used node
classification and graph classification tasks to evaluate the effectiveness of
our proposed model. The selected graph classification benchmarks include seven
datasets: D&D, PROTEINS, NCI1, Mutagenicity, COLLAB, OGBG-MOLHIV, and QM7.
These datasets are commonly used in the field and have been extensively
studied for tasks such as protein categorization, chemical compound
identification, and molecular property prediction. Notably, the OGBG-MOLHIV
dataset was used for large-scale molecule classification.

To validate our proposed model, we compared it with state-of-the-art (SOTA)
baseline MPNN models, specifically GCN and GIN, which are popular GNN models.
In our experimental design, we prioritized fairness and comprehensiveness
rather than solely focusing on achieving the best possible performance on each
dataset. To ensure consistency and eliminate hyperparameter tuning as a
potential source of performance gain, we used the same GNN and optimization
hyperparameters across the same search space for each task and baseline model.
Each dataset follows the standard public split and processing rules in
PyTorch-Geometric (Fey and Lenssen, 2019b) or OGB (Hu et al., 2020). The
training process was stopped when the validation loss ceased to improve for 20
consecutive epochs or reached a maximum number of epochs.

To obtain reliable performance estimates, we performed 10 random trials for
each experiment and reported the mean test accuracy along with the variance
value. For the graph-level classification tasks using TUDataset benchmarks, we
reported the mean test accuracy, while for the OGBG-MOLHIV dataset, we
reported the ROC-AUC scores. The regression task on the QM7 dataset was
evaluated using the mean square error (MSE) metric. Further details of the
experiments are reported in Table 5. These experiments provide a comprehensive
evaluation of our proposed model’s performance compared to the baseline
models, demonstrating its effectiveness across various node and graph
classification tasks.

#### Setup

We train DSMP with $2$ levels and $2$ scales, resulting in $3$ scattering
coefficient matrices feeding into one message passing layer. The output graph
embedding is processed by a softmax activation function for label prediction.
Regarding hyperparameters, a grid search is applied to the learning rate,
weight decay, hidden units, and dropout ratio to obtain the optimal
combination from the search space as shown in Table 1.

The model is trained with the Adam optimizer (Kingma and Ba, 2015). The
maximum number of training epochs is 500 for ogbg-molhiv and 200 for other
benchmarks. All benchmarks follow standard public data splitting and
preprocessing procedures. The average test accuracy and its standard deviation
are reported from 10 repetitive runs.

Table 1: Grid search space for the hyperparameters.

Hyperparameters | Search Space
---|---
Learning Rate | 1$\mathrm{e}$-4, 5$\mathrm{e}$-4, 1$\mathrm{e}$-3, 5$\mathrm{e}$-3, 1$\mathrm{e}$-2
Hidden Size |  $16$, $32$, $64$, $128$
Weight Decay ($L_{2}$) | 1$\mathrm{e}$-4, 5$\mathrm{e}$-4, 1$\mathrm{e}$-3, 5$\mathrm{e}$-3
Batch Size |  $8$, $64$, $128$, $1024$

###  5.2 Evolution of Dirichlet Energy

In order to assess the smoothness of graph embeddings, we employ the Dirichlet
energy metric. Figure LABEL:fig:dirichlet illustrates that the proposed DSMP
can preserve the Dirichlet energy of graph signal even if the layer goes deep,
where we sample a graph from the PROTEINS dataset and feed it into message
passing models with varying numbers of layers.

###  5.3 Experiments and Analysis on Over-squashing Issue

In recent research (Karhadkar et al., 2023), the severity of the over-
squashing issue within a graph structure is quantitatively assessed using the
spectral gap, as outlined in section 2.4. The spectral gap is intimately tied
to the graph’s connectedness, as the magnitude of its eigenvalues provides
crucial insights. Specifically, the spectral gap can offer a proxy for
structural bottlenecks through the Cheeger constant (Chung, 1997; Cheeger,
2015; Alon and Milman, 1984). In essence, a larger spectral gap implies a
graph with fewer structural bottlenecks, thereby indicating a reduced
likelihood of encountering the over-squashing problem in MPNNs.

To empirically validate this, we computed the spectral gap for each graph
instance in several publicly available datasets commonly used for graph
classification tasks. Figure LABEL:fig:gap_range presents a histogram of the
spectral gaps, along with their mean values, providing a visual representation
of the distribution and typical magnitudes observed across the datasets.

Numerically, we compare the proposed DSMP model to no graph rewiring, as well
as three SOTA baselines: SDRF (Topping et al., 2022), FoSR (Karhadkar et al.,
2023), and BORF (Nguyen et al., 2023). We also evaluate the graph transformer
model U2GNN (Nguyen et al., 2022) to verify the impact of a fully connected
graph structure.

We statistically measure the magnitude of the over-squashing issue on some
public graph datasets and verify the effectiveness of DSMP, where the baseline
rewiring methods have been previously tested (pmlr-v202-nguyen23c; Nguyen et
al., 2022). Compared to no rewiring GNN model, the baselines generally yield
relatively significant gains except on IMDB as Table 2 shows. By analyzing the
“bottlenecked-ness” from Figure LABEL:fig:gap_range and Table 2 results, we
conclude these baseline rewiring methods are more impactful when inherent
graph bottlenecks are severe. However, DSMP achieves the top performance
across three graph classification tasks despite the over-squashing issue among
these datasets. This indicates that DSMP more strongly addresses the over-
squashing problem.

Table 2: The comparison with the graph rewiring methods on solving the over-
squashing problem.

| ENZYMES | IMDB | MUTAG | PROTEINS
---|---|---|---|---
# graphs | $600$ | $1,000$ | $188$ | $1,113$
# classes | $6$ | $2$ | $2$ | $2$
# nodes | $2-126$ | $12-136$ | $10-28$ | $4-620$
Avg. # nodes | $32.63$ | $19.77$ | $17.93$ | $39.06$
Avg. # edges | $124.27$ | $193.062$ | $39.58$ | $145.63$
U2GNN |  $37.5$$\pm 1.8$ |  $53.6$$\pm 3.5$ |  $\bf{88.5}$$\pm 3.2$ |  $78.5$$\pm 4.1$
GCN |  $25.5$$\pm 1.3$ |  $49.3$$\pm 1.0$ |  $68.8$$\pm 2.1$ |  $70.6$$\pm 1.0$
GCN+SDRF |  $26.1$$\pm 1.1$ |  $49.1$$\pm 0.9$ |  $70.5$$\pm 2.1$ |  $71.4$$\pm 0.8$
GCN+FoSR |  $27.4$$\pm 1.1$ |  $49.6$$\pm 0.8$ |  $75.6$$\pm 1.7$ |  $72.3$$\pm 0.9$
GCN+BORF |  $24.7$$\pm 1.0$ |  $50.1$$\pm 0.9$ |  $75.8$$\pm 1.9$ |  $71.0$$\pm 0.8$
GIN |  $31.3$$\pm 1.2$ |  $69.0$$\pm 1.3$ |  $75.5$$\pm 2.9$ |  $69.7$$\pm 1.0$
GIN+SDRF |  $33.5$$\pm 1.3$ |  $68.6$$\pm 1.2$ |  $77.3$$\pm 2.3$ |  $72.2$$\pm 0.9$
GIN+FoSR |  $25.3$$\pm 1.2$ |  $69.5$$\pm 1.1$ |  $75.2$$\pm 3.0$ |  $74.2$$\pm 0.8$
GIN+BORF |  $35.5$$\pm 1.2$ |  $71.3$$\pm 1.5$ |  $80.8$$\pm 2.5$ |  $71.3$$\pm 1.0$
DSMP |  $\bf{40.6}$$\pm 0.5$ |  $\bf{77.3}$$\pm 0.6$ |  $85.0$$\pm 0.25$ |  $\bf{79.6}$$\pm 0.7$

1 0.38,0.75- $100\%$ 0.33,0.25- $50\%$ 0.28,-0.35- $0\%$ 0.4,-0.95- $-50\%$
0.45,-1.45- $-100\%$

###  5.4 Stability of DSMP

To assess the stability of our proposed DSMP model in the semi-supervised node
classification task, we conducted experiments on eight widely used benchmark
datasets. The characteristics of these datasets are summarized in Table 3. To
introduce perturbations to the node attributes, we employed the _attribute
injection_ method. This approach involves swapping attributes between distinct
nodes in a random subgraph, similar to the technique described in (Ding et
al., 2019).

A certain number of targeted nodes are randomly selected from the input graph.
For each selected node $v_{i}$, we randomly pick $k$ nodes from the graph and
select the node $v_{j}$ whose attributes deviate the most from node $v_{i}$
among the $k$ nodes by maximizing the Euclidean distance of node attributes,
i.e., $\left\|x_{i}-x_{j}\right\|_{2}$. Then, we substitute the attributes
$x_{i}$ of node $v_{i}$ with $x_{j}$. Specifically, we set the size of
candidates $k=100$ for small datasets, i.e., Cora and Citeseer, and $k=500$
for relatively larger datasets, i.e., PubMed, Coauthor-CS, and Wiki-CS.

By introducing these controlled perturbations, we aimed to evaluate the DSMP
model’s resilience to attribute noise and assess its stability in node
classification tasks. The results of these experiments will provide insights
into the model’s robustness and its ability to handle real-world scenarios
with noisy or perturbed data.

Table 3: Summary of the datasets for node classification tasks.

| Cora | Citeseer | PubMed | Wiki-CS | Coauthor-CS | Wisconsin | Texas | OGB-arxiv
---|---|---|---|---|---|---|---|---
# Nodes | $2,708$ | $3,327$ | $19,717$ | $11,701$ | $18,333$ | $251$ | $183$ | $169,343$
# Edges | $5,429$ | $4,732$ | $44,338$ | $216,123$ | $100,227$ | $499$ | $309$ | $1,166,243$
# Features | $1,433$ | $3,703$ | $500$ | $300$ | $6,805$ | $1,7033$ | $1,703$ | $128$
# Classes | $7$ | $6$ | $3$ | $10$ | $15$ | $5$ | $5$ | $40$
# Training Nodes | $140$ | $120$ | $60$ | $580$ | $300$ | $120$ | $87$ | $90,941$
# Validation Nodes | $500$ | $500$ | $500$ | $1769$ | $200$ | $80$ | $59$ | $29,799$
# Test Nodes | $1,000$ | $1,000$ | $1,000$ | $5847$ | $1000$ | $51$ | $37$ | $48,603$
Label Rate | $0.052$ | $0.036$ | $0.003$ | $0.050$ | $0.016$ | $0.478$ | $0.475$ | $0.537$
Feature Scale | $\\{0,1\\}$ | $\\{0,1\\}$ | $[0,1.263]$ | $[-3,3]$ | $\\{0,1\\}$ | $\\{0,1\\}$ | $\\{0,1\\}$ | $[-1.389,1.639]$

To verify the stability of the proposed model, we compare the DSMP to three
popular graph denoising models with different design philosophies: APPNP
(Klicpera et al., 2018) avoids global smoothness with residual connections;
GNNGuard (Zhang and Zitnik, 2020) modifies the neighbor relevance in message
passing to mitigate local corruption; ElasticGNN (Liu et al., 2021b) and
AirGNN (Liu et al., 2021a) pursue local smoothness with sparseness
regularizers. We also train a $2$-layer GCN (Kipf and Welling, 2017) as the
baseline method, which smooths graph signals by minimizing the Dirichlet
energy in the spatial domain.

Table 4: Average performance for node classification over $10$ repetitions.
Due to a large number of nodes, the APPNP model encountered an out-of-memory
issue on OGB-arxiv dataset.

Module | Cora | Citeseer | PubMed | Coauthor-CS | Wiki-CS | Wisconsin | Texas | OGB-arxiv
---|---|---|---|---|---|---|---|---
clean |  $81.26$$\pm 0.65$ |  $71.77$$\pm 0.29$ |  $79.01$$\pm 0.44$ |  $90.19$$\pm 0.48$ |  $77.62$$\pm 0.26$ |  $56.47$$\pm 5.26$ |  $65.14$$\pm 1.46$ |  $71.10$$\pm 0.21$
GCN |  $69.06$$\pm 0.74$ |  $57.58$$\pm 0.71$ |  $67.69$$\pm 0.40$ |  $82.41$$\pm 0.23$ |  $65.44$$\pm 0.23$ |  $48.24$$\pm 3.19$ |  $58.92$$\pm 2.02$ |  $68.42$$\pm 0.15$
APPNP |  $68.46$$\pm 0.81$ |  $60.04$$\pm 0.59$ |  $68.70$$\pm 0.47$ |  $71.14$$\pm 0.54$ |  $56.53$$\pm 0.72$ |  $61.76$$\pm 5.21$ |  $59.46$$\pm 0.43$ | OOM
GNNGuard |  $61.96$$\pm 0.30$ |  $54.94$$\pm 1.00$ |  $68.50$$\pm 0.38$ |  $80.67$$\pm 0.88$ |  $65.69$$\pm 0.32$ |  $46.86$$\pm 1.06$ |  $59.19$$\pm 0.81$ |  $65.75$$\pm 0.32$
ElasticGNN |  $77.74$$\pm 0.79$ |  $64.61$$\pm 0.85$ |  $71.23$$\pm 0.21$ |  $79.91$$\pm 1.39$ |  $64.18$$\pm 0.53$ |  $53.33$$\pm 2.45$ |  $59.77$$\pm 3.24$ |  $41.34$$\pm 0.38$
AirGNN | $76.22{\scriptsize\pm 3.75}$ | $62.14{\scriptsize\pm 0.82}$ | $74.73{\scriptsize\pm 0.43}$ | $80.18{\scriptsize\pm 0.31}$ | $71.36{\scriptsize\pm 0.20}$ | $61.56{\scriptsize\pm 0.72}$ | $59.46{\scriptsize\pm 1.24}$ | $52.32{\scriptsize\pm 0.58}$
DSMP | $72.23{\scriptsize\pm 0.34}$ | $62.79{\scriptsize\pm 0.55}$ | $72.53{\scriptsize\pm 0.82}$ | $\bm{82.70}{\scriptsize\pm 0.37}$ | $\bm{75.54{\scriptsize\pm 0.27}}$ | $52.94{\scriptsize\pm 1.9}$ | $\bm{60.36{\scriptsize\pm 2.67}}$ | $66.73{\scriptsize\pm 0.56}$

###  5.5 Graph Classification Experiments

In this section, we compare the DSMP model with baseline models of graph
representation learning. Following the convention, we report the percentage
value of mean test accuracy for the classification tasks with TUDatasets and
ROC-AUC score for the ogbg-molhiv dataset. The mean performance scores are
averaged over 10 repetitions, with their standard deviations provided after
the $\pm$ signs. As indicated in Table 5, the DSMP model outperforms the other
baseline models in general.

Compared with GCN, the DSMP model outperforms by 12.8% at most on the PROTEINS
dataset, and by 0.9% at least on the QM7 dataset. When compared with GIN, our
model outperforms by 14.2% at most on the PROTEINS dataset, and by 0.24% at
least on the COLLAB dataset. The DSMP model consistently demonstrates superior
performance than the HDS-GAT model, further validating its effectiveness and
reliability. Our method shows a noticeable improvement compared to GCN on most
datasets, except for QM7 and COLLAB. As illustrated in Figure
LABEL:fig:freq_range, the observed performance difference can be attributed to
the fact that these two datasets contain limited high-frequency information.
Consequently, methods designed to extract high-frequency features may not be
as effective in this context.

Table 5: Summary of the datasets for graph-level classification tasks. QM7
dataset is not a typical dataset for graph-level classification and we use $R$
to denote it for a regression task.

| QM7 | PROTEINS | D&D | Mutagenicity | COLLAB | OGBG-MOLHIV | NCI1
---|---|---|---|---|---|---|---
# graphs | $7,165$ | $1,113$ | $1,178$ | $4,337$ | $5,000$ | $41,127$ | $4,110$
# classes | 1 ($R$) | $2$ | $2$ | $2$ | $3$ | $2$ | $2$
# nodes | $4-23$ | $4-620$ | $30-5,748$ | $30-417$ | $32-492$ | $2-222$ | $3-111$
Avg. # nodes | $15$ | $39$ | $284$ | $30$ | $74$ | $26$ | $30$
Avg. # edges | $123$ | $73$ | $716$ | $31$ | $2,458$ | $28$ | $32$
# Features | $0$ | $3$ | $89$ | $14$ | $0$ | $9$ | $37$
GCN |  $41.9$$\pm 1.1$ |  $70.6$$\pm 1.0$ |  $78.9$$\pm 3.3$ |  $80.7$$\pm 3.3$ |  $81.8$$\pm 1.3$ |  $74.9$$\pm 2.6$ |  $76.9$$\pm 1.7$
GIN |  $43.0$$\pm 1.5$ |  $69.7$$\pm 1.0$ |  $78.6$$\pm 1.3$ |  $81.3$$\pm 1.1$ |  $82.3$$\pm 0.9$ |  $77.4$$\pm 1.2$ |  $78.8$$\pm 1.5$
UFG |  $41.7$$\pm 0.8$ |  $77.8$$\pm 2.6$ |  $80.9$$\pm 1.7$ |  $81.6$$\pm 1.4$ |  $\bf{82.8}$$\pm 0.3$ |  $78.8$$\pm 0.5$ |  $77.9$$\pm 1.2$
HDS-GAT |  $42.2$$\pm 0.5$ |  $78.5$$\pm 1.5$ |  $79.8$$\pm 1.1$ |  $81.2$$\pm 0.5$ |  $82.2$$\pm 0.6$ |  $78.3$$\pm 0.8$ |  $78.2$$\pm 1.5$
DSMP |  $\bf{41.5}$$\pm 1.0$ |  $\bf{79.6}$$\pm 0.7$ |  $\bf{82.5}$$\pm 1.1$ |  $\bf{82.0}$$\pm 0.9$ |  $82.5$$\pm 0.3$ |  $\bf{79.2}$$\pm 1.7$ |  $\bf{79.8}$$\pm 1.0$

###  5.6 Computational Cost

#### Time complexity

Let $E=|\mathcal{E}|$ be the number of edges in the graph, $N$ be the number
of nodes and $d$ be the feature dimension of nodes. According to (11), the
computational cost for one propagation is $O\left(ENd\right)$. For the
propagation of the SMP module, as (15) shows, the computational cost is also
$O\left(ENd\right)$.

We also compare the computational cost of the DSMP model with some popular
MPNNs, which are considered as baseline models, in terms of testing time
complexity. In Table 6, we present the time cost for each model. We randomly
select a batch of data with $8$ graph samples from the PROTEINS dataset and
increase the number of layers in the message passing frameworks.

Table 6: Comparison of the time cost on PROTEINS dataset.

Layer Num | $1$ | $3$ | $50$
---|---|---|---
GCN | $0.0847$ | $0.0868$ | $0.1151$
GIN | $0.0857$ | $0.088$ | $0.1444$
DSMP | $0.0992$ | $0.1226$ | $0.3048$

#### Space complexity

In conventional message passing frameworks, the sparse edge connection matrix
and dense node feature matrix are stored in memory, resulting in a memory cost
of $2E+O(Nd)$. According to (16), the propagation of each SMP module requires
one low-pass framelet matrix and $K\times J$ high-pass framelet matrices. In
our numerical experiments, we set $K=1$ and $J=2$ to utilize scaling
information up to $2$, which connects only 2-hop neighbors as shown in Figure
3. Let $R$ be the number of unconnected 2-hop neighbors, then the memory cost
for SMP is $6E+2R+O(Nd)$.

##  6 Conclusion and Future Work

In conclusion, this work introduces the DSMP neural network as a solution to
address the stability, over-smoothing, and over-squashing issues commonly
observed in current GNN models. By integrating the multiscale and multilevel
scattering transform with a learnable message passing scheme, the DSMP model
offers a novel approach to processing graph-structured data. The proposed
message passing framework, viewed from the perspective of spectral methods,
has been theoretically proven to effectively mitigate the aforementioned
problems. Through extensive experiments, we have demonstrated that the DSMP
model outperforms existing GNNs in various graph reasoning tasks, highlighting
its superior performance and effectiveness.

In terms of future work, there are several potential directions to explore.
Firstly, further investigation can be conducted to enhance the scalability and
efficiency of the DSMP model, particularly in large-scale graph datasets.
Additionally, exploring the applicability of the DSMP model in other domains,
such as social network analysis or recommendation systems, could provide
valuable insights into its versatility and generalizability. Furthermore, the
interpretability of the DSMP model can be explored, enabling a better
understanding of its decision-making process. This could involve analyzing the
importance of different features or nodes in the graph and providing
explanations for the model’s predictions.

Overall, the proposed DSMP model presents a promising direction for improving
the performance and addressing the limitations of GNNs. The future work
mentioned above will contribute to further advancing the field of GNNs and
their applications.

## References

  * Alon and Milman (1984) Noga Alon and Vitali D Milman.  Eigenvalues, expanders and superconcentrators.  In _25th Annual Symposium on Foundations of Computer Science, 1984._ , pages 320–322, 1984.
  * Alon and Yahav (2021) Uri Alon and Eran Yahav.  On the bottleneck of graph neural networks and its practical implications.  In _ICLR_ , 2021.
  * Banerjee et al. (2022) Pradeep Kr Banerjee, Kedar Karhadkar, Yu Guang Wang, Uri Alon, and Guido Montúfar.  Oversquashing in gnns through the lens of information contraction and graph expansion.  In _2022 58th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , pages 1–8, 2022.
  * Bodnar et al. (2021) Cristian Bodnar, Fabrizio Frasca, Nina Otter, Yuguang Wang, Pietro Lio, Guido F Montufar, and Michael Bronstein.  Weisfeiler and Lehman go cellular: CW networks.  In _NeurIPS_ , volume 34, pages 2625–2640, 2021.
  * Brody et al. (2022) Shaked Brody, Uri Alon, and Eran Yahav.  How attentive are graph attention networks?  _ICLR_ , 2022.
  * Bruna and Mallat (2013) Joan Bruna and Stéphane Mallat.  Invariant scattering convolution networks.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 35(8):1872–1886, 2013.
  * Bruna et al. (2014) Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  In _ICLR_ , 2014.
  * Bumgardner et al. (2021) Bri Bumgardner, Farhan Tanvir, Khaled Mohammed Saifuddin, and Esra Akbas.  Drug-drug interaction prediction: a purely smiles based approach.  In _2021 IEEE International Conference on Big Data (Big Data)_ , pages 5571–5579, 2021.
  * Cai and Wang (2020) Chen Cai and Yusu Wang.  A note on over-smoothing for graph neural networks.  In _ICML 2020 Graph Representation Learning workshop_ , 2020.
  * Cheeger (1970) Jeff Cheeger.  A lower bound for the smallest eigenvalue of the Laplacian, problems in analysis (papers dedicated to salomon bochner, 1969), 1970.
  * Cheeger (2015) Jeff Cheeger.  A lower bound for the smallest eigenvalue of the Laplacian.  In _Problems in Analysis: A Symposium in Honor of Salomon Bochner (PMS-31)_ , pages 195–200, 2015.
  * Chen et al. (2020) Deli Chen, Yankai Lin, Wei Li, Peng Li, Jie Zhou, and Xu Sun.  Measuring and relieving the over-smoothing problem for graph neural networks from the topological view.  In _AAAI_ , volume 34, pages 3438–3445, 2020.
  * Chew et al. (2022) Joyce Chew, Holly Steach, Siddharth Viswanath, Hau-Tieng Wu, Matthew Hirn, Deanna Needell, Matthew D Vesely, Smita Krishnaswamy, and Michael Perlmutter.  The manifold scattering transform for high-dimensional point cloud data.  In _Topological, Algebraic and Geometric Learning Workshops 2022_ , pages 67–78, 2022.
  * Chung (1997) Fan RK Chung.  _Spectral graph theory_ , volume 92.  American Mathematical Soc., 1997.
  * Di Giovanni et al. (2022) Francesco Di Giovanni, James Rowbottom, Benjamin P Chamberlain, Thomas Markovich, and Michael M Bronstein.  Graph neural networks as gradient flows.  _arXiv preprint arXiv:2206.10991_ , 2022.
  * Di Giovanni et al. (2023) Francesco Di Giovanni, Lorenzo Giusti, Federico Barbero, Giulia Luise, Pietro Lio, and Michael M Bronstein.  On over-squashing in message passing neural networks: The impact of width, depth, and topology.  In _ICML_ , pages 7865–7885, 2023.
  * Ding et al. (2019) Kaize Ding, Jundong Li, Rohit Bhanushali, and Huan Liu.  Deep anomaly detection on attributed networks.  In _SIAM International Conference on Data Mining (SDM)_ , 2019.
  * Dong (2017) Bin Dong.  Sparse representation on graphs by tight wavelet frames and applications.  _Applied and Computational Harmonic Analysis_ , 42(3):452–479, 2017.
  * Fey and Lenssen (2019a) Matthias Fey and Jan E. Lenssen.  Fast graph representation learning with PyTorch Geometric.  In _ICLR Workshop on Representation Learning on Graphs and Manifolds_ , 2019a.
  * Fey and Lenssen (2019b) Matthias Fey and Jan Eric Lenssen.  Fast graph representation learning with PyTorch Geometric.  In _ICLR Workshop on Representation Learning on Graphs and Manifolds_ , 2019b.
  * Gama et al. (2019a) Fernando Gama, Alejandro Ribeiro, and Joan Bruna.  Diffusion scattering transforms on graphs.  In _ICLR_ , 2019a.
  * Gama et al. (2019b) Fernando Gama, Alejandro Ribeiro, and Joan Bruna.  Stability of graph scattering transforms.  In _NeurIPS_ , volume 32, 2019b.
  * Gama et al. (2020) Fernando Gama, Joan Bruna, and Alejandro Ribeiro.  Stability properties of graph neural networks.  _IEEE Transactions on Signal Processing_ , 68:5680–5695, 2020.
  * Gao et al. (2019) Feng Gao, Guy Wolf, and Matthew Hirn.  Geometric scattering for graph data analysis.  In _ICML_ , pages 2122–2131, 2019.
  * Gauthier et al. (2022) Shanel Gauthier, Benjamin Thérien, Laurent Alsene-Racicot, Muawiz Chaudhary, Irina Rish, Eugene Belilovsky, Michael Eickenberg, and Guy Wolf.  Parametric scattering networks.  In _CVPR_ , pages 5749–5758, 2022.
  * Gilmer et al. (2017) Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In _ICML_ , 2017.
  * Hu et al. (2020) Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec.  Open graph benchmark: Datasets for machine learning on graphs.  In _NeurIPS_ , 2020.
  * Karhadkar et al. (2023) Kedar Karhadkar, Pradeep Kr Banerjee, and Guido Montúfar.  FoSR: First-order spectral rewiring for addressing oversquashing in gnns.  In _ICLR_ , 2023.
  * Kingma and Ba (2015) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, _ICLR_ , 2015.
  * Kipf and Welling (2017) Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _ICLR_ , 2017.
  * Klicpera et al. (2018) Johannes Klicpera, Aleksandar Bojchevski, and Stephan Günnemann.  Predict then propagate: Graph neural networks meet personalized pagerank.  In _ICLR_ , 2018.
  * Krizhevsky et al. (2017) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  _Communications of the ACM_ , 60(6):84–90, 2017.
  * Levie et al. (2019) Ron Levie, Elvin Isufi, and Gitta Kutyniok.  On the transferability of spectral graph filters.  In _2019 13th International conference on Sampling Theory and Applications_ , pages 1–5, 2019.
  * Li et al. (2009) Guo-Dong Li, Shiro Masuda, Daisuke Yamaguchi, and Masatake Nagai.  The optimal gnn-pid control system using particle swarm optimization algorithm.  _International Journal of Innovative Computing, Information and Control_ , 5(10):3457–3469, 2009.
  * Li et al. (2018) Qimai Li, Zhichao Han, and Xiao-Ming Wu.  Deeper insights into graph convolutional networks for semi-supervised learning.  In _AAAI_ , 2018.
  * Liu et al. (2021a) Xiaorui Liu, Jiayuan Ding, Wei Jin, Han Xu, Yao Ma, Zitao Liu, and Jiliang Tang.  Graph neural networks with adaptive residual.  In _NeurIPS_ , 2021a.
  * Liu et al. (2021b) Xiaorui Liu, Wei Jin, Yao Ma, Yaxin Li, Hua Liu, Yiqi Wang, Ming Yan, and Jiliang Tang.  Elastic graph neural networks.  In _ICML_ , 2021b.
  * Mallat (2012) Stéphane Mallat.  Group invariant scattering.  _Communications on Pure and Applied Mathematics_ , 65(10):1331–1398, 2012.
  * Maskey et al. (2022) Sohir Maskey, Ron Levie, Yunseok Lee, and Gitta Kutyniok.  Generalization analysis of message passing neural networks on large random graphs.  In _NeurIPS_ , volume 35, pages 4805–4817, 2022.
  * Min et al. (2021) Shengjie Min, Zhan Gao, Jing Peng, Liang Wang, Ke Qin, and Bo Fang.  Stgsn—a spatial–temporal graph neural network framework for time-evolving social networks.  _Knowledge-Based Systems_ , 214:106746, 2021.
  * Min et al. (2020) Yimeng Min, Frederik Wenkel, and Guy Wolf.  Scattering GCN: Overcoming oversmoothness in graph convolutional networks.  _NeurIPS_ , 33:14498–14508, 2020.
  * Morris et al. (2019) Christopher Morris, Martin Ritzert, Matthias Fey, William L Hamilton, Jan Eric Lenssen, Gaurav Rattan, and Martin Grohe.  Weisfeiler and Leman go neural: Higher-order graph neural networks.  In _AAAI_ , volume 33, pages 4602–4609, 2019.
  * Nguyen et al. (2022) Dai Quoc Nguyen, Tu Dinh Nguyen, and Dinh Phung.  Universal graph transformer self-attention networks.  In _The Web Conference_ , 2022.
  * Nguyen et al. (2023) Khang Nguyen, Nong Minh Hieu, Vinh Duc Nguyen, Nhat Ho, Stanley Osher, and Tan Minh Nguyen.  Revisiting over-smoothing and over-squashing using Ollivier-Ricci curvature.  In _ICML_ , pages 25956–25979, 2023.
  * Oono and Suzuki (2019) Kenta Oono and Taiji Suzuki.  Graph neural networks exponentially lose expressive power for node classification.  In _ICLR_ , 2019.
  * Oyallon et al. (2017) Edouard Oyallon, Eugene Belilovsky, and Sergey Zagoruyko.  Scaling the scattering transform: Deep hybrid networks.  In _ICCV_ , pages 5618–5627, 2017.
  * Rampášek et al. (2022) Ladislav Rampášek, Michael Galkin, Vijay Prakash Dwivedi, Anh Tuan Luu, Guy Wolf, and Dominique Beaini.  Recipe for a general, powerful, scalable graph transformer.  In _NeurIPS_ , volume 35, pages 14501–14515, 2022.
  * Ruiz et al. (2021) Luana Ruiz, Fernando Gama, and Alejandro Ribeiro.  Graph neural networks: Architectures, stability, and transferability.  _Proceedings of the IEEE_ , 109(5):660–682, 2021.
  * Topping et al. (2022) Jake Topping, Francesco Di Giovanni, Benjamin Paul Chamberlain, Xiaowen Dong, and Michael M Bronstein.  Understanding over-squashing and bottlenecks on graphs via curvature.  In _ICLR_ , 2022.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In _NeurIPS_ , volume 30, 2017.
  * Veličković et al. (2018) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  In _ICLR_ , 2018.
  * Verma and Zhang (2019) Saurabh Verma and Zhi-Li Zhang.  Stability and generalization of graph convolutional neural networks.  In _KDD_ , 2019.
  * Wang and Zhuang (2019) Yu Guang Wang and Xiaosheng Zhuang.  Tight framelets on graphs for multiscale data analysis.  In _Wavelets and Sparsity XVIII_ , volume 11138, page 111380B, 2019.
  * Wang et al. (2023) Yuelin Wang, Kai Yi, Xinliang Liu, Yu Guang Wang, and Shi Jin.  ACMP: Allen-Cahn message passing with attractive and repulsive forces for graph neural networks.  In _ICLR_ , 2023.
  * Wu et al. (2022) Shiwen Wu, Fei Sun, Wentao Zhang, Xu Xie, and Bin Cui.  Graph neural networks in recommender systems: a survey.  _ACM Computing Surveys_ , 55(5):1–37, 2022.
  * Wu et al. (2021) Zhanghao Wu, Paras Jain, Matthew Wright, Azalia Mirhoseini, Joseph E Gonzalez, and Ion Stoica.  Representing long-range context for graph neural networks with global attention.  In _NeurIPS_ , volume 34, pages 13266–13279, 2021.
  * Xie et al. (2020) Yu Xie, Yuanqiao Zhang, Maoguo Gong, Zedong Tang, and Chao Han.  Mgat: Multi-view graph attention networks.  _Neural Networks_ , 132:180–189, 2020.
  * Xiong et al. (2023) Yuehan Xiong, Wenrui Dai, Wen Fei, Shaohui Li, and Hongkai Xiong.  Anisotropic spherical scattering networks via directional spin wavelet.  _IEEE Transactions on Signal Processing_ , 2023.
  * Xu et al. (2018) Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In _ICLR_ , 2018.
  * Yehudai et al. (2021) Gilad Yehudai, Ethan Fetaya, Eli Meirom, Gal Chechik, and Haggai Maron.  From local structures to size generalization in graph neural networks.  In _ICML_ , pages 11975–11986, 2021.
  * Ying et al. (2021) Chengxuan Ying, Tianle Cai, Shengjie Luo, Shuxin Zheng, Guolin Ke, Di He, Yanming Shen, and Tie-Yan Liu.  Do transformers really perform badly for graph representation?  In _NeurIPS_ , volume 34, pages 28877–28888, 2021.
  * Zhang et al. (2022) Ke Zhang, Xinyan Pu, Jiaxing Li, Jiasong Wu, Huazhong Shu, and Youyong Kong.  Hierarchical diffusion scattering graph neural network.  In _IJCAI_ , volume 7, pages 3737–3743, 2022.
  * Zhang and Zitnik (2020) Xiang Zhang and Marinka Zitnik.  Gnnguard: Defending graph neural networks against adversarial attacks.  In _NeurIPS_ , 2020.
  * Zheng et al. (2021) Xuebin Zheng, Bingxin Zhou, Junbin Gao, Yu Guang Wang, Pietro Lio, Ming Li, and Guido Montúfar.  How framelets enhance graph neural networks.  In _ICML_ , 2021.
  * Zhou (2014) Houqing Zhou.  On some trace inequalities for positive definite Hermitian matrices.  _Journal of Inequalities and Applications_ , 2014(1):1–6, 2014.
  * Zou and Lerman (2019) Dongmian Zou and Gilad Lerman.  Encoding robust representation for graph generation.  In _International Joint Conference on Neural Networks_ , pages 1–9, 2019.
  * Zou and Lerman (2020) Dongmian Zou and Gilad Lerman.  Graph convolutional neural networks via scattering.  _Applied and Computational Harmonic Analysis_ , 49(3):1046–1074, 2020.

##  A Table of notations

Symbol |  Meaning
---|---
$\mathcal{G}$ |  graph
$\mathcal{V}$ |  set of nodes (vertices)
$\mathcal{E}$ |  set of edges
$\mathcal{A}$ |  adjacency matrix
$\mathcal{N}^{m}(i)$ |  $m$-hop neighbors of $i$-th node
$\mathcal{N}(i)$ |  $1$-hop neighbors of $i$-th node
$|\mathcal{V}|$ |  cardinality of set $\mathcal{V}$
$\operatorname{vol}(\mathcal{G})$ |  volume of a graph
$\bm{d}(v)$ |  degree of vertex $v$
$\mathbb{R}$ |  real vector space
$l_{2}(\mathcal{G})$ |  Hilbert space on graph $\mathcal{G}$
$L$ |  graph Laplacian matrix
$\mathcal{L}$ |  normalized graph Laplacian matrix
$\mathcal{D}$ |  degree matrix
$\delta_{\ell,\ell^{\prime}}$ |  Kronecker delta function
$a$ |  low-pass filter
$b^{(r)}$ |  $r$-th high-pass filter
$\widehat{\alpha},\widehat{\beta}$ |  Fourier transformation of $a,b$
$\bm{\eta}$ |  filter bank ( set of filtering functions )
$\Psi$ |  set of scaling functions
$R$ |  dilation scale
$\\{(\lambda_{\ell},{\bm{u}}_{\ell})\\}_{\ell=1}^{n}$ |  eigenpairs of the graph Laplacian $\mathcal{L}$
$\bm{\varphi}_{l,p}(v)$, $\bm{\psi}_{l,p}^{(r)}(v)$ |  low-pass and the $r$-th high-pass framelet basis of node $v$ on level $l$ and node $p$
$\bm{{\mathcal{W}}}_{k,l}$ |  framelet decomposition operator of scale $k$ and level $l$
$\bm{{\mathcal{W}}}_{k,l}^{\natural}$ |  approximated framelet decomposition operator
$\bm{X}^{\text{in }}$ |  input graph feature matrix
$\bm{X}^{\text{out }}$ |  embedding of graph feature matrix
$\hat{\bm{X}}$ |  framelet coefficients
$\Phi_{\mathcal{G}}({\bm{X}})$ |  graph scattering transformation of $\mathcal{G}$
$\phi_{m}(\mathcal{G},{\bm{X}})$ |  $m$-th layer of scattering transformation
$\varphi$ |  average operator
$\hat{\sigma}$ |  pointwise absolute value function
$\sigma$ |  pointwise ReLU function
$\square(\cdot)$ |  aggregation function
$\Gamma(\cdot)$, $\phi(\cdot)$ |  linear functions
$E_{\mathcal{G}}({\bm{X}})$ |  Dirichlet energy of graph $\mathcal{G}$ with feature matrix ${\bm{X}}$
$h(\mathcal{G})$ |  Cheeger constant
$\partial\mathcal{S}$ |  edge boundary of subset $\mathcal{S}$
$\bm{Z}^{(t)}$ |  scattering coefficients at layer $t$
$\bm{\theta}_{0},\bm{\theta}_{r,l}$ |  learnable parameter matrices
$\mathrm{tr}({\bm{X}})$ |  trace of matrix ${\bm{X}}$
$\mu_{i}(\mathcal{A})$ |  $i$-th largest eigenvalue of matrix $\mathcal{A}$
$\mathcal{R}(\mathcal{G})$ |  rewired graph of $\mathcal{G}$