Overcoming Oversmoothness in Graph Convolutional Networks via Hybrid Scattering Networks

Let $G=(V,E,w)$ be a weighted graph, characterized by a set of nodes (also called vertices) $V\coloneqq\{v_{1},\dots,v_{n}\}$, a set of undirected edges $E\subseteq V\times V$, and a function $w:E\to(0,\infty)$ assigning positive edge weights to the edges. Let $\operatorname{\boldsymbol{X}}\in\operatorname{\mathbb{R}}^{n\times d_{X}}$ be a node features matrix. We shall interpret the $i^{\text{th}}$ row of $\operatorname{\boldsymbol{X}}$ as representing the features of the node $v_{i}$, and therefore, we shall denote these rows by either $\operatorname{\boldsymbol{X}}[v_{i}]\in\operatorname{\mathbb{R}}^{1\times d_{X}}$ or $\operatorname{\boldsymbol{X}}_{v_{i}}\in\operatorname{\mathbb{R}}^{1\times d_{X}}$. The columns of $\operatorname{\boldsymbol{X}}$, on the other hand, will be denoted by $\operatorname{\boldsymbol{x}}_{i}$. Each of these columns may be naturally identified with a graph signal, i.e., a function $x_{i}:V\rightarrow\mathbb{R}$, $x_{i}(v_{j})\coloneqq\operatorname{\boldsymbol{x}}_{i}[j]$. In what follows, for simplicity, we will not distinguish between the vectors $\operatorname{\boldsymbol{x}}_{i}$ and the functions $x_{i}$ and will refer to both as graph signals.

We define the weighted adjacency matrix $\operatorname{\boldsymbol{W}}\in\operatorname{\mathbb{R}}^{n\times n}$ of $G$ by $\operatorname{\boldsymbol{W}}[v_{i},v_{j}]\coloneqq w(v_{i},v_{j})$ if $\{v_{i},v_{j}\}\in E$, and set all other entries to zero. We further define the degree matrix $\operatorname{\boldsymbol{D}}\in\operatorname{\mathbb{R}}^{n\times n}$ as the diagonal matrix $\operatorname{\boldsymbol{D}}\coloneqq\operatorname{diag}(\operatorname{deg}(v_{1}),\dots,\operatorname{deg}(v_{n}))$ with each diagonal element $\operatorname{deg}(v_{i})\coloneqq\sum_{j=1}^{n}\operatorname{\boldsymbol{W}}[v_{i},v_{j}]$ being the degree of a node $v_{i}$. In the following, we will also use the shorthand $d_{v_{i}}$ to denote the degree of $v_{i}$. We consider the combinatorial graph Laplacian matrix $\operatorname{\boldsymbol{L}}\coloneqq\operatorname{\boldsymbol{D}}-\operatorname{\boldsymbol{W}}$ and the symmetric normalized Laplacian given by

It is well known that this matrix is symmetric, positive semi-definite, and admits an orthonormal basis of eigenvectors such that $\operatorname{\boldsymbol{\mathcal{L}}}\operatorname{\boldsymbol{q}}_{i}=\lambda_{i}\operatorname{\boldsymbol{q}}_{i}.$ Therefore, we may write

where $\boldsymbol{\Lambda}\coloneqq\operatorname{diag}(\lambda_{1},\dots,\lambda_{n})$ and $\operatorname{\boldsymbol{Q}}$ is the orthogonal matrix whose $i$-th column is $\operatorname{\boldsymbol{q}}_{i}$.

