Title: UNSUPERVISED HIERARCHICAL GRAPH REPRESEN-TATION LEARNING WITH VARIATIONAL BAYES

Abstract: Hierarchical graph representation learning is an emerging subject owing to the increasingly popular adoption of graph neural networks in machine learning and applications. Loosely speaking, work under this umbrella falls into two categories: (a) use a predefined graph hierarchy to perform pooling; and (b) learn the hierarchy for a given graph through differentiable parameterization of the coarsening process. These approaches are supervised; a predictive task with ground-truth labels is used to drive the learning. In this work, we propose an unsupervised approach, BAYESPOOL, with the use of variational Bayes. It produces graph representations given a predefined hierarchy. Rather than relying on labels, the training signal comes from the evidence lower bound of encoding a graph and decoding the subsequent one in the hierarchy. Node features are treated latent in this variational machinery, so that they are produced as a byproduct and are used in downstream tasks. We demonstrate a comprehensive set of experiments to show the usefulness of the learned representation in the context of graph classification.

1. INTRODUCTION:

Graph representation learning has attracted a surge of interest recently, inspired by the widespread success of representation learning in the image and language domains through the use of deep neural networks for parameterization. A substantial number of graph neural network (GNN) architectures (Bruna et al., 2014; Henaff et al., 2015; Duvenaud et al., 2015; Defferrard et al., 2016; Kipf & Welling, 2017; Hamilton et al., 2017; Chen et al., 2018; Velicković et al., 2018; Ying et al., 2018a; Liao et al., 2019; Xu et al., 2019) extend the convolution filters for a regular grid of data (e.g., image pixels, time series, and sequences) to irregularly connected graph neighborhoods. This extension naturally stimulates the quest of also extending the pooling operation in convolutional neural networks (CNN) to graphs. The challenge lies in the irregular connections as opposed to a regular grid structure, whereby partitioning is straightforward.

Graph pooling is used in at least two scenarios. One is global pooling: it pools the vector representations of all nodes to a single vector as the graph representation. Simple operators such as max or mean are applied. A slightly more complex operator is the weighted sum, wherein the weights are computed through attention (Velicković et al., 2018; Lee et al., 2019) . A recently proposed operator is top-k pooling (Zhang et al., 2018) , whereby a fixed number of node representations at the top of a sorted list is retained so that convolutions or feed-forward transformations are applied.

The second use of pooling is the creation of a graph hierarchy. In this scenario, pooling is local, more similar to that in CNNs. It is interfaced with graph coarsening (also called graph reduction or graph compression), a generic form of which is to cluster the nodes of the original graph into a node of the coarse graph. Then, the representations of the nodes in the cluster are pooled. The clustering may be obtained by using existing graph coarsening or graph clustering approaches, as in Bruna et al. (2014) ; Defferrard et al. (2016) ; Simonovsky & Komodakis (2017) ; or learned through parameterization as in Ying et al. (2018b) ; Gao & Ji (2019) ; Lee et al. (2019) . In either approach, the result include both a hierarchy of graphs and the accompanying node representations.

Representation learning in these local pooling approaches is supervised, with the training signal coming from labels of the downstream task. In this work, we propose an unsupervised learning approach, named BAYESPOOL, through the use of variational Bayes. We use an existing coarsening Figure 1 : Coarsening sequence of a graph by using the method of Loukas (2019) .

method to obtain the graph hierarchy, an example of which is shown in Figure 1 . Then, our contribution is the learning of node representations for all graphs in this sequence. The high-level idea is to employ an encoder-decoder architecture: the encoder takes the graph and its node features as input and produces node features for the next graph in the hierarchy, whereas the decoder uses these produced node features to construct the graph. The objective is to obtain a decoding result as close to the given next graph as possible. The tool we use is variational Bayes. It is, however, slightly different from variational autoencoders (Kingma & Welling, 2014), because our decoder does not intend to reconstruct the input graph.

