Year: 2021

Abstract:
Graph representation learning is a fundamental problem for modeling relational data and benefits a number of downstream applications. Traditional Bayesian-based random graph models and recent deep learning based methods are complementary to each other in interpretability and scalability. To take the advantages of both models, some combined methods have been proposed. However, existing models are mainly designed for undirected graphs, while a large portion of real-world graphs are directed. The focus of this paper is on the more challenging directed graphs where both the existences and directions of edges need to be learned. We propose a Deep Latent Space Model (DLSM) for directed graphs to incorporate the traditional latent space random graph model into deep learning frameworks via a hierarchical variational auto-encoder architecture. To adapt to directed graphs, our model generates multiple highly interpretable latent variables as node representations, and the interpretability of representing node influences is theoretically proved. The experimental results on real-world graphs demonstrate that our proposed model achieves the state-of-the-art performances on link prediction and community detection tasks while generating interpretable node representations of directed graphs.

heading: 1 INTRODUCTION
text: Learning representations on graphs is a fundamental problem for graph analysis and benefits a number of downstream applications, typically including link prediction, community detection and influential node identification. Traditionally, a plethora of Bayesian-based random graph models have been proposed for learning graph representations (Holland et al., 1983; Hoff et al., 2002; Karrer & Newman, 2011; Sewell & Chen, 2015). Despite the ideal interpretability of these methods, they are unscalable to model large-scale graphs due to the expensive iterative inference procedures. Boosted by the powerful representation learning ability, graph neural networks (GNNs) have been proposed to learn the topology of graph-structured data (Kipf & Welling, 2017; Hamilton et al., 2017; Veličković et al., 2018; Zeng et al., 2020). However, these GNN-based methods usually bring about the interpretability issues, which restrict the practical applications of these methods.
As an effective deep learning method based on the Bayesian theory, the variational auto-encoders (VAEs) (Kingma & Welling, 2014b) have been proposed for representation learning on graphstructured data (Kipf & Welling, 2016; Grover et al., 2019; Salha et al., 2019). To capture the mutual facilitation of random graph models and deep learning based methods, some recent research attempts to combine them using VAEs, and build the deep generative random graph models (Mehta et al., 2019; Sarkar et al., 2020). However, most of these VAE-based methods are designed only for undirected graphs. In contrast, a large portion of real-world graphs, such as social and citation networks, are directed graphs, of which undirected graphs are special cases (i.e. an undirected edge can be regarded as a bidirectional edge). The representation learning on directed graphs is particularly challenging since not only the existences but also the directions of edges need to be learned. Although directed graphs have been extensively studied by traditional research on random graph models (Hoff et al., 2002; Krivitsky et al., 2009; Sewell & Chen, 2015), there are still few deep learning based methods and no combined models existed to address them.
In this paper, we present the first VAE-based generative model for directed graphs and propose the Deep Latent Space Model (DLSM), which combines graph convolutional networks (GCNs) with the
Bayesian random graph model using a hierarchical VAE architecture. Within our model, a deep GCN is leveraged to encode the asymmetric adjacency matrix and node features as hidden states, which are layer-wise transformed to variational parameters to generate latent random variables as interpretable node representations. Then, the hierarchical decoder network reconstructs the adjacency matrix using the learned representations.
To better capture the characteristics of graphs for deep random graph models, our method generates three types of interpretable representations, including the latent positions, node random factors and community memberships. The latent positions measure the distances between nodes and enable our method to model the link reciprocity of directed graphs, which means that both bidirectional and unidirectional edges exist. To characterize the degree heterogeneity, we introduce a pair of node random factors, namely the social activity and popularity, thus the asymmetric structure of directed edges can be modeled. Besides, we use a binary latent variable, named the community membership, to accommodate an overlapping community structure. We further prove that these node-specific variables, including the latent positions and random factors, are interpretable for representing the existence and strength of node influences.
For fast inference, we leveraged the stochastic gradient variational Bayes (SGVB) algorithm (Kingma & Welling, 2014a), which is far more efficient and scalable compared to the traditional iterative inference methods. We conduct experiments on link prediction task to verify the effectiveness of directed graph modeling with several real-world graphs, and further demonstrate the interpretability of latent variables on community detection task.
The main contributions of this work are summarized as follows:
1. We propose the DLSM model for directed graphs to incorporate the latent space model into deep learning frameworks using a hierarchical VAE architecture, which can generate highly interpretable representations of graph nodes for multiple downstream tasks.
2. To adapt our method to directed graphs, we introduce the social activity and popularity factors of each node, with theoretical proof of their capabilities in representing node influences.
3. We conduct experiments on five real-world graph datasets, and the results show that our proposed method can well support downstream tasks and achieve the state-of-the-art performances in link prediction and community detection.

heading: 2 RELATED WORK
text: In this section, we briefly review the traditional Bayesian-based random graph models and deep learning based generative graph representation methods.

heading: 2.1 BAYESIAN-BASED RANDOM GRAPH MODELS
text: Classic Bayesian-based random graph models have developed for decades and are still valued in modeling and generating relational data. This methods regard graph nodes as random variables following a prior distribution, and learns low-dimensional embeddings via the posterior distributions of nodes. One of the most well-known methods is the stochastic blockmodel (SBM) (Holland et al., 1983), which generates a latent variable indicating the community membership of each node. Following SBM, a large number of variants have been proposed, such as allowing an overlapping community structure (Miller et al., 2009) and involving degree heterogeneity (Karrer & Newman, 2011).
Another important random graph model is the latent space model (LSM) (Hoff et al., 2002), which assumes each node as a position in an unobservable latent space and use distances to measure the relationships between nodes. Such method has soon been extended to model directed graphs by involving the degree heterogeneity of nodes (Krivitsky et al., 2009; Sewell & Chen, 2015). Although these above methods provide good theoretical properties and learn interpretable node embeddings, they rely on either the MCMC posterior sampling or mean-field variational inference with dramatically high computational complexity, and thus are usually unfeasible when modeling large-scale graphs.