We will use the eigendecomposition of $\operatorname{\boldsymbol{\mathcal{L}}}$ to define the graph Fourier transform, with the eigenvectors $\operatorname{\boldsymbol{q}}_{1},\dots,\operatorname{\boldsymbol{q}}_{n}$ being interpreted as Fourier modes. The notion of oscillation on irregular domains like graphs is delicate, but can be reframed in terms of increasing variation of the modes, with the eigenvalues $0=\lambda_{1}\leqslant\lambda_{2}\leqslant\dots\leqslant\lambda_{n}\leqslant 2$ interpreted as (squared) frequencies.²²2This interpretation originates from motivating the graph Fourier transform via the combinatorial graph Laplacian $\operatorname{\boldsymbol{L}}=\sum_{i=1}^{n}\operatorname{\tilde{\lambda}}_{i}\operatorname{\tilde{\operatorname{\boldsymbol{q}}}}_{i}\operatorname{\tilde{\operatorname{\boldsymbol{q}}}}_{i}^{T}$ with $\operatorname{\tilde{\lambda}}_{i}=\operatorname{\tilde{\operatorname{\boldsymbol{q}}}}_{i}^{T}\operatorname{\boldsymbol{L}}\operatorname{\tilde{\operatorname{\boldsymbol{q}}}}_{i}=\sum_{i=1}^{n}\operatorname{\boldsymbol{W}}[i,j](\operatorname{\tilde{\operatorname{\boldsymbol{q}}}}_{i}[i]-\operatorname{\tilde{\operatorname{\boldsymbol{q}}}}_{i}[j])^{2}$ the variation of $\operatorname{\tilde{\operatorname{\boldsymbol{q}}}}_{i}$. The Fourier transform of a graph signal $\operatorname{\boldsymbol{x}}\in\operatorname{\mathbb{R}}^{n}$ is defined by $\boldsymbol{\hat{x}}[i]\coloneqq\langle\operatorname{\boldsymbol{x}},\operatorname{\boldsymbol{q}}_{i}\rangle$ for $1\leq i\leq n$ and the inverse Fourier transform may be computed by $\operatorname{\boldsymbol{x}}=\sum_{i=1}^{n}\boldsymbol{\hat{x}}[i]\operatorname{\boldsymbol{q}}_{i}$. It will frequently be convenient to write these equations in matrix form as $\boldsymbol{\hat{x}}=\operatorname{\boldsymbol{Q}}^{T}\operatorname{\boldsymbol{x}}$ and $\operatorname{\boldsymbol{x}}=\operatorname{\boldsymbol{Q}}\boldsymbol{\hat{x}}$.

We recall that in the Euclidean setting, the convolution of two signals in the spatial domain corresponds to the pointwise multiplication of their Fourier transforms. Therefore, we may define the convolution of a signal $\operatorname{\boldsymbol{x}}$ with a filter $\operatorname{\boldsymbol{g}}$ by the rule that $\operatorname{\boldsymbol{g}}\star\operatorname{\boldsymbol{x}}$ is the unique vector such that $(\widehat{\operatorname{\boldsymbol{g}}\star\operatorname{\boldsymbol{x}}})[i]=\boldsymbol{\hat{g}}[i]\boldsymbol{\hat{x}}[i]$ for $1\leq i\leq n$. Applying the inverse Fourier transform, one may verify that

where $\boldsymbol{\widehat{G}}\coloneqq\operatorname{diag}(\boldsymbol{\hat{g}})=\operatorname{diag}(\boldsymbol{\hat{g}}[1],\dots,\boldsymbol{\hat{g}}[n])$. Hence, convolutional graph filters can be parameterized by considering the Fourier coefficients in $\boldsymbol{\widehat{G}}$.

A graph filter is a function $F:\mathbb{R}^{n\times d_{X}}\rightarrow\mathbb{R}^{n\times d_{Y}}$ that transforms a node feature matrix $\operatorname{\boldsymbol{X}}\in\mathbb{R}^{n\times d_{X}}$ into a new feature matrix $\operatorname{\boldsymbol{Y}}\in\mathbb{R}^{n\times d_{Y}}$. GNNs typically feature several layers each of which produces a new set of features by filtering the output of the previous layer. We will usually let $\operatorname{\boldsymbol{X}}^{0}$ denote the initial node feature matrix, which is the input to the network, and let $\operatorname{\boldsymbol{X}}^{\ell}$ denote the node feature matrix after the $\ell$-th layer.