A clear benefit of unsupervised learning is that the learned representation is not tailored to a specific downstream task and hence may be more generalizable. Moreover, the coarsening method we adopt is a recent development that holds spectral guarantee (Loukas & Vandergheynst, 2018; Loukas, 2019) on the quality of the coarse graphs. We demonstrate the effectiveness of such a combination of hierarchical production and node representation learning in the context of graph classification. In particular, the classification performance is rather competitive with state-of-the-art supervised GNN approaches.

2. RELATED WORK:

This work is in part based on graph coarsening that produces a hierarchy for a given graph. Denote by G = (V, E) a graph with the vertex set V and the edge set E. Graph coarsening is concerned with computing a smaller (coarse) graph G c = (V c , E c ) with |V c | < |V | that retains the structure of G. A multilevel coarsening technique recursively coarsens the graph, yielding a hierarchy. Graph coarsening has been studied in the context of graph partitioning (Karypis & Kumar, 1998; Dhillon et al., 2007) , graph visualization (Harel & Koren, 2000) , machine learning (Lafon & Lee, 2006; Gavish et al., 2010; Ubaru & Saad, 2019) , and pooling in graph neural networks (Bruna et al., 2014; Defferrard et al., 2016; Simonovsky & Komodakis, 2017) . A variety of heuristic coarsening techniques have been proposed in different disciplines, including matching (Hendrickson & Leland, 1995; Ubaru & Saad, 2019) , first choice (Cong & Shinnerl, 2013) , contraction-based schemes (Dhillon et al., 2005; Sanders & Schulz, 2011) , and algebraic multigrid (AMG)-inspired schemes (Sharon et al., 2000; Hu & Scott, 2001; Ron et al., 2011; Chen & Safro, 2011) . Many well-known software packages exist for graph coarsening, e.g., Jostle (Walshaw & Cross, 2007) , Metis (Karypis & Kumar, 1998), and DiBaP (Meyerhenke et al., 2008) .

Recently, a few graph coarsening techniques achieving certain theoretical guarantees were presented (Moitra, 2011; Dorfler & Bullo, 2012; Loukas & Vandergheynst, 2018) . Loukas (2019) presented variational approaches for graph coarsening with spectral guarantees. In particular, it was shown that the coarse graphs preserve the top eigenspace (whose dimension is an input to the method) within a predefined error tolerance. Here, we use this variational approach to obtain the graph hierarchy. This work is concerned with unsupervised graph representation learning. Recent literature has focused on generative models to achieve the same; see. e.g., Kipf & Welling (2016) ; Li et al. (2018) ; Ma et al. (2018) ; Simonovsky & Komodakis (2018) ; Zhang et al. (2019) . For learning hierarchical representations of graphs, most of the works that we are aware of are based on supervised learning, including Bruna et al. (2014) ; Defferrard et al. (2016) ; Simonovsky & Komodakis (2017) ; Ying et al. (2018b) ; Gao & Ji (2019) ; Lee et al. (2019) . Methods most relevant to our work include: DIFFPOOL (Ying et al., 2018b) , where the coarsening matrices are learned in an end-to-end fashion; GRAPH U-NET (Gao & Ji, 2019) , where graph pooling is achieved using a learnable vector and node ranking; and SAGPOOL (Lee et al., 2019) , which is similar to GRAPH U-NET but uses graph self-attention to compute the ranking.

3. METHOD:

The proposed method BAYESPOOL is an extension of variational autoencoders. As the name suggests, the goal of an autoencoder is to reconstruct the original input object after encoding it in the latent space. Our approach does not reconstruct the original input, but rather, aims as decoding an output faithful to another prescribed object. To this end, we first revisit variational Bayes and justify the use of variational lower bound for learning. Then, the machinery is applied to the graph context.

3.1. VARIATIONAL BAYES:

Let x be the observed (data) variable and z be the unobserved (latent) variable. A core subject of Bayesian inference is concerned with estimating the posterior distribution p(z|x). It is related to the prior p(z) and the likelihood p(x|z) through the Bayes theorem

p(z|x) = p(x|z)p(z) p(x, z) dz .