heading: 2.2 DEEP LEARNING BASED MODELS FOR UNDIRECTED GRAPHS
text: The intriguing achievements of deep learning models for learning representations of Euclidean data have encouraged efforts to employ them on graph-structured data. The earliest attempts include DeepWalk (Perozzi et al., 2014) and node2vec (Grover & Leskovec, 2016), both of which encode local relations as node representations by conducting random walks on a graph.
GNN-based graph models More recently, many GNN-based models have been proposed for graph representation learning. Kipf & Welling (2017) first leverages the spectral GCN to learn node embeddings using the global topology of a graph. GraphSAGE (Hamilton et al., 2017) samples a fixed-size neighborhood of each node for the inductive representation learning, which allows unseen nodes to be excluded during training. The graph attention networks (Veličković et al., 2018) introduce the attention mechanism to aggregate the information from neighbors. GraphSAINT (Zeng et al., 2020) constructs minibatches by sampling subgraphs across GCN layers. (Zhang et al., 2021a) proposes a multi-node representation learning method for link prediction using the node labeling trick. These methods usually learn representations for a single or a set of nodes and only focus on some specific task of graph modeling, such as node classification or link prediction.
Deep generative graph models Another stream of research is the deep generative models, which attempt to generate the full structures of graphs and are typically able to be applied for multiple downstream tasks. These methods usually consist of a inferential model to capture the inner probabilistic distribution of graphs and an unsupervised generative model to generate the similar data. The variational graph auto-encoders (VGAE) (Kipf & Welling, 2016), for example, sample node embeddings from variational Normal distributions to reconstruct the adjacency matrix. GraphGAN Wang et al. (2018) unifies a generator and a discriminator to learn node embeddings by playing a minimax game. GraphRNN You et al. (2018) generates small graphs through a sequence of node and edge formations. Recently, some work tries to incorporate traditional random graph models into deep learning frameworks. Mehta et al. (2019) first combines the classic SBM with GCN to learn sparse node embeddings. Furthermore, Sarkar et al. (2020) builds a Gamma ladder VAE (Sønderby et al., 2016) architecture to discover the communities at multiple levels of granularities.

heading: 2.3 DEEP LEARNING BASED MODELS FOR DIRECTED GRAPHS
text: Despite a plethora of literature for undirected graph modeling, the deep learning based methods for directed graphs are still relatively rare. Existing methods are mostly designed for some special cases and are not generalizable enough. For instance, Salha et al. (2019) and Zhang et al. (2021b) focus on the unidirectional graphs and neglect the link reciprocity of generalized directed graphs. Funke et al. (2020) learns the low-dimensional statistical manifold embedding for unsupervised directed graph modeling. Zhu et al. (2021) learns a pair of source and target embeddings for each node using the generative adversarial network. In this paper, we propose a novel deep generative model and learn interpretable node representations for generalized directed graphs where both uni- and bi-directional edges exist.

heading: 3 DEEP LATENT SPACE MODEL
text: Consider a directed graph containing n nodes. The input data include an asymmetric adjacency matrix A = (Aij) ∈ {0, 1}n×n, where each binary element Aij denotes the presence (1) or not (0) of the directed edge from node i to j, and a node attribute matrix X ∈ Rn×p. Throughout this paper we assume all edges to be conditionally independent and satisfy Aij |Θ ∼ Bernoulli(pij), where Θ is the collection of latent variables (representations) and pij ∈ (0, 1) is the posterior probability to form an edge. The objective is to learn the representations which can best reconstruct the adjacency matrix.
We propose a deep VAE architecture composed of a GCN encoder and a hierarchical latent space model (HLSM) decoder, as shown in Fig. 1. The encoder block takes the adjacency A and attribute matrices X, if available, as inputs and learns a hidden state for each node using a directed GCN. Then, the decoder block recursively generates interpretable latent variables from variational posteriors via Monte Carlo sampling and reconstructs the adjacency matrix as outputs. For inference, the
hierarchical architecture enables the model to learn variables using not only the priors passed from previous layers of the decoder, but also the approximate likelihood learned by the encoder.