In light of Eq. 1, a natural way to construct learned graph filters would be to directly learn the Fourier coefficients in $\boldsymbol{\hat{g}}$. Indeed this was the approach used in the pioneering work of Bruna et al. [27]. However, this approach has several drawbacks. First, it results in $n$ learnable parameters in each convolutional filter. Therefore, a network using such filters would not scale well to large data sets due to the computational cost. At the same time, such filters are not necessarily well-localized in space and are prone to overfitting [28]. Moreover, networks of the form introduced in [27] typically cannot be generalized to different graphs [1]. However, recent work [29] has shown that this latter issue can be overcome by formulating the Fourier coefficients $\boldsymbol{\hat{g}}[i]$ as smooth functions of the Laplacian eigenvalues $\lambda_{i}$, $1\leq i\leq n$. In particular, this will be the case for the filters used in the networks considered in this work.

A common approach (e.g., used in [28, 5, 30, 31, 18]) to formulate such filters is by using polynomials of the Laplacian eigenvalues to set $\boldsymbol{\hat{g}}[i]\coloneqq\sum_{j=0}^{k}\theta_{j}\lambda_{i}^{j}$ (or equivalently $\boldsymbol{\widehat{G}}=\sum_{j=0}^{k}\theta_{j}\boldsymbol{\Lambda}^{j}$) for some $k\ll n$. It can be verified [28] that this approach yields convolutional filters that are $k$-localized in space and that can be written as $\operatorname{\boldsymbol{g}}_{\operatorname{\boldsymbol{\theta}}}\star\operatorname{\boldsymbol{x}}=\sum_{j=0}^{k}\theta_{j}\operatorname{\boldsymbol{\mathcal{L}}}^{j}\operatorname{\boldsymbol{x}}.$ This reduces the number of trainable parameters in each filter from $n$ to $k+1$ and allows one to perform convolution without explicitly computing the spectral decomposition of the Laplacian, which is expensive for large graphs.

One particularly noteworthy network that uses this method is [28], which writes the filters in terms of the Chebyshev polynomials defined by $T_{0}(x)=1$, $T_{1}(x)=x$ and $T_{j}(x)\coloneqq 2xT_{j-1}(x)-T_{j-2}(x)$. They first renormalize the eigenvalue matrix $\boldsymbol{\tilde{\Lambda}}\coloneqq 2\operatorname{\boldsymbol{\Lambda}}/\lambda_{n}-\operatorname{\boldsymbol{I}}_{n}$ and then define $\boldsymbol{\widehat{G}}\coloneqq\sum_{j=0}^{k}\theta_{j}T_{j}(\boldsymbol{\tilde{\Lambda}})$. Letting $\boldsymbol{\tilde{\mathcal{L}}}\coloneqq 2\boldsymbol{\mathcal{L}}/\lambda_{n}-\operatorname{\boldsymbol{I}}_{n}$, yields

One of the most widely used GNNs is the Graph Convolutional Network (GCN) [5]. This network is derived from the Chebyshev construction[28] mentioned above by setting $k=1$ in Eq. 2 and approximating $\lambda_{n}\approx 2$, which yields

To further reduce the number of trainable parameters, the authors then set $\theta\coloneqq\theta_{0}=-\theta_{1}$. The resulting convolutional filter has the form

One may verify that $\operatorname{\boldsymbol{I}}_{n}+\boldsymbol{D}^{-1/2}\boldsymbol{W}\boldsymbol{D}^{-1/2}=2\operatorname{\boldsymbol{I}}_{n}-\boldsymbol{\mathcal{L}}$, and therefore, Eq. 3 essentially corresponds to setting $\boldsymbol{\hat{g}}[i]\coloneqq\theta(2-\lambda_{i})$ in Eq. 1. The eigenvalues of $\operatorname{\boldsymbol{I}}_{n}+\boldsymbol{D}^{-1/2}\boldsymbol{W}\boldsymbol{D}^{-1/2}$ take values in $[0,2].$ Thus, to avoid vanishing or exploding gradients, the authors use a renormalization trick