The challenge lies in the marginalization over z in the denominator, which is generally computationally intractable. Hence, various approximations were developed. Typically one adopts a surrogate model q(z) independent of data; and recently in the context of VAEs, the data dependent distribution q(z|x) is often used. In our setting, we introduce a new variable x and consider q(z| x).

The difference, in terms of the Kullback-Leibler divergence, between the surrogate (variational) posterior q(z| x) and the true posterior p(z|x) may be decomposed as

D KL q(z| x) || p(z|x) = q(z| x) log q(z| x) p(z|x) dz = q(z| x) log q(z| x) p(z) dz DKL q(z| x) || p(z)

+ q(z| x) log p(x) dz

log p(x)

-q(z| x) log p(x|z) dz

E q(z| x) log p(x|z)

.

It consists of three terms: the KL divergence between the variational posterior and the prior p(z), the log-evidence log p(x), and the marginal log-likelihood log p(x|z) under the surrogate distribution.

Because the KL divergence is nonnegative, the log-evidence is lower bounded by the combination of the other two terms:

EQUATION

The better the surrogate, the tigher the lower bound.

One sees that the right-hand side of ( 1) is almost the same as the usual log-evidence lower bound (ELBO), except that the surrogate q(z| x) appears in place of q(z|x). This observation is not surprising, because the marginalization is over the latent variable z and has nothing to do with x and x. We thus conclude that the usual machinery of VAE applies, with only a notational change of the variational posterior. In the usual VAE setting, the first term of the right-hand side of (1) is considered the decoding accuracy, whereas the second term is a regularization in the latent space. Our setting follows this interpretation.

3.2. GRAPH REPRESENTATION LEARNING WITH VARIATIONAL BAYES:

In our setting, a pair of graphs-the original one and the coarse one-is given. Let A ∈ R n×n and A c ∈ R m×m be the corresponding graph adjacency matrices, respectively. Similarly, denote by X ∈ R n×d and X c ∈ R m×d the corresponding node feature matrices. We apply the encoderdecoder formalism, whereby the encoder encodes A and X into the coarse graph features X c that we seek, such that the decoder can use X c to decode a coarse graph as similar to A c as possible. See Figure 2 for an illustration.

A, X X c A c encoding decoding:

Figure 2 : Encoder-decoder architecture.

Specifically, in the language of generative modeling, the encoder is the parameterized inference model that produces the parameters of q(X c |A, X) and the decoder is the parameterized generative model that produces the parameters of p(A c |X c ). Following (1), A c plays the role of x, X c plays the role of z, and (A, X) plays the role of x. The variational lower bound for model learning is thus

EQUATION

Maxmizing the ELBO amounts to maximizing the likelihood of decoding the coarse A c (given coarse node features X c resulting from the encoder), while minimizing a regularization of the variational posterior q(X c |A, X) that departs from the latent distribution p(X c ).

3.3. MODELING AND PARAMETERIZATION:

Generally, the latent space may be kept simple and unparameterized, with more emphasis placed on the encoder and the decoder. Thus, we let the prior p(X c ) be the standard matrix normal MN (X c | 0 m×d , I m , I d ). Occasionally, specifying simple Gaussian structures on p(X c ) may improve performance, such as letting p(X c ) be the factored Gaussian, with the mean and the diagonal variance being parameters to learn. We have not yet, however, obtained strong empirical evidence of the benefit of using a parameterized factored Gaussian in this case.

For the decoder (generative model), a natural choice is to treat each element of the coarse graph adjacency matrix as a Bernoulli variable (scaled to the magnitude of the corresponding edge weight), with the success probability parameterized by a function of the corresponding coarse node features.

For notational convenience, let the column vector x c i ≡ X c (i, :) T and let A c ij ≡ A c (i, j). Then, p(A c |X c ) is the product of independent Bernoulli distributions:

p(A c |X c ) = i =j p(A c ij | x c i , x c j ) = i =j Bernoulli(1 A c ij | p ij ),

where 1 A c ij is the indicator function that returns 1 if A c ij = 0 and 0 otherwise, and p ij is a parameterized function that computes the success probability. A simple choice of the probability is an (unparameterized) dot product:

EQUATION

but there exist several other straightforward parameterized variants. For example,

EQUATION

where W and w are a parameter matrix and vector, respectively. One may also replace them by an MLP. Note that all these functions are symmetric with respect to i and j because the graph is undirected.

For the encoder (inference model), we treat the variational posterior q(X c |A, X) as a factored Gaussian: each coarse node x c i is an independent Gaussian with vector mean µ i and diagonal variance diag(σ 2 i ). Then,

q(X c |A, X) = i q(x c i |A, X) = i N (x c i | µ i , diag(σ 2 i )),

where µ i and σ i are parameterized functions of A and X. By doing so, the KL term in the ELBO (2) admits a closed form. We transpose the column vectors µ i and stack them to form a matrix M . Similarly, we proceed with the σ i 's and form a matrix S. In what follows, we model the parameterized function for M . The one for S is analogous.

Let C be the set of coarse nodes; hence A(C, :) keeps only the rows of A that correspond to the coarse nodes. We let M = σ(HXW 1 ), (5) where W 1 is a parameter matrix, H has the same nonzero pattern as A(C, :), and σ is an activation function. This expression is in form similar to one graph convolution layer in GCN (Kipf & Welling, 2017) , except that the square normalized adjacency matrix is replaced by a fat rectangular matrix H. Rather than basing H on the original adjacency matrix, we designate its nonzero elements to be self-attention weights computed from the node features. Specifically, the nonzero elements of the i-th row of H is computed as

EQUATION

where w 2 is a parameter vector and W 3 and W 4 are parameter matrices. As an alternative, one may replace the attention calculation by that in GAT (Velicković et al., 2018) :

EQUATION

)

To further enhance the representational power, in the parameterization (5) one may replace the original feature matrix X by the node embedding matrix Z output from GCN: M = σ(HZW 1 ) with Z = GCN(A, X; W 5 ).

(8) The GCN introduces additional parameters W 5 that may be useful for large data sets.

In Section 4, we experiment with the different variants and suggest a default choice that works generally well.

3.4. MULTILEVEL LEARNING:

Coarsening may be done recursively, forming a sequence of increasingly coarse graphs. Let the adjacency matrices of this sequence be A 0 , A 1 , . . . , A L , where A 0 = A corresponds to the initial given graph. Given this sequence and the initial feature matrix X 0 = X, the goal is to obtain the subsequent feature matrices X 1 , . . . , X L .

To this end, model learning is conducted through maximizing the evidence of observing A 1 , . . . , A L , treated independently. That is, we want to optimize log p(A 1 ) + • • • + log p(A L ). Following the argument as before, the actual quantity to optimize is the evidence lower bound. Inserting the layer index into (2), we have

ELBO = E q(X |A -1 ,X -1 ) log p(A |X ) -D KL q(X |A -1 , X -1 ) || p(X ) .

Then, the log-evidence lower bound is the sum

ELBO = L =1 ELBO .

The encoder and decoder parameters across coarsening levels differ but they are jointly learned through maximizing the ELBO.

3.5. DOWNSTREAM TASKS:

The learned features X 1 , . . . , X L , together with the original X 0 , may be used for predictive tasks through learning a separate predictive model. We follow a common practice and define the model as

y p = MLP concat readout(X 0 ), readout(X 1 ), • • • , readout(X L ) ,

where readout is a global pooling across graph nodes (e.g., a concatenation of the max pooling and the mean pooling), concat denotes vector concatenation, MLP is self-explanatory, and y p is the class probability vector. In this paper, we consider the graph classification task.

4. EXPERIMENTS:

In this section, we evaluate the performance of BAYESPOOL through the task of graph classification. Note again that BAYESPOOL is an unsupervised method; but as we will see, it is rather competitive with recently proposed supervised methods, even outperforming them on several data sets. We first present the details of the experimented data sets, the compared methods, and the training procedure. Then, we compare the graph classification accuracies. We also perform sensitivity analysis regarding the number of coarsening levels and compare the performance of several variants in the implementation of BAYESPOOL.

Data sets: We consider the same data sets used by Lee et al. (2019) . They are standard benchmarks publicly available from Kersting et al. (2016) . Table 1 summarizes the information. (Wale et al., 2008) . Here, the graphs correspond to chemical compounds, with the atoms and the chemical bonds represented as the nodes and edges, respectively. The FRANKEN-STEIN (Orsini et al., 2015) data set contains molecular graphs with 780 node features. The labels indicate if a molecule is mutagen or non-mutagen. Data sets DD, NCI1, and NCI109 do not come with node attributes for use as features. Hence, we employ transformations of node degrees as features, following the practice of DIFFPOOL.

Compared methods: We compare with three supervised hierarchical graph representation learning methods, namely DIFFPOOL, GRAPH U-NET, and SAGPOOL. We duplicate the test accuracies of these methods reported by Lee et al. (2019) .

DIFFPOOL (Ying et al., 2018b) computes hierarchical representations of graphs by using end-to-end trainable pooling based on soft clustering. This method is expensive because of the computation of the dense projection matrix. Lee et al. (2019) report that the method ran out of memory for a pooling ratio greater than 0.5. GRAPH U-NET (Gao & Ji, 2019 ) uses a learnable scoring vector to rank the graph nodes and selects top-ranked nodes for pooling. SAGPOOL (Lee et al., 2019) uses a similar approach as GRAPH U-NET for pooling, but incorporates a self-attention layer for learning the scoring vector. In both methods, the pooling ratio was set to be 0.5.

Training procedure: We follow Lee et al. (2019) and perform several rounds of random splits with 80% training, 10% validation and 10% test. The learning rate is tuned over the range [1e-2, 5e-2, 1e-3, 5e-3, 1e-4]. The hidden dimensions are tuned over [10, 20, 32, 48] . For LeakyReLU, the slope coefficient is set to 0.01. The coarsening ratio ρ is experimented with [0.25, 0.5, 0.75]; see Section 4.2 for sensitivity analysis. The dimension of the top eigenspace to be preserved by the coarsening procedure is set to K = 5. We implement the method in PyTorch and use Adam as the optimizer. The training employs early stopping with a patience of 50 epochs. The training setting for other methods are reported in Lee et al. (2019) .

For a fair comparison, we use the same MLP classifier as in Lee et al. (2019) . For readout, we use mean pooling and max pooling and concatenate the two outputs. We then use 3 feedforward layers along with softmax for classification. The classifier is trained for 150 epochs.

For architecture variation, the following combination consistently achieves the best results: the decoder uses an unparameterized dot product (3), the encoder uses parameterization (5) with the attention matrix H computed by using the GAT form (7). The classification results reported in Section 4.1 below follow this choice. The results of other variants are reported in the subsections that follow.

The code is available at https://anonymous.4open.science/r/ a50d6411-55f7-4e24-8f6c-6eecee118ea0/.

4.1. GRAPH CLASSIFICATION:

We compare the performance of BAYESPOOL with several high-performing supervised methods recently proposed. These methods are all hierarchical methods. Table 2 lists the average test accuracies. The pooling/coarsening ratio is 0.5 in all cases. BAYESPOOL outperforms other methods on three out of five data sets. Its performance is also on par with DIFFPOOL on the other two data sets, although infereior to GRAPH U-NET and SAGPOOL. The graphs in DD and PROTEINS are relatively large; with a coarsening ratio 50% there does not seem to cause information loss. Hence, BAYESPOOL works rather appealingly. On the other hand, although the graphs are smaller in FRANKENSTEIN, the data set contains a large number (780) of features, which possibly dwarf the graph structure information. Therefore, all methods yield similar results (with ours slightly outperforming the others). Encouragingly, BAYESPOOL is an unsupervised method; hence, these results show that the method is highly competitive for downstream tasks such as graph classification. It enables incorporating sophisticated coarsening techniques such as Loukas (2019) for graph representation learning.

Note that the adjacency matrices of the graphs are typically sparse. BAYESPOOL leverages sparse matrix computation, as opposed to DIFFPOOL where the projection matrix does not have an a priori sparsity structure. The coarsening procedure used by BAYESPOOL is implemented in sparse matrix format, along with the calculations of the neural network. This implementation results in a lower time cost and space complexity.

4.2. EFFECT OF COARSENING RATIO:

One of the key factors that affects the performance of BAYESPOOL is the amount of graph reduction (the ratio of the coarse graph size to the initial size), or equivalently, the number of coarsening levels. This ratio is an input parameter to the coarsening method of Loukas (2019) that we use. In Table 3 , we evaluate the performance of BAYESPOOL on two data sets with respect to the coarsening ratio. These data sets have relatively large graphs so that aggressive coarsening is possible. In the table, we report the results for three different levels of coarsening with ratio ρ = m/n, where m is the number nodes in the coarse graph and n the original graph. We observe that the performance of BAYESPOOL is relatively stable when ρ ≥ 0.5, but degrades as ρ becomes smaller. As we lower ρ, more and more nodes and edges are removed, causing significant loss of information. However,

:



Figures and Tables:

Data set statistics. The first two data sets are related to protein structures. The graphs in DD (Dobson & Doig, 2003) have different amino acids as nodes; and the edges correspond to the distance between the nodes. Labels indicate if the protein is an enzyme or not. The PROTEINS (Dobson & Doig, 2003) graphs have secondary structure elements of proteins as nodes. The edges indicate whether the nodes are in amino acids. NCI1 and NCI109 are biological data sets popularly used for anticancer activity classification

Table content: <table><tr><td>Data set</td><td colspan="5">No. graphs Nodes (max) Nodes (avg) Edges (avg) Classes</td></tr><tr><td>DD</td><td>1178</td><td>5748</td><td>284.32</td><td>715.66</td><td>2</td></tr><tr><td>PROTEINS</td><td>1113</td><td>620</td><td>39.06</td><td>72.82</td><td>2</td></tr><tr><td>NCI1</td><td>4110</td><td>111</td><td>29.87</td><td>32.30</td><td>2</td></tr><tr><td>NCI109</td><td>4127</td><td>111</td><td>29.68</td><td>32.13</td><td>2</td></tr><tr><td>FRANKENSTEIN</td><td>4337</td><td>146</td><td>16.90</td><td>17.88</td><td>2</td></tr></table>

Graph classification accuracies.

Table content: <table><tr><td>Data set</td></tr></table>

Graph classification accuracies for different coarsening ratios.

Table content: <table><tr><td>Coarsening ratio</td><td>DD</td><td>PROTEINS</td></tr><tr><td>ρ = 0.75</td><td>77.87</td><td>77.27</td></tr><tr><td>ρ = 0.5</td><td>78.12</td><td>76.85</td></tr><tr><td>ρ = 0.25</td><td>68.75</td><td>66.96</td></tr></table>

Peer review (iclr2020):
The authors propose in this paper a new unsupervised graph representation learning method. The method leverages recent advances in graph coarsening, mainly Loukas* method. The key idea of the method consists in using a reconstruction target that is not the classical one in an auto-encoder setting. More precisely, the encoder takes as an input the original adjacency matrix and node features but the decode only aims at reconstructing the coarse adjacency matrix (obtained via Loukas* method). The experimental evaluation is quite thorough and shows that the method performs quite well, especially considering it is unsupervised but is compared to supervised representation methods. It would be nice to include statistical tests to assess the significance of the differences in cases were accuracies are very close one to another. A missing part would be to explore the relevance of the learned representation for other tasks (i.e. to use a multi task data set). Of course as the representation is learned in an unsupervised way, one can argue that the current evaluation is already providing an answer. Overall, I find the paper clear, but the variational bayes part could be much clearer. In fact I*m not sure why this is presented as variational bayes and not only variational. I do not see any prior distribution over parameters, for instance. I understand that the recent *tradition* in variational auto-encoder is to use this terminology, but as a (part time) bayesian, this is a bit annoying.