Diffusion Improves Graph Learning

Introduction

When people started using graphs for evaluating chess tournaments in the middle of the 19th century they only considered each player’s direct opponents, i.e. their first-hop neighbors. Only later was the analysis extended to recursively consider higher-order relationships via ${\bm{A}}^{2}$, ${\bm{A}}^{3}$, etc. and finally generalized to consider all exponents at once, using the adjacency matrix’s dominant eigenvector . The field of Graph Neural Networks (GNNs) is currently in a similar state. Graph Convolutional Networks (GCNs) , also referred to as Message Passing Neural Networks (MPNNs) are the prevalent approach in this field but they only pass messages between neighboring nodes in each layer. These messages are then aggregated at each node to form the embedding for the next layer. While MPNNs do leverage higher-order neighborhoods in deeper layers, limiting each layer’s messages to one-hop neighbors seems arbitrary. Edges in real graphs are often noisy or defined using an arbitrary threshold , so we can clearly improve upon this approach.

Since MPNNs only use the immediate neigborhod information, they are often referred to as spatial methods. On the other hand, spectral-based models do not just rely on first-hop neighbors and capture more complex graph properties . However, while being theoretically more elegant, these methods are routinely outperformed by MPNNs on graph-related tasks and do not generalize to previously unseen graphs. This shows that message passing is a powerful framework worth extending upon. To reconcile these two separate approaches and combine their strengths we propose a novel technique of performing message passing inspired by spectral methods: Graph diffusion convolution (GDC). Instead of aggregating information only from the first-hop neighbors, GDC aggregates information from a larger neighborhood. This neighborhood is constructed via a new graph generated by sparsifying a generalized form of graph diffusion. We show how graph diffusion is expressed as an equivalent polynomial filter and how GDC is closely related to spectral-based models while addressing their shortcomings. GDC is spatially localized, scalable, can be combined with message passing, and generalizes to unseen graphs. Furthermore, since GDC generates a new sparse graph it is not limited to MPNNs and can trivially be combined with any existing graph-based model or algorithm in a plug-and-play manner, i.e. without requiring changing the model or affecting its computational complexity. We show that GDC consistently improves performance across a wide range of models on both supervised and unsupervised tasks and various homophilic datasets. In summary, this paper’s core contributions are:

Proposing graph diffusion convolution (GDC), a more powerful and general, yet spatially localized alternative to message passing that uses a sparsified generalized form of graph diffusion. GDC is not limited to GNNs and can be combined with any graph-based model or algorithm.

Analyzing the spectral properties of GDC and graph diffusion. We show how graph diffusion is expressed as an equivalent polynomial filter and analyze GDC’s effect on the graph spectrum.

Comparing and evaluating several specific variants of GDC and demonstrating its wide applicability to supervised and unsupervised learning on graphs.

Generalized graph diffusion

with the weighting coefficients $\theta_{k}$, and the generalized transition matrix ${\bm{T}}$. The choice of $\theta_{k}$ and ${\bm{T}}^{k}$ must at least ensure that Eq. 1 converges. In this work we will consider somewhat stricter conditions and require that $\sum_{k=0}^{\infty}\theta_{k}=1$, $\theta_{k}\in$, and that the eigenvalues of ${\bm{T}}$ are bounded by $\lambda_{i}\in$, which together are sufficient to guarantee convergence. Note that regular graph diffusion commonly requires ${\bm{T}}$ to be column- or row-stochastic.

Weighting coefficients. We compute the series defined by Eq. 1 either in closed-form, if possible, or by restricting the sum to a finite number $K$. Both the coefficients defined by PPR and the heat kernel give a closed-form solution for this series that we found to perform well for the tasks considered. Note that we are not restricted to using ${\bm{T}}_{\text{rw}}$ and can use any generalized transition matrix along with the coefficients $\theta^{\text{PPR}}_{k}$ or $\theta^{\text{HK}}_{k}$ and the series still converges. We can furthermore choose $\theta_{k}$ by repurposing the graph-specific coefficients obtained by methods that optimize coefficients analogous to $\theta_{k}$ as part of their training process. We investigated this approach using label propagation and node embedding models . However, we found that the simple coefficients defined by PPR or the heat kernel perform better than those learned by these models (see Fig. 7 in Sec. 6).

Graph diffusion convolution