where $\boldsymbol{\tilde{W}}\coloneqq\operatorname{\boldsymbol{I}}_{n}+\boldsymbol{W}$ and $\boldsymbol{\tilde{D}}$ is a diagonal matrix with $\boldsymbol{\tilde{D}}[v_{i},v_{i}]\coloneqq\sum_{j=1}^{n}\boldsymbol{\tilde{W}}[v_{i},v_{j}]$ for $i\in[n]$. Setting $\boldsymbol{A}\coloneqq\boldsymbol{\tilde{D}}^{-1/2}\boldsymbol{\tilde{W}}\boldsymbol{\tilde{D}}^{-1/2},$ and using multiple channels we obtain a layer-wise propagation rule of the form $\operatorname{\boldsymbol{x}}_{j}^{\ell}=\sigma\big{(}\sum_{i=1}^{N_{\ell-1}}\theta_{ij}^{\ell}\operatorname{\boldsymbol{A}}\operatorname{\boldsymbol{x}}_{i}^{\ell-1}\big{)},$ where $N_{\ell}$ is the number of channels used in the $\ell$-th layer and $\sigma(\cdot)$ is an elementwise nonlinearity. In matrix form we write

We interpret the matrix $\operatorname{\boldsymbol{A}}$ as computing a localized average of each channel $\operatorname{\boldsymbol{x}}_{i}^{\ell-1}$ around each mode and the matrix $\operatorname{\boldsymbol{\Theta}}$ as sharing information across channels. This filter can also be observed at the node level as

where $\operatorname{\mathcal{N}}_{v}$ denotes the one-step neighborhood of node $v$ and $\operatorname{\mathcal{N}}_{\underline{v}}\coloneqq\operatorname{\mathcal{N}}_{v}\cup\{v\}$. This process can be split into three steps:

which we refer to as the transformation step (Eq. 6a), the aggregation step (Eq. 6b) and the activation step (Eq. 6c).

As discussed earlier, the GCN filter described above may be viewed as a low-pass filter that suppresses high-frequency information. For simplicity, we focus on the convolution in Eq. 3 before the renormalization. This convolution essentially corresponds to point-wise Fourier multiplication by $\boldsymbol{\hat{g}}[i]=\theta(2-\lambda_{i})$, which is strictly decreasing in $\lambda_{i}$. Therefore, repeated applications of this filter effectively zero-out the higher frequencies. This is consistent with the oversmoothing problem discussed in [9].

Another popular network that is widely used for node classification tasks is the graph attention network (GAT) [17], which uses an attention mechanism to guide and adjust the aggregation of features from adjacent nodes. First, the node features $\operatorname{\boldsymbol{X}}^{\ell-1}$ are linearly transformed to $\boldsymbol{\bar{X}}^{\ell}\coloneqq\operatorname{\boldsymbol{X}}^{\ell-1}\operatorname{\boldsymbol{\Theta}}$ using a learned weight matrix $\operatorname{\boldsymbol{\Theta}}\in\mathbb{R}^{d\times d^{\prime}}$. Then, the aggregation coefficients are learned via

where $\operatorname{\boldsymbol{a}}\in\operatorname{\mathbb{R}}^{2d^{\prime}}$ is a shared attention vector and $\|$ denotes horizontal concatenation. The output feature corresponding to a single attention head is given by $\operatorname{\boldsymbol{X}}^{\ell}[v]=\sigma(\sum_{u\in\operatorname{\mathcal{N}}_{v}}\alpha_{v\leftarrow u}\boldsymbol{\bar{X}}^{\ell}[u])$. To increase the expressivity of this network, the authors then use a multi-headed attention mechanism, with $\Gamma$ heads, to generate concatenated features

Many GNN models, including GCN [5] and GAT [17], are subject to the so-called oversmoothing problem [9], caused by aggregation steps (such as Eq. 6b) that essentially consist of localized averaging operations. As discussed in Sec. 3.3 and also [10], from a signal processing point of view, this corresponds to a low-pass filtering of graph signals. Moreover, as discussed in [32], these networks are also subject to underreaching. Most GNNs (including GCN and GAT) can only relate information from nodes within a distance equal to the number of GNN layers, and because of the oversmoothing problem, they typically use a small number of layers in practice. Therefore, the oversmoothing and underreaching problems combine to significantly limit the ability of GNNs to capture long-range dependencies. In Sec. 4, we will introduce a hybrid network, which aims to address these challenges by using both GCN-style channels and channels based on the geometric scattering transform discussed below.