heading: 3.1 HLSM DECODER
text: The HLSM decoder learns interpretable node representations based on the classic LSM approach and generates an adjacency matrix. Assume each node to correspond to a latent position in an unobservable D-dimensional space, denoted as zi ∈ RD. The distances between latent positions indicate the relationships of nodes. The closer two latent positions are, the more possible an edge will exist. Formally, an directed edge is generated by following
P (Aij = 1|X,Θ) = σ (β0 − βout‖γi (zi − zj) ‖ − βin‖δj (zi − zj) ‖) , (1)
where σ(·) is the sigmoid function. The latent positions zi are involved as Euclidean distances, indicating the relationships between nodes. γi and δi ∈ RD are the activity and popularity factors, the reverse of which represent the tendencies for a node to send and receive an edge, respectively. The global weights βout and βin are to measure the importance of activity and popularity, respectively, and β0 is a bias. To be mentioned, though Eq. (1) is designed for directed graphs, it can be easily degraded to an undirected variant by setting γi = δi and βout = βin.
We build an HLSM decoder to generate three types of latent variables, i.e. the latent positions z (l) i ∈ RGl , the community membership s (l) i ∈ RGl and the node random factors γi, δ (l) i ∈ RGl , where Gl is the size of the l−th decoder layer. All of these variables are randomly sampled from the variational distributions. At the first layer (l = 1), the parameters of variational distributions are defined by priors solely, while for other layers (l = 2, . . . , L− 1), the parameters are obtained by the feedforward representations generated from previous layers as well as priors.
Latent position The latent positions z(l)i represent the location of each node in a latent space, generated as
z (l) i ∼ Normal ( s (l) i f(W (l−1) z z (l−1) i ), diag(σ (l) i 2 ) ) , (2)
where f(·) is a nonlinear activation function (e.g. leaky ReLU), denotes the element-wise multiplication, diag(·) denotes a diagonal matrix, σ(l)i 2 ∈ RGl is the prior variance to be specified, and W(l−1)z ∈ RGl×Gl−1 is a weight matrix to transform the variable from dimension Gl−1 to Gl.
Community membership The latent positions in Eq. (2) are separated into different communities by the sparse binary variables s(l)i = (si1, . . . , siGl) ′, thus an overlapping community structure can
be adapted. Referring to Miller et al. (2009), we employ the Indian buffet process (IBP) prior on s(l)i to learn the effective number of communities given a sufficiently large truncation parameter Gl, i.e.
π (l) ig = logit
( g∏
h=1
v (l) ih
) , s
(l) ig ∼ Bernoulli
( σ(π (l) ig ) ) . (3)
Typically, the log odds π(l)ig is generated using a stick-breaking construction, where v (l) ih is drawn from a Beta distribution (Teh et al., 2007). In our model we simplify such hierarchical prior structure by specifying a global v for all nodes. At each layer, the community membership s(l)ig denotes whether node i belongs to community g, hence the size of each decoder layer Gl can be explained as the number of communities. Additionally, the proposed HVAE architecture enables our model to detect community structures at multiple levels of granularities. Letting the layer sizes Gl be downward increasing, the top layers indicate the coarse-grained communities and the bottom layers indicate the fine-grained communities.
Node random factors To model the prevalent power-law of node degrees, we propose the Dirichlet random factors γ(l)i , δ (l) i ∈ RGl , i.e.
γ (l) i ∼ Dirichlet ( ξ (l) i + s (l) i f(W (l−1) γ γ (l−1) i ) ) , (4)
δ (l) i ∼ Dirichlet
( ψ
(l) i + s (l) i f(W (l−1) δ δ (l−1) i )
) , (5)
where ξ(l)i and ψ (l) i are prior parameters to be specified, and W (l−1) γ , W (l−1) δ ∈ RGl×Gl−1 are weight matrices. Note that the node random factors are dependent on the community membership s(l)i as well, which can be explained as only the random effects within the community are involved in the representations.
At the output layer (l = L), the latent positions zi and node random factors γi, δi are obtained by linear transformation from dimension GL−1 to D. Each column of the weight matrices belongs to a (D − 1)-simplex. Such transformation also changes the interpretation of the layer size, from the number of communities to the dimension of the latent space. Last, the adjacency matrix is reconstructed with the latent variables by Eq. (1).

heading: 3.2 GCN ENCODER
text: The proposed DLSM employs a deep encoder as a non-iterative recognition model to infer the parameters of posterior distributions. Assuming the mean-field approximation of the variational distributions, the true joint posterior of the latent variables pθ(Θ|A,X) can be approximated by a variational posterior qφ(Θ), where θ and φ denote the generative (decoder) and inference (encoder) parameters to be trained, respectively. Then, the variational posterior is given as
qφ(z, s,γ, δ) = n∏ i=1 L∏ l=1 qφ(z (l) i |h (l) i , z (l−1) i )qφ(si|h (l) i ,π (l−1) i )
qφ(γi|h(l)i ,γ (l−1) i )qφ(δi|h (l) i , δ (l−1) i ), (6)
where h(l)i ∈ RKl is the output of the encoder and Kl denotes the size of the l-th layer. GCN has been proved effective in learning the topology of non-Euclidean data, and thus is an ideal choice for the encoder of our model. Referring to Kipf & Welling (2017), we propose the directed GCN operator as
H(l+1) = f ( D̃outÃD̃inH(l)W(l) ) . (7)
Here Ã = A + In and In is the n-dimensional identity matrix. D̃out and D̃in are diagonal matrices with elements as the out- and in-degrees of Ã, respectively. H(0) = X if X is available and H(0) = In if not. During inference, the hidden states are passed to the corresponding layer of the decoder and then combined with the prior information from previous layers to generate the parameters of variational distributions. Note that the vanilla GCN used here can be substituted by any other GNN for directed graphs.

heading: 4 INFERENCE
text: We now introduce our fast inference method using the SGVB algorithm (Kingma & Welling, 2014a). Compared with the iterative methods such as MCMC adopted by traditional Bayesian random graph approaches, SGVB is much more efficient and scalable. Such method requires differential Monte Carlo expectations to perform backpropagation, thus the reparameterization trick for each of the latent variables is leveraged.
Denoting z∗ ∈ RGl a vector with standard Normal elements, the latent positions are reparametrized as z (l) i = µ̂ (l) i (z (l−1) i , s (l−1) i ,h (l) i ) + σ̂ (l) i (z (l−1) i , s (l−1) i ,h (l) i ) z∗, where µ̂ (l) i (z (l−1) i , s (l−1) i ,h (l) i ), σ̂ (l) i (z (l−1) i , s (l−1) i ,h (l) i ) are variational posterior parameters.
Following Mehta et al. (2019), the Bernoulli posterior pθ(s (l) ig ) is approximated by the Binary Concrete distribution (Maddison et al., 2017), i.e.
s̃ (l) ig =
1
λ
( π̂ (l) ig (s (l−1) i ,h (l) i ) + logit(u) ) , (8)
where u ∼ U(0, 1), π̂(l)ig (s (l−1) i ,h (l) i ) is a variational posterior parameter, λ ∈ R+ is the temperature to be specified, and then s(l)ig = σ(s̃ (l) ig ).
Referring to the recent literature (Joo et al., 2020), we use multiple normalized Gamma variables with unified rate parameters to compose the Dirichlet distributions pθ(γ (l) i ) and pθ(δ (l) i ), i.e.
γ̃ (l) i ∼ Gamma ( ξ̂i (l) (γ (l−1) i , s (l−1) i ,h (l) i ),1 ) , (9)
δ̃ (l) i ∼ Gamma ( ψ̂i (l) (δ (l−1) i , s (l−1) i ,h (l) i ),1 ) , (10)
where ξ̂i (l) (γ (l−1) i , s (l−1) i ,h (l) i ) and ψ̂i (l) (δ (l−1) i , s (l−1) i ,h (l) i ) are variational posterior parameters. For the convenience of notations, here the parameters of Gamma distribution are symbolized as vectors, meaning the element-wise operations. The node random factors are then derived by γ(l)i = γ̃ (l) i / ∑n j γ̃ (l) j , δ (l) i = δ̃ (l) i / ∑n j δ̃ (l) j . In practice, the Dirichlet variables are magnified by n times to avoid too small values of γ(l)i and δ (l) i when n is large.
In particular, the initial variational posterior parameters (l = 0) are set as nonlinear combinations of h (1) i and h (L−1) i , as illustrated in Fig. 1.
The loss function is defined by minimizing the negative evidence lower bound (ELBO), i.e.
L = n∑ i=1 L∑ l=1 ( KL [ qφ(z (l) i ) ∥∥∥pθ(z(l)i )]+ KL [qφ(s(l)i )∥∥∥pθ(s(l)i )]+ KL [qφ(γ(l)i )∥∥∥pθ(γ(l)i )]
+KL [ qφ(δ (l) i ) ∥∥∥pθ(δ(l)i )])− n∑
i=1 n∑ j=1 Eq [ logpθ(Aij |Θ(L)) ] , (11)
where KL[q(·)||p(·)] denotes the Kullback-Leibler (KL) divergence between q(·) and p(·) (of which the full expression is given in Appendix A). The second term is the cross entropy of adjacency matrix reconstruction. Here all of the true posteriors pθ(Θ(l)) and variational posteriors qφ(Θ(l)), except for the input layer (l = 1), are conditioned on Θ(l−1), which is omitted for simplification.

heading: 5 MODEL INTERPRETABILITY
text: In this section, we show that the proposed latent positions and node random factors are interpretable for representing the existence and strength of node influences, by hypothesizing that the true posterior of zi can be approximated by the variational Normal distribution.
Imagine that there is a new node i with ni ≥ 2 neighbors entering a graph and we want to calculate zi using other existed nodes. It is rational to suppose that the new node should be posited around its
neighbors and thus can be denoted as a residual to the mean of its neighbors, i.e.
zi ∼ Normal (z̄i + i, diag (σi)) , (12) where z̄i = 1ni ∑ j∈Ni zj is the mean position of node i’s neighbors, and i ∈ R
D is the residual for zi to the mean position. Intuitively, the residual should depend on all of the neighbors to various extent, of which the influential ones are likely to contribute more. Formally, we decompose the residual as
i = ∑ j∈Ni rijeij , (13)
where rij ≥ 0 is a scalar and eij ∈ RD is a unit vector at the direction of zj − z̄i. The non-negative scalar denotes the strength of influences for neighbor j exerting on node i (see Fig. 2 in Appendix for an example in a two-dimensional latent space). An extreme case is rij = 0, which means the latent position of i is not influenced by neighbor j at all.
To further analyze the relation between the latent variables and node influences, we define the prior of rij as a mixture of a point mass on 0 and an exponential distribution, i.e.
fr(rij) =
{ p0 if rij = 0
(1− p0)fExpλ(rij) if rij > 0, (14)
where 0 < p0 < 1 and fExpλ(·) denotes the density function of the exponential distribution with parameter λ. What is of our interest is the posterior odds ratio (OR) for rij > 0, defined as OR(rij > 0|A,X) = ∫∞ 0 fr(rij = v|A,X)/fr(rij = 0|A,X)dv, which measures the ratio for node influences to exist.
Theorem 1 Let Θ denote a collection of all parameters except zi. The posterior OR for rij > 0 satisfies
OR(rij > 0|A,X) = E [ g (zi,Θ) ∣∣∣rij = 0,A,X] , (15) The full expression and proof are provided in Appendix B.1. It is interesting that in the proof (zi − z̄i − ∑ k∈Ni,k 6=j rikeik)
′eij is a scalar projection of zi’s residual component onto the unit vector whose direction is determined by zj . In a word, it shows that the posterior latent positions are able to model the existence of node influences via the residual component projections.
Now we consider the interpretation of node random factors, which explains the strengths of node influences. Formally, we have the following conclusion, i.e.
Theorem 2 Assume the vectors from z̄i to all influential neighbors (rij > 0) are orthogonal to each other, and the random factors are the same at all dimensions of the latent space, denoted as γi and δi, respectively. The strength of influences for neighbor j exerting on node i can be denoted as
rij =h(γi, δj , θ, o)‖zj − zi‖, (16)
where θ is the angle between zi − z̄i and zj − z̄i, o is the angle between z̄i − zj and zi − zj .
Here h(γi, δj , θ, o) is a function of the random factors as well as θ and o. The full expression and proof are given in Appendix B.2. Theorem 2 reflects that given the direction of eij and zj − zi, the node random factors can be regarded as coefficients for measuring the ratio of rij to ‖zj − zi‖, which means that the strength of influences for a neighbor delivering to the node can be strengthened or undermined via shrinking or enlarging the distance between them by the random factors.

heading: 6 EXPERIMENTS
text: To evaluate the effectiveness and interpretability of the proposed model, we conduct a series of experiments, including link prediction and community detection, on several real-world graphs. Ablation study for the components of our method, including the community membership, node random factors and the HLSM decoder architecture are also conducted.

heading: 6.1 DATASETS
text: The experiments of link prediction are conducted on five real-world benchmark datasets, namely Emails (Leskovec et al., 2007), Political blogs (Adamic & Glance, 2005), Cora (Sen et al., 2008), WikiVote (Leskovec et al., 2010) and Google (Palla et al., 2007). All of the graphs have been preprocessed by omitting the isolated nodes and loops. In our experiments, the edges are randomly splitted as 85% for training, 10% for testing, and 5% for validation. More details and descriptive statistics of the datasets are given in Appendix C.

heading: 6.2 BASELINES
text: We compare our proposed DLSM1 with four recent deep generative methods for graph representation learning, i.e. the variational graph auto-encoder (VGAE) (Kipf & Welling, 2016), the deep generative latent feature relational model (DGLFRM) (Mehta et al., 2019), the gravity graph variational autoencoder (GGVAE) (Salha et al., 2019), the ladder Gamma variational auto-encoder for graphs (LGVG) (Sarkar et al., 2020), and the directed graph generative adversarial network (DGGAN) (Zhu et al., 2021). For link prediction, we also consider a non-generative GNN-based model, i.e. SEAL (Zhang et al., 2021a), which is the state-of-the-art method for this task. The methods that are not particularly designed for directed graphs, including VGAE, DGLFRM and LGVG, are slightly modified by altering the original encoder with our proposed directed GCN given in Eq (7).

heading: 6.3 EXPERIMENTAL RESULTS
text: Link Prediction For link prediction, we employ the widely used area under the ROC curve (AUC) and average precision (AP) as evaluation metrics. Hyperparameter settings are provided in Appendix D. The experimental results of DLSM and the baselines are presented in Table 1, where the reported results are the means and standard deviations of 10 independent random splits. On all datasets, our model achieves better or comparable performances to the non-generative method SEAL. However, since it requires to materialize a subgraph for each link, SEAL is much more unscalable than our DLSM and takes about 20 multiples longer for running one epoch. The traditional LSM model is unpractical to fit large graphs such as WikiVote and Google because of the MCMC-based inference methods, while the SGVB method adopted by our model is much more computationally efficient. In addition, it seems that GGVAE performs well on approximately unidirectional graphs (i.e. the reciprocal rate is close to 0 as shown in Table 4) such as Cora, whereas on other more generalized directed graphs where more bidirectional edges exist, our method is significantly superior to GGVAE. This is because our proposed DLSM considers the link reciprocity of directed graphs, which is neglected by GGVAE due to its absolutely asymmetric decoding scheme. We also conduct experiments of link prediction on three real-world undirected graphs, given in Appendix E, which show that our proposed model achieves comparable performances with the baselines and verify the capability for modeling undirected graphs as special cases.
Community Detection We further conduct community detection on three real-world datasets with ground-truth community labels, i.e. Emails, political blogs and Cora. In particular, for the Emails
1The source code is available at https://github.com/upperr/DLSM.
VGAE 31.5 ± 1.6 42.5 ± 1.5 34.9 ± 1.2 52.3 ± 1.4 23.2 ± 1.0 24.8 ± 0.4 DGLFRM 53.3 ± 1.3 60.7 ± 1.0 42.4 ± 1.1 50.5 ± 1.3 18.2 ± 0.5 22.2 ± 0.5 GGVAE 62.8 ± 1.5 67.7 ± 1.4 92.0 ± 1.1 92.5 ± 0.8 22.7 ± 0.7 24.0 ± 0.4 LGVG 67.0 ± 1.3 63.8 ± 1.2 72.8 ± 1.3 73.9 ± 1.3 20.0 ± 0.5 23.9 ± 0.7 DGGAN 62.1 ± 1.6 66.7 ± 1.3 88.9 ± 1.1 89.5 ± 1.3 18.7 ± 1.0 21.4 ± 1.2
DLSM 72.0 ± 1.0 77.0 ± 1.2 92.5 ± 0.7 92.6 ± 0.7 24.6 ± 0.8 25.7 ± 0.5
DLSM 95.2 ± 0.2 94.1 ± 0.6 94.8 ± 0.4 93.7 ± 0.8 91.6 ± 0.2 92.5 ± 0.2 97.0 ± 0.2 97.0 ± 0.2 98.8 ± 0.2 98.8 ± 0.1
dataset, communities (departments) with less than 30 nodes are excluded and finally 10 communities are retained. We utilize the node embeddings learned by each model (the latent positions zi for DLSM) to conduct K-means clustering with ground-truth community numbers. The commonly used macro and micro F1 scores are leveraged as evaluation metrics, which are computed using the most likely mappings between true and predicted clusters. Table 2 shows that our proposed model significantly outperforms all other baselines for graph representation learning. Such results verify the advantages of leveraging the Euclidean distance rather than inner product to compute the divergences between nodes for community detection. (See Appendix G for visualizations of node embeddings using a 2D t-SNE projectionthe distributions of the reverse node random factors, which shows a consistent power-law with node degrees and offers more interpretability.)
Ablation study To justify the effectiveness of each component in our framework, including the community membership si, node random factors, δi, γi, and the HLSM decoder architecture, we compare the proposed DLSM with several variants of our method, which are named by the corresponding ablated components. For example, -∆-Γ indicates the variant removed node random factors (see Appendix F for more details). The link prediction results are presented in Table 3 and the community detection results are in Table 8 of Appendix. Our DLSM performs best on all datasets for both tasks, verifying that all components are useful under the proposed framework. Interestingly, from the first two lines of Table 3, we find that the community membership seems ineffective when the node random factors are not involved. This reflects that the binary si mainly benefits our model by shrinking the redundant latent variables to 0 when there are too many parameters to optimize.

heading: 7 CONCLUSION
text: We establish a deep generative model to learn multiple highly interpretable node representations for directed graphs. Our proposed DLSM, comprised of a deep GCN encoder and a HLSM decoder, combines the traditional random graphs models with deep learning based methods to take the complementary advantages of interpretability and scalability. Series of experiments have shown that the model is effective for fitting directed graphs and achieves the state-of-the-art performance on link prediction and community detection. In addition, the interpretable node representations learned by the model can naturally represent both the community structure and degree heterogeneity of complex directed graphs. In the future, we shall extend our model for the more complicated scenes such as weighted or dynamic graphs. While the former with multi-valued edges can be simply achieved, the latter demands for a more efficient method to learn the evolutionary topology of graphs.

heading: A LOSS FUNCTIONS
text: The KL divergence between Normal distributions Kingma & Welling (2014b) for z(l)i is
KL [ qφ ( z (l) i )∥∥∥pθ (z(l)i )] = Gl∑ g=1 −1 2 ( 1 + log ( σ̂ (l) ig )2 − ( µ̂ (l) ig )2 − ( σ̂ (l) ig )2) . (17)
The KL divergence between Binary Concrete distributions Maddison et al. (2017) for s(l)i is
KL [ qφ ( s (l) i )∥∥∥pθ (s(l)i )] = Gl∑ g=1 log π̂ (l) ig
( 1 + π
(l) ig e −λs̃(l)ig )2 π (l) ig ( 1 + π̂ (l) ig e −λs̃(l)ig )2 , (18)
where s̃(l)ig is the log odds of s (l) ig .
The KL divergence between Gamma distributions Joo et al. (2020) for γ(l)i (the same for δ (l) i ) is
KL [ qφ ( γ (l) i )∥∥∥pθ (γ(l)i )] = Gl∑ g=1 log Γ (ξig) Γ ( ξ̂ig ) + (ξ̂ig − ξig)ψ (ξ̂ig) , (19) where Γ(·) and ψ(·) denote the Gamma and digamma function, respectively.

heading: B DETAILS OF THEOREMS
text: B.1 PROOF OF THEOREM 1
For any positive value v > 0, the posterior probability density of rij at v is
f ( rij = v ∣∣A,X) = ∫∫ f ( rij = v,zi,Θ
∣∣A,X) dzidΘ = ∫∫ f ( A,X ∣∣zi,Θ) f (rij = v,zi,Θ) f (A,X) dzidΘ
=
∫∫ f ( A,X ∣∣zi,Θ) f (rij = 0, zi,Θ) f (A,X) · f (rij = v, zi,Θ) f (rij = 0, zi,Θ) dzidΘ
=
∫∫ f (rij = v,zi,Θ)
f (rij = 0, zi,Θ) · f ( rij = 0, zi,Θ ∣∣A,X) dzidΘ = ∫∫ f ( zi ∣∣rij = v,Θ) f (rij = v) f (Θ)
f ( zi ∣∣rij = 0,Θ) f (rij = 0) f (Θ)
· f ( rij = 0 ∣∣A,X) f (zi,Θ∣∣rij = 0,A,X) dzidΘ. (20)
Assume that the true posterior of zi can be approximated by the variational Normal distribution. Then, according to Eq. (12) and (13), we have
f ( zi ∣∣rij = v,Θ) f (rij = v)
f ( zi ∣∣rij = 0,Θ) f (rij = 0)
= 1− p0 p0 exp −12 D∏ d 1 σ2id zi − z̄i − veij − ∑ k∈Ni,k 6=j rikeik ′zi − z̄i − veij −
∑ k∈Ni,k 6=j rikeik
− zi − z̄i − ∑
k∈Ni,k 6=j
rikeik ′zi − z̄i − ∑ k∈Ni,k 6=j rikeik − λv 
= 1− p0 p0 exp −12 D∏ d 1 σ2id v2 − 2v zi − z̄i − ∑ k∈Ni,k 6=j rikeik ′ eij − λv  =
1− p0 p0 exp −12 D∏ d 1 σ2id v − zi − z̄i − ∑ k∈Ni,k 6=j rikeik ′ eij − λ D∏ d σ2id 2
+ 1
2 D∏ d 1 σ2id zi − z̄i − ∑ k∈Ni,k 6=j rikeik ′ eij − λ D∏ d σ2id 2  .
By dividing f(rij = 0|A,X) on both sides of Eq. (20) and integrating over v on (0,+∞), we can obtain
OR(rij > 0|A,X)
=
∫∫ f ( zi ∣∣rij = v,Θ) f (rij = v)
f ( zi ∣∣rij = 0,Θ) f (rij = 0) · f
( zi,Θ ∣∣∣rij = 0,A,X) dzidΘ = ∫∫ g (zi,Θ) f ( zi,Θ
∣∣∣rij = 0,A,X) dzidΘ =E [ g (zi,Θ)
∣∣∣rij = 0,A,X] , where
g (zi,Θ) =
(1− p0)Φ (∏D d σ −1 id [( zi − z̄i − ∑ k∈Ni,k 6=j rikeik )′ eij − λ ∏D d σ 2 id ]) p0φ (∏D d σ −1 id [( zi − z̄i − ∑ k∈Ni,k 6=j rikeik )′ eij − λ ∏D d σ 2 id
]) , and Φ(·) and φ(·) denote the cumulative distribution function and probability density function of the standard Normal distribution, respectively.
The conclusion is proved.
B.2 PROOF OF THEOREM 2
We start our proof by giving two definitions, the standard distance (SD) and actual distance (AD) between node i and neighbor j. Consider a simple case, where the influence discrepancies of all neighbors are neglected and the latent position zi is only dependent on the distances between the node and neighbors (see Fig. 2(b) for an illustration). In this case, the residual i = 0 and zi will coincide with z̄i, thus the distance from neighbor j to node i is equal to ‖zj − z̄i‖, which we refer as the standard distance.
Now we reconsider the influence discrepancies of neighbors upon the node. According to Eq. (1), the actual distance between node i and neighbor j when involving node influences is node influences in
our DLSM are involved by rescaling the standard distance. Assume the vectors from the z̄i to all influential neighbors (rij > 0) are orthogonal to each other, and the random factors are the same at all dimensions of the latent space, denoted as γi and δi, respectively. Then, the actual distance between zi and zj is
ADij = (βoutγi + βinδj)‖ (zi − zj) ‖. (21)
The coefficients βoutγi + βinδj enlarge or shrink SDij = ‖zj − z̄i‖ to ADij using the node activity and popularity factors. Let θ ∈ (0, π/2) be the angle between zi − z̄i and zj − z̄i, i.e.
θ = arccos ( (zi − z̄i)′(z̄j − z̄i) ‖zi − z̄i‖‖zj − z̄i‖ ) . (22)
Since eij is orthogonal to other residual components, the modulus of i can be denoted as ‖ i‖ = rij/ cos θ. According to the law of cosines, we have
AD2ij =SD 2 ij + ‖ i‖2 − 2SD‖ i‖ cos θ
=SD2ij + r2ij
cos2 θ − 2SDijrij . (23)
By combining Eq. (21) and Eq. (23), we obtain
‖zj − z̄i‖2 + r2ij
cos2 θ − 2‖zj − z̄i‖rij =(βoutγi + βinδj)2‖ (zi − zj) ‖2( rij
cos θ − cos θ‖zj − z̄i‖
)2 = [ (βoutγi + βinδj) 2 − sin2 θ ] ‖ (zi − zj) ‖2
rij cos θ
− cos θ‖zj − z̄i‖ =± [ (βoutγi + βinδj) 2 − sin2 θ ] 1 2 ‖zj − zi‖.
Given θ ∈ (0, π/2) (when θ ∈ (π/2, π), rij = 0), the sign for the right side depends on the positive or negative correlation for ‖ i‖ and AD given θ, which is determined by the angle between z̄i − zj and zi − zj , denoted as
o = arccos
( (βoutγi + βinδj) (zi − zj)′(z̄i − zj)
(βoutγi + βinδj) ‖zi − zj‖2 ) = arccos
( z̄i − zj zi − zj ) . (24)
The sign is positive when o ∈ (0, π/2− θ) and negative when o ∈ (π/2− θ, π). In summary, we have the following conclusion, i.e.
rij =h(γi, δj , θ, o)‖zj − zi‖,
where
h(γi, δj , θ, o) =  cos 2 θ − cos θ [ (βoutγi + βinδj) 2 − sin2 θ ] 1 2 , if o ∈ ( 0, π2 − θ ) cos2 θ + cos θ [ (βoutγi + βinδj) 2 − sin2 θ ] 1 2 , if o ∈ ( π 2 − θ, π ) ,
and θ, o is defined by Eq. (22) and Eq. (24), respectively.
The conclusion is proved.

heading: C DIRECTED GRAPH DATASETS
text: Here we provide more details about the directed graph datasets employed in our experiments. Emails network consists of the members from 42 departments of a European research institution (Leskovec et al., 2007). An edge exists if a person sends at least one email to another. Political blogs is a well-studied social network composed of U.S. political blog nodes (Adamic & Glance, 2005), which are labeled as “liberal” and “conservative” clusters. Two blog pages are connected if one is referenced by another. Cora is a citation network of scientific publications which are classified into seven
classes (Sen et al., 2008). Each node and directed edge represent a paper and a citation, respectively. WikiVote is a network of Wikipedia users (Leskovec et al., 2010), where each edge represents a user voting on another to become the administrator. Google is a network of web pages which are connected by hyperlinks (Palla et al., 2007).
The descriptive statistics of the directed graph datasets for link prediction and community detection are summarized in Table 4.

heading: D HYPERPARAMETER SETTING AND SENSITIVITY ANALYSIS
text: The proposed DLSM employs a hierarchical VAE architecture composed of a GCN encoder and HLSM decoder. In practice, both of the two parts have three layers, and the layer sizes are fixed to be 32/64/128 and 50/100/50 for the encoder and decoder, respectively. Note that the last layer of the decoder (l = L) is a linear full connection layer which transforms the sparse latent variables to dense for adjacency reconstruction. Fig. 3(a) shows the link prediction results for increasing numbers of layers, and the performance of our method is continuously rising until more than 4 layers.
The HLSM stochastic decoder requires prior distributions to generate latent variables. Specifically, we employ a standard Normal prior for the latent positions z(l)i , i.e. the variance σ (l) i = 1. For the community membership s(l)i , we use the IBP prior to infer community numbers given the truncation Gl (i.e. layer sizes of the decoder) and a stick-breaking parameters v ∈ (0, 1). Typically, a smaller v results in fewer active communities (i.e. non-zero elements of s(l)i ). Fig. 3(b) presents the link prediction results of our method with different values of v on the Political blogs dataset, and our method performs best when v ∈ (0.9, 1). Therefore, we set v = 0.9 to ensure the model fully exploring the potential community structure of graphs. Last, we adopt flat priors for the Dirichlet node random factors γi and δi, where all parameters are set as ξ (l) ig = ψ (l) ig = 1/Gl.
For inference, the temperature parameter λ of the Binary Concrete distribution is fixed to 1. All models are trained by 1,000-2,000 iterations with a learning rate of 0.01 on an RTX 2080 Ti GPU.

heading: E LINK PREDICTION ON UNDIRECTED GRAPHS
text: Undirected graphs can be regarded as special cases in our framework, as an undirectional edge is treated as a bidirectional one. In this section we provide some experimental results for link prediction on three real-world undirected graph datasets. The baselines are the same as those in the experiments of directed graphs.
E.1 DATASETS
We consider three real-world undirected graph datasets, namely NIPS12, Cora and Pubmed. NIPS12 is a coauthor network of the authors in NIPS papers from volumes 1-12 Zhou (2015). Cora is the same dataset used in the experiments for directed link prediction, except that the directions of edges are neglected here. Pubmed is a citation network of scientific publications Sen et al. (2008). Both Cora and Pubmed datasets contain sparse bag-of-words feature matrices, which are used as node attributes. The descriptive statistics of the datasets are summarized in Table 5.
E.2 LINK PREDICTION RESULTS
The results of link prediction on undirected graphs are presented in Table 6. On all datasets, our DLSM achieves the best or comparative performances with the state-of-the-art baselines for undirected graphs.

heading: F MORE DETAILS FOR ABLATION STUDY
text: The community detection results of ablation study is given in Table 8. Generally, models including the community membership perform better than those without this component, and the full version DLSM reaches the best performances on all datasets.
G VISUALIZATIONS OF NODE REPRESENTATIONS
Latent position We leverage a 2D t-SNE projection (Van der Maaten & Hinton, 2008) to visualize the learned latent positions for the Emails dataset. As comparisons, we also illustrate the node representations learned by DGLFRM and LGVG, both of which also consider the community structure but overlook the degree heterogeneity of graphs. The transformed latent variables learned by these three models are plotted in Fig. 4. It is clearly seen that our model performs best in fitting such directed graphs.
DLSM 72.0 ± 1.0 77.0 ± 1.2 92.5 ± 0.7 92.6 ± 0.7 24.6 ± 0.8 25.7 ± 0.5
Community membership The IBP prior enables the community membership s(l)i to freely detect the true number of communities. In the experiments we set the truncation Kl to be 50/100, which are far larger than the true numbers of communities. Therefore, the hyperparameter (truncation) actually exert little impact for community detection. In Fig. 5 we illustrate the learned s(1)i with IBP and independent Bernoulli priors for fitting the Emails network, of which most nodes belong to one of the 10 communities. A lighter color indicate a larger probability for the corresponding dimension (community) of latent variables to be activated ( 6= 0). It shows that IBP prior condenses the activated communities to be the first 10-15 dimensions, which is relatively approximate to the true number of communities, while the independent Bernoulli variables are much more chaotic and uniformly distributed at all dimensions. Therefore, the proposed IBP prior provide better interpretability for the community membership, in addition to the improvement of model performances.
Reverse node random factors The pairwise node random factors γi and δi are supposed to measure the heterogeneity of out-degrees and in-degrees, respectively, which typically follow the power-law (Barabási & Albert, 1999). Fig. 6(a) and (b) present the probability density distributions (PDD) of the node degrees and reverse random factors learned by DLSM on the political blogs network. It seems that the degree distributions are well fitted by the reverse node random factors, indicating that our DLSM can well represent the degree heterogeneity via these latent variables. Furthermore, Fig. 6(c) illustrates the complementary cumulative distributions (CCD) of the random factors. As can be seen, the logarithm CCD of both γi and δi are approximately linear, with different slopes though. This shows that the proposed Dirichlet latent variables are flexible enough to accommodate the power-law distribution of degrees.

Review summaries:
This paper proposed a directed graphical model called deep latent space model for modeling the graph adjacency matrix generation process, where the latent variables are organized in a hierarchical way with three meaningful components: “latent position”, “community membership” and “node random factors” for each node. The learning leverages the variational framework originated from VAE, with GCN as the core component for posterior parameterization. Experiments on link prediction and community detection show that the proposed method is able to achieve better empirical results.  

This paper studies a deep latent space model for directed graph representation. The model is basically a hierarchical variational auto-encoder architecture and some latent variables like social activity and popularity factors of node account for the edge directions. Some latent variables are interpreted as positions in the latent space or community membership. The model, therefore, is claimed to be interpretable and superior in experiments.  ----Update after rebuttal  The authors basically answered my questions and I'd like to increase the recommendation score. I still think this submission has incremental technical contributions so my new recommendation score is 6 (marginally above the acceptance threshold).

This paper proposes a VAE-based graph representation learning model for directed graphs. The model is composed of a GCN encoder and a hierarchical latent space model decoder. The decoder explicitly models three kinds of node representations: the node latent position vectors $z$, the community membership vectors $s$ and node random factors $\gamma,\delta$. Specifically, the community membership vectors $s$ are used to gate control which of $z$ and $\gamma,\delta$ pass to the next layer. And finally the adjacency matrix is reconstructed by the distance between $z$'s filtered by node influences $\gamma,\delta$. Experiments show better link prediction and community detection performance than plain VGAE and other VGAE methods.

In this paper, the authors proposed a new directed graph representation model. The proposed model is developed from a Bayesian viewpoint and is implemented in the framework of the variational autoencoder. Additionally, the authors justified the interpretability of the model. They indicate that the proposed model quantitatively captures the influence of each node on its neighbors. Experimental results show that the proposed model performs well in commonly-used datasets.