Intuition. The general intuition behind GDC is that graph diffusion smooths out the neighborhood over the graph, acting as a kind of denoising filter similar to Gaussian filters on images. This helps with graph learning since both features and edges in real graphs are often noisy. Previous works also highlighted the effectiveness of graph denoising. Berberidis & Giannakis showed that PPR is able to reconstruct the underlying probability matrix of a sampled stochastic block model (SBM) graph. Kloumann et al. and Ragain showed that PPR is optimal in recovering the SBM and DCSBM clusters in the space of landing probabilities under the mean field assumption. Li et al. generalized this result by analyzing the convergence of landing probabilities to their mean field values. These results confirm the intuition that graph diffusion-based smoothing indeed recovers meaningful neighborhoods from noisy graphs.

Limitations. GDC is based on the assumption of homophily, i.e. “birds of a feather flock together” . Many methods share this assumption and most common datasets adhere to this principle. However, this is an often overlooked limitation and it seems non-straightforward to overcome. One way of extending GDC to heterophily, i.e. “opposites attract”, might be negative edge weights . Furthermore, we suspect that GDC does not perform well in settings with more complex edges (e.g. knowledge graphs) or graph reconstruction tasks such as link prediction. Preliminary experiments showed that GDC indeed does not improve link prediction performance.

Spectral analysis of GDC

Even though GDC is a spatial-based method it can also be interpreted as a graph convolution and analyzed in the graph spectral domain. In this section we show how generalized graph diffusion is expressed as an equivalent polynomial filter and vice versa. Additionally, we perform a spectral analysis of GDC, which highlights the tight connection between GDC and spectral-based models.

Spectral graph theory. To employ the tools of spectral theory to graphs we exchange the regular Laplace operator with either the unnormalized Laplacian ${\bm{L}}_{\text{un}}={\bm{D}}-{\bm{A}}$, the random-walk normalized ${\bm{L}}_{\text{rw}}={\bm{I}}_{N}-{\bm{T}}_{\text{rw}}$, or the symmetric normalized graph Laplacian ${\bm{L}}_{\text{sym}}={\bm{I}}_{N}-{\bm{T}}_{\text{sym}}$ . The Laplacian’s eigendecomposition is ${\bm{L}}={\bm{U}}{\bm{\Lambda}}{\bm{U}}^{T}$, where both ${\bm{U}}$ and ${\bm{\Lambda}}$ are real-valued. The graph Fourier transform of a vector ${\bm{x}}$ is then defined via $\hat{{\bm{x}}}={\bm{U}}^{T}{\bm{x}}$ and its inverse as ${\bm{x}}={\bm{U}}\hat{{\bm{x}}}$. Using this we define a graph convolution on $\mathcal{G}$ as ${\bm{x}}*_{\mathcal{G}}{\bm{y}}={\bm{U}}(({\bm{U}}^{T}{\bm{x}})\odot({\bm{U}}^{T}{\bm{y}}))$, where $\odot$ denotes the Hadamard product. Hence, a filter $g_{\xi}$ with parameters $\xi$ acts on ${\bm{x}}$ as $g_{\xi}({\bm{L}}){\bm{x}}={\bm{U}}\hat{{\bm{G}}}_{\xi}({\bm{\Lambda}}){\bm{U}}^{T}{\bm{x}}$, where $\hat{{\bm{G}}}_{\xi}({\bm{\Lambda}})=\text{diag}(\hat{g}_{\xi,1}({\bm{\Lambda}}),\ldots,\hat{g}_{\xi,N}({\bm{\Lambda}}))$. A common choice for $g_{\xi}$ in the literature is a polynomial filter of order $J$, since it is localized and has a limited number of parameters :

Graph diffusion as a polynomial filter. Comparing Eq. 1 with Eq. 2 shows the close relationship between polynomial filters and generalized graph diffusion since we only need to exchange ${\bm{L}}$ by ${\bm{T}}$ to go from one to the other. To make this relationship more specific and find a direct correspondence between GDC with $\theta_{k}$ and a polynomial filter with parameters $\xi_{j}$ we need to find parameters that solve

To find these parameters we choose the Laplacian corresponding to ${\bm{L}}={\bm{I}}_{n}-{\bm{T}}$, resulting in (see App. A)