In this section, we review the geometric scattering transform constructed in [12] for graph classification and show how it may be adapted for node-level tasks. As we shall see, this node-level geometric scattering will address the challenges discussed above in Sec. 3.5, by using band-pass filters that capture high-frequency information and have wider receptive fields.

The geometric scattering transform uses wavelets based upon raising the lazy random walk matrix

to dyadic powers $2^{j}$, which can be interpreted as differing degrees of resolution. The right eigenvectors and eigenvalues of $\operatorname{\boldsymbol{P}}$ are $\operatorname{\boldsymbol{u}}_{i}\coloneqq\operatorname{\boldsymbol{D}}^{1/2}\operatorname{\boldsymbol{q}}_{i}$ and $\omega_{i}\coloneqq 1-\lambda_{i}/2,1\leq i\leq n$, respectively. Entrywise, we note that

Thus, $\operatorname{\boldsymbol{P}}$ may be viewed as a localized averaging operation operator analogous to those used in, e.g., GCN, and the powers $\operatorname{\boldsymbol{P}}^{2^{j}}$ may be viewed as low-pass filters which suppress high-frequencies. In order to better retain this high-frequency information, we define multiscale graph diffusion wavelets by subtracting these low-pass filters at different scales [16]. Specifically, for $k\in\mathbb{N}_{0}$, we define a wavelet $\boldsymbol{\Psi}_{k}\in\mathbb{R}^{n\times n}$ at scale $2^{k}$ by

We may interpret each $\boldsymbol{\Psi}_{k}$ as capturing information at a different frequency band. From a spatial perspective, we may view $\operatorname{\boldsymbol{\Psi}}_{k}$ as encoding information on how a $2^{k}$-step neighborhood differs from a smaller $2^{k-1}$-step one. Such wavelets are usually organized in a filter bank

along with a low-pass filter $\boldsymbol{\Phi}_{K}\coloneqq\boldsymbol{P}^{2^{K}}$. Proposition 2.2 of  [22] (restated in here as Proposition 1) shows that this filter bank is a self-adjoint, non-expansive frame on a suitably weighted inner-product space.

The geometric scattering transform is a multi-layered architecture that iteratively applies wavelet convolutions and nonlinear activation functions in an alternating cascade as illustrated in Fig. 1. It is parameterized by paths $p\coloneqq(k_{1},\dots,k_{m})$. Formally, we define

where $\sigma$ is a nonlinear activation function.³³3In a slight deviation from previous work, here $\boldsymbol{U}_{p}$ does not include the outermost nonlinearity in the cascade. We note that in other work focusing on graph-classification, the scattering features are frequently aggregated into $q^{\text{th}}$-order moments,

We also note that the nonlinearity $\sigma$ might vary in each step of the cascade. However, we will suppress this possible dependence to avoid cumbersome notation. We also note that in our theoretical results, if we assume, for example, that $\sigma$ is strictly increasing, this assumption is intended to apply to all nonlinearities used in the cascade. In our experiments, we use the absolute value, i.e., $\sigma(\operatorname{\boldsymbol{x}}[v_{i}])=\lvert\operatorname{\boldsymbol{x}}[v_{i}]\rvert$.

The original formulations of geometric scattering were fully designed networks without any learned convolutions between channels. Here, we will incorporate learning by defining the following scattering propagation rule similar to the one used by GCN:

Analogously to GCN, we note that for $v\in V$, the scattering propagation rule can be decomposed into three steps:

Importantly, we note that the scattering transform addresses the underreaching problem as wavelets $\boldsymbol{\Psi}_{k}=\operatorname{\boldsymbol{P}}^{2^{k-1}}-\operatorname{\boldsymbol{P}}^{2^{k}}$ that are leveraged in the aggregation step (Eq. 14b) have larger receptive fields than most traditional GNNs. However, scattering does not result in oversmoothing because the subtraction results in band-pass filters rather than low-pass filters. In this manner, the scattering transform addresses the challenges discussed in the previous subsection.

The Scattering GCN (Sc-GCN) network was first introduced in [33]. Here, the aggregation module concatenates the filter responses horizontally yielding wide node representations $\operatorname{\boldsymbol{X}}^{\ell}$ of the form

Sc-GCN then learns relationships between the channels via the graph residual convolution.

The paths $p$ and the parameter $\alpha$ from the graph residual convolution are tuned as hyperparameters of the network.

An important observation in Sc-GCN above is that the model decides globally about how to combine different channel information. The network first concatenates all of the features from the low-pass and band-pass channels in Eq. 18 and then combines these features via multiplication with the weight matrix $\operatorname{\boldsymbol{\Theta}}_{\text{res}}$. However, for complex tasks or datasets, important regularity patterns may vary significantly across different graph regions. In such settings, a model should ideally attend locally over the aggregation and adaptively combine filter responses at different nodes.

This observation inspired the design of the Geometric Scattering Attention Network (GSAN) [34]. Drawing inspiration from  [35], GSAN uses an aggregation module based on a multi-head node-level attention framework. However, the attention mechanism used in GSAN differs from [35] by attending over the combination of the different filter responses rather than over the combination of node features from neighboring nodes. We will first focus on the processing performed independently by each attention head and then discuss a multi-head configuration.

In a slight deviation from Sc-GCN, the weight matrices in Eq. 15 and Eq. 16 are shared across the filters of each attention head. Therefore, both aggregation steps (Eq. 6b and Eq. 14b) take the same input denoted by $\boldsymbol{\bar{X}}^{\ell}\coloneqq\operatorname{\boldsymbol{X}}^{\ell-1}\operatorname{\boldsymbol{\Theta}}^{\ell}$. Next, we define $\boldsymbol{\bar{X}}_{\text{low},i}\coloneqq\operatorname{\boldsymbol{A}}^{r_{i}}\boldsymbol{\bar{X}}^{\ell}$ and $\boldsymbol{\bar{X}}_{\text{band},i}\coloneqq|\operatorname{\boldsymbol{U}}_{p_{i}}(\boldsymbol{\bar{X}}^{\ell})|$ to be the outputs of the aggregation steps, Eq. 6b and Eq. 14b, with the bias terms set to zero. We then compute attention score vectors $\operatorname{\boldsymbol{e}}_{\text{low},i}^{\ell},\operatorname{\boldsymbol{e}}_{\text{band},i}^{\ell}\in\operatorname{\mathbb{R}}^{n}$ that will be used to determine the importance of each of the filter responses for every single node. We calculate

with analogous $\operatorname{\boldsymbol{e}}_{\text{band},i}^{\ell}$ and $\operatorname{\boldsymbol{a}}\in\operatorname{\mathbb{R}}^{2d_{\ell}}$ being a learned attention vector that is shared across all filters of the attention head. These attention scores are then normalized across all filters using the softmax function. Specifically, we define

where the exponential function is applied elementwise, and define $\operatorname{\boldsymbol{\alpha}}_{\text{band},j}^{\ell}$ analogously. Finally, for every node, the filter responses are summed together, weighted by the corresponding (node-to-filter) attention score. We also normalize by the number of filters $C\coloneqq C_{\text{low}}+C_{\text{band}}$, which gives

where $\boldsymbol{\alpha}\odot\operatorname{\boldsymbol{X}}\coloneqq\operatorname{diag}(\boldsymbol{\alpha})\operatorname{\boldsymbol{X}}$ denotes the Hadamard (elementwise) product of $\boldsymbol{\alpha}$ with each column of $\operatorname{\boldsymbol{X}}$. We further use $\tilde{\sigma}(\cdot)=\text{ReLU}(\cdot)$ in the equation above.