which shows the direct correspondence between graph diffusion and spectral methods. Note that we need to set $J=K$. Solving Eq. 4 for the coefficients corresponding to the heat kernel $\theta_{k}^{\text{HK}}$ and PPR $\theta_{k}^{\text{PPR}}$ leads to

showing how the heat kernel and PPR are expressed as polynomial filters. Note that PPR’s corresponding polynomial filter converges only if $\alpha>0.5$. This is caused by changing the order of summation when deriving $\xi^{\text{PPR}}_{j}$, which results in an alternating series. However, if the series does converge it gives the exact same transformation as the equivalent graph diffusion.

2. Sum over ${\bm{T}}^{k}$. Summation does not affect the eigenvectors of the original matrix, since ${\bm{T}}^{k}{\bm{v}}_{i}=\lambda_{i}{\bm{T}}^{k-1}{\bm{v}}_{i}=\lambda_{i}^{k}{\bm{v}}_{i}$, for the eigenvector ${\bm{v}}_{i}$ of ${\bm{T}}$ with associated eigenvalue $\lambda_{i}$. This also shows that the eigenvalues are transformed as

Limitations of spectral-based models. While there are tight connections between GDC and spectral-based models, GDC is actually spatial-based and therefore does not share their limitations. Similar to polynomial filters, GDC does not compute an expensive eigenvalue decomposition, preserves locality on the graph and is not limited to a single graph after training, i.e. typically the same coefficients $\theta_{k}$ can be used across graphs. The choice of coefficients $\theta_{k}$ depends on the type of graph at hand and does not change significantly between similar graphs. Moreover, the hyperparameters $\alpha$ of PPR and $t$ of the heat kernel usually fall within a narrow range that is rather insensitive to both the graph and model (see Fig. 8 in Sec. 6).

Related work

Graph diffusion and random walks have been extensively studied in classical graph learning , especially for clustering , semi-supervised classification , and recommendation systems . For an overview of existing methods see Masuda et al. and Fouss et al. .

The first models similar in structure to current Graph Neural Networks (GNNs) were proposed by Sperduti & Starita and Baskin et al. , and the name GNN first appeared in . However, they only became widely adopted in recent years, when they started to outperform classical models in many graph-related tasks . In general, GNNs are classified into spectral-based models , which are based on the eigendecomposition of the graph Laplacian, and spatial-based methods , which use the graph directly and form new representations by aggregating the representations of a node and its neighbors. However, this distinction is often rather blurry and many models can not be clearly attributed to one type or the other. Deep learning also inspired a variety of unsupervised node embedding methods. Most models use random walks to learn node embeddings in a similar fashion as word2vec and have been shown to implicitly perform a matrix factorization . Other unsupervised models learn Gaussian distributions instead of vectors , use an auto-encoder , or train an encoder by maximizing the mutual information between local and global embeddings .

There have been some isolated efforts of using extended neighborhoods for aggregation in GNNs and graph diffusion for node embeddings. PPNP propagates the node predictions generated by a neural network using personalized PageRank, DCNN extends node features by concatenating features aggregated using the transition matrices of $k$-hop random walks, GraphHeat uses the heat kernel and PAN the transition matrix of maximal entropy random walks to aggregate over nodes in each layer, PinSage uses random walks for neighborhood aggregation, and MixHop concatenates embeddings aggregated using the transition matrices of $k$-hop random walks before each layer. VERSE learns node embeddings by minimizing KL-divergence from the PPR matrix to a low-rank approximation. Attention walk uses a similar loss to jointly optimize the node embeddings and diffusion coefficients $\theta_{k}$. None of these works considered sparsification, generalized graph diffusion, spectral properties, or using preprocessing to generalize across models.

Experimental results

Datasets and models. We evaluate GDC on six datasets: The citation graphs Citeseer , Cora , and PubMed , the co-author graph Coauthor CS , and the co-purchase graphs Amazon Computers and Amazon Photo . We only use their largest connected components. We show how GDC impacts the performance of 9 models: Graph Convolutional Network (GCN) , Graph Attention Network (GAT) , jumping knowledge network (JK) , Graph Isomorphism Network (GIN) , and ARMA are supervised models. The degree-corrected stochastic block model (DCSBM) , spectral clustering (using ${\bm{L}}_{\text{sym}}$) , DeepWalk , and Deep Graph Infomax (DGI) are unsupervised models. Note that DGI uses node features while other unsupervised models do not. We use $k$-means clustering to generate clusters from node embeddings. Dataset statistics and hyperparameters are reported in App. B.

Clustering. We highlight GDC’s ability to be combined with any graph-based model by reporting the performance of a diverse set of models that use a wide range of paradigms. Fig. 4 shows the unsupervised accuracy obtained by matching clusters to ground-truth classes using the Hungarian algorithm. Accuracy consistently and significantly improves for all models and datasets. Note that spectral clustering uses the graph’s eigenvectors, which are not affected by the diffusion step itself. Still, its performance improves by up to $30$ percentage points. Results in tabular form are presented in App. B.2.

In this work we concentrate on node-level prediction tasks in a transductive setting. However, GDC can just as easily be applied to inductive problems or different tasks like graph classification. In our experiments we found promising, yet not as consistent results for graph classification (e.g. $2.5$ percentage points with GCN on the DD dataset ). We found no improvement for the inductive setting on PPI , which is rather unsurprising since the underlying data used for graph construction already includes graph diffusion-like mechanisms (e.g. regulatory interactions, protein complexes, and metabolic enzyme-coupled interactions). We furthermore conducted experiments to answer five important questions:

How to choose the transition matrix ${\bm{T}}$? We found ${\bm{T}}_{\text{sym}}$ to perform best across datasets. More specifically, Fig. 7 shows that the symmetric version on average outperforms the random walk transition matrix ${\bm{T}}_{\text{rw}}$. This figure also shows that GCN accuracy is largely insensitive to self-loops when using ${\bm{T}}_{\text{sym}}$ – all changes lie within the estimated uncertainty. However, we did find that other models, e.g. GAT, perform better with self-loops (not shown).

How to choose the coefficients $\theta_{k}$? We found the coefficients defined by PPR and the heat kernel to be effective choices for $\theta_{k}$. Fig. 8 shows that their optimal hyperparameters typically fall within a narrow range of $\alpha\in[0.05,0.2]$ and $t\in$. We also tried obtaining $\theta_{k}$ from models that learn analogous coefficients . However, we found that $\theta_{k}$ obtained by these models tend to converge to a minimal neighborhood, i.e. they converge to $\theta_{0}\approx 1$ or $\theta_{1}\approx 1$ and all other $\theta_{k}\approx 0$.

This is caused by their training losses almost always decreasing when the considered neighborhood shrinks. We were able to control this overfitting to some degree using strong regularization (specifically, we found $L_{2}$ regularization on the difference of neighboring coefficients $\theta_{k+1}-\theta_{k}$ to perform best). However, this requires hand-tuning the regularization for every dataset, which defeats the purpose of learning the coefficients from the graph. Moreover, we found that even with hand-tuned regularization the coefficients defined by PPR and the heat kernel perform better than trained $\theta_{k}$, as shown in Fig. 7.

Which nodes benefit from GDC? Our experiments showed no correlation of improvement with most common node properties, except for the distance from the training set. Nodes further away from the training set tend to benefit more from GDC, as demonstrated by Fig. 10. Besides smoothing out the neighborhood, GDC also has the effect of increasing the model’s range, since it is no longer restricted to only using first-hop neighbors. Hence, nodes further away from the training set influence the training and later benefit from the improved model weights.

Conclusion

We propose graph diffusion convolution (GDC), a method based on sparsified generalized graph diffusion. GDC is a more powerful, yet spatially localized extension of message passing in GNNs, but able to enhance any graph-based model. We show the tight connection between GDC and spectral-based models and analyzed GDC’s spectral properties. GDC shares many of the strengths of spectral methods and none of their weaknesses. We conduct extensive and rigorous experiments that show that GDC consistently improves the accuracy of a wide range of models on both supervised and unsupervised tasks across various homophilic datasets and requires very little hyperparameter tuning. There are many extensions and applications of GDC that remain to be explored. We expect many graph-based models and tasks to benefit from GDC, e.g. graph classification and regression. Promising extensions include other diffusion coefficients $\theta_{k}$ such as those given by the methods presented in Fouss et al. and more advanced random walks and operators that are not defined by powers of a transition matrix.

This research was supported by the German Federal Ministry of Education and Research (BMBF), grant no. 01IS18036B, and by the Deutsche Forschungsgemeinschaft (DFG) through the Emmy Noether grant GU 1409/2-1 and the TUM International Graduate School of Science and Engineering (IGSSE), GSC 81. The authors of this work take full responsibilities for its content.

References

Appendix A Graph diffusion as a polynomial filter

We want to find a direct correspondence between graph diffusion with $\theta_{k}$ and a polynomial filter with parameters $\xi_{j}$, i.e.

To do so, we first expand ${\bm{T}}={\bm{I}}_{N}-{\bm{L}}$ and use the binomial equation, i.e.

where we recognize the coefficients $\xi_{j}$ and see that we need to set $J=K$. Note that we reordered the summation indices by recognizing the triangular sum, i.e. the sum over index pairs $(j,k)$ with $j\leq k$. The equation for conversion in the opposite direction is obtained in the same way since ${\bm{L}}={\bm{I}}_{N}-{\bm{T}}$. To obtain a more convenient form for $K\to\infty$ we shift the summation index using $m=k-j$, i.e.

To find corresponding coefficients for the heat kernel, we let $K\to\infty$, set $\theta_{k}=e^{-t}\frac{t^{k}}{k!}$, and use the exponential series to obtain

To obtain the coefficients for PPR, we let $K\to\infty$, set $\theta_{k}=\alpha(1-\alpha)^{k}$, and recognize the series expansion $\frac{1}{(1-x)^{j+1}}=\sum_{m=0}^{\infty}\binom{m+j}{m}x^{m}$, resulting in

Appendix B Experiments

For optimizing the hyperparameters for node classification the data is split into a development and a test set. The development set contains 1500 nodes for all datasets but for Coauthor CS, where it contains 5000 nodes. All remaining nodes are part of the test set and only used once for testing. The development set is split into a training set containing 20 nodes per class and a validation set with the remaining nodes. For every run the accuracy is determined using 100 different random splits of the development set using fixed seeds. Different seeds are used for validation and test splits. Early stopping patience is set to 100 epochs with a maximum limit of 10000 epochs, which is never reached. The patience is reset after an increase in accuracy on the validation set. For the test runs we select the hyperparameter configurations that showed the highest average accuracy on the validation splits.

We use the same development set for optimizing the hyperparameters for clustering. The test set is only once for generating test results. Clustering results are averaged over 20 randomly initialized runs.

Confidence intervals are calculated by bootstrapping the accuracy results from 100 or 20 runs, respectively, with 1000 samples. All implementations for node classification as well as DGI are based on PyTorch and PyTorch Geometric . The remaining experiments are based on NumPy , SciPy , graph-tool , and gensim . For $k$-means clustering we use the implementation by scikit-learn . All datasets are included in PyTorch Geometric, available at https://github.com/rusty1s/pytorch_geometric. Experiments using PyTorch are run on Nvidia GPUs using CUDA and the remaining experiments are run on Intel CPUs.

For all experiments the largest connected component of the graph is selected. Dropout probability is set to $p=0.5$ for all experiments and performed after every application of the activation function. PPR preprocessing is done with $\alpha\in[0.05,0.30]$, heat kernel preprocessing with $t\in$. For top-$k$ matrix sparsification $k$ is set to either 64 or 128 and for $\epsilon$-thresholding $\epsilon$ is chosen from $[0.00001,0.01]$. We do not choose $\epsilon$ directly but rather calculate which $\epsilon$ corresponds to a chosen average degree. For node classification we use the Adam optimizer with a learning rate of 0.01. The hidden dimension of GNNs is kept fixed at 64 with the exception of ARMA, where the dimensionality of a single stack is chosen from 16 or 32. For ARMA, up to three stacks and two layers are tested. GCN and GAT are run with up to 4 layers, JK and GIN with up to six layers. $L_{\text{2}}$-regularization is performed on the parameters of the first layer of every model with $\lambda_{L_{\text{1}}}\in[0.001,10]$. Unsupervised models use a node embedding dimension of 128. DGI uses the Adam optimizer with a learning rate of 0.001. For a full list of final hyperparameters per model, diffusion, and dataset see Sec. B.3.

B.2 Results

To support our claim of achieving state-of-the-art node classification performance we also include results (and hyperparameters) of APPNP, which has been shown to be the current state of the art for semi-supervised node classification and uses graph diffusion internally .

B.3 Hyperparameters