Joint embedding of structure and features via graph convolutional networks

TL;DR

This paper introduces AN2VEC, a graph convolutional network-based method that disentangles and captures shared and individual information from node features and network structure, improving understanding of social and biological networks.

Contribution

It presents a novel multitask GCN Variational Autoencoder that separates feature, structure, and shared information in node embeddings, enhancing interpretability and predictive power.

Findings

Shared information improves embedding performance in correlated networks.

Embedding captures joint feature-structure information effectively.

Performance increases with higher feature-structure correlation.

Abstract

The creation of social ties is largely determined by the entangled effects of people's similarities in terms of individual characters and friends. However, feature and structural characters of people usually appear to be correlated, making it difficult to determine which has greater responsibility in the formation of the emergent network structure. We propose \emph{AN2VEC}, a node embedding method which ultimately aims at disentangling the information shared by the structure of a network and the features of its nodes. Building on the recent developments of Graph Convolutional Networks (GCN), we develop a multitask GCN Variational Autoencoder where different dimensions of the generated embeddings can be dedicated to encoding feature information, network structure, and shared feature-network information. We explore the interaction between these disentangled characters by comparing the embedding reconstruction performance to a baseline case where no shared information is extracted. We use synthetic datasets with different levels of interdependency between feature and network characters and show (i) that shallow embeddings relying on shared information perform better than the corresponding reference with unshared information, (ii) that this performance gap increases with the correlation between network and feature structure, and (iii) that our embedding is able to capture joint information of structure and features. Our method can be relevant for the analysis and prediction of any featured network structure ranging from online social systems to network medicine.

Figures

Tables

Table 1. Table 1: Link prediction task in citation networks. SC, DW, GAE and VGAE values are from [ 30 ] . Error values indicate the sample standard deviation.

Method	Cora		CiteSeer		PubMed
	AUC	AP	AUC	AP	AUC	AP
SC	84.6 $\pm$ 0.01	88.5 $\pm$ 0.00	80.5 $\pm$ 0.01	85.0 $\pm$ 0.01	84.2 $\pm$ 0.02	87.8 $\pm$ 0.01
DW	83.1 $\pm$ 0.01	85.0 $\pm$ 0.00	80.5 $\pm$ 0.02	83.6 $\pm$ 0.01	84.4 $\pm$ 0.00	84.1 $\pm$ 0.00
GAE	91.0 $\pm$ 0.02	92.0 $\pm$ 0.03	89.5 $\pm$ 0.04	89.9 $\pm$ 0.05	96.4 $\pm$ 0.00	96.5 $\pm$ 0.00
VGAE	91.4 $\pm$ 0.01	92.6 $\pm$ 0.01	90.8 $\pm$ 0.02	92.0 $\pm$ 0.02	94.4 $\pm$ 0.02	94.7 $\pm$ 0.02
AN2VEC-0	89.5 $\pm$ 0.01	90.6 $\pm$ 0.01	91.2 $\pm$ 0.01	91.5 $\pm$ 0.02	91.8 $\pm$ 0.01	93.2 $\pm$ 0.01
AN2VEC-16	89.4 $\pm$ 0.01	90.2 $\pm$ 0.01	91.1 $\pm$ 0.01	91.3 $\pm$ 0.02	92.1 $\pm$ 0.01	92.8 $\pm$ 0.01
AN2VEC-S-0	92.9 $\pm$ 0.01	93.4 $\pm$ 0.01	94.3 $\pm$ 0.01	94.8 $\pm$ 0.01	95.1 $\pm$ 0.01	95.4 $\pm$ 0.01
AN2VEC-S-16	93.0 $\pm$ 0.01	93.5 $\pm$ 0.00	94.9 $\pm$ 0.00	95.1 $\pm$ 0.00	93.1 $\pm$ 0.01	93.1 $\pm$ 0.01

Equations

$$\mathbf{\mu }={\mathrm{GCN}}_{\mathbf{\mu }}(\mathbf{X},\mathbf{A})\text{and}\log \mathbf{\sigma }={\mathrm{GCN}}_{\mathbf{\sigma }}(\mathbf{X},\mathbf{A}).$$

$${\mathrm{GCN}}_{\alpha }(\mathbf{X},\mathbf{A})=\hat{\mathbf{A}}\mathrm{ReLU}(\hat{\mathbf{A}}\mathbf{X}{\mathbf{W}}_0^{enc}){\mathbf{W}}_{1,\alpha }^{enc}$$

$$\hat{\mathbf{A}}$$

$$\tilde{\mathbf{A}}$$

$${\tilde{D}}_{ii}$$

$${\mathbf{\xi }}_i|\mathbf{A},\mathbf{X}\sim \mathcal{N}({\mathbf{\mu }}_i,\mathrm{diag}({\mathbf{\sigma }}_i^2)).$$

$$q(\mathbf{\xi }|\mathbf{X},\mathbf{A})=\prod _{i=1}^{N}q({\mathbf{\xi }}_i|\mathbf{A},\mathbf{X}).$$

$$A_{ij}|{\mathbf{\xi }}_i,{\mathbf{\xi }}_j$$

$$\mathrm{MLB}(\mathbf{\xi })$$

$$\mathbf{\gamma }$$

$${\mathbf{X}}_i|{\mathbf{\xi }}_i$$

$$\mathrm{MLP}(\mathbf{\xi })$$

$$p(\mathbf{X},\mathbf{A}|\mathbf{\xi })$$

$$p(\mathbf{A}|\mathbf{\xi })$$

$$p(\mathbf{X}|\mathbf{\xi })$$

$$-\log p(\mathbf{A},\mathbf{X})$$

$$=D_{KL}(q(\mathbf{\xi }|\mathbf{A},\mathbf{X})||\mathcal{N}(0,{\mathbf{I}}_F))$$

$$-\text{mathds}{E}_{q(\mathbf{\xi }|\mathbf{A},\mathbf{X})}[\log (p(\mathbf{A},\mathbf{X}|\mathbf{\xi },\mathbf{\theta })p(\mathbf{\theta }))]$$

$$={\mathcal{L}}_{KL}+{\mathcal{L}}_{\mathbf{A}}+{\mathcal{L}}_{\mathbf{X}}+{\mathcal{L}}_{\mathbf{\theta }}$$

$${\mathcal{L}}_{\mathbf{A}}$$

$$\simeq -\frac{1}{K}\sum _{k=1}^K\sum _{i,j=1}^N[A_{ij}\log (\mathrm{MLB}({\mathbf{\xi }}^{(k)}{)}_{ij})$$

$$+(1-A_{ij})\log (1-\mathrm{MLB}({\mathbf{\xi }}^{(k)}{)}_{ij})]$$

$${\mathcal{L}}_{\mathbf{X}}$$

$$\simeq -\frac{1}{K}\sum _{k=1}^K\sum _{i=1}^N\sum _{j=1}^D{X}_{ij}\log (\mathrm{MLP}({\mathbf{\xi }}^{(\mathbf{k})}{)}_{ij})$$

$${\mathcal{L}}_{KL}=\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^F{\mu }_{ij}^2+{\sigma }_{ij}^2-2\log {\sigma }_{ij}-1$$

$${\tilde{\mathcal{L}}}_{\mathbf{A}}$$

$$+\frac{1-A_{ij}}{1-d}\log (1-\mathrm{MLB}({\mathbf{\xi }}^{(k)}{)}_{ij})]$$

$$\mathcal{L}=\frac{{\tilde{\mathcal{L}}}_{\mathbf{A}}}{N^2\log 2}+\frac{{\mathcal{L}}_{\mathbf{X}}}{N\log D}+\frac{{\mathcal{L}}_{KL}}{NF{\kappa }_{KL}}+\frac{||\mathbf{\theta }|{|}_2^2}{2{\kappa }_{\mathbf{\theta }}}$$

$${\mathbf{\xi }}_{i,a:b}$$

$${\mathbf{\xi }}_{a:b}$$

$$A_{ij}|{\mathbf{\xi }}_{i,1:F_{\mathbf{A}}+F_{\mathbf{AX}}},{\mathbf{\xi }}_{j,1:F_{\mathbf{A}}+F_{\mathbf{AX}}}\sim \mathrm{Ber}(\mathrm{MLB}({\mathbf{\xi }}_{1:F_{\mathbf{A}}+F_{\mathbf{AX}}}{)}_{ij})$$

$${\mathbf{X}}_i|{\mathbf{\xi }}_{i,F_{\mathbf{A}}+1:F}\sim \mathrm{Multinomial}(1,\mathrm{MLP}({\mathbf{\xi }}_{F_{\mathbf{A}}+1:F}{)}_i)$$

$$\tilde{\mathbf{\mu }}={\mathbf{\mu }}_{1:F_{\mathbf{A}}}$$

$$\displaystyle{}\operatorname*{\scalebox{1.0}[1.5]{$\parallel$}}\bm{\mu}_{F_{\mathbf{A}}+2F_{\mathbf{AX}}+1:F_{\mathbf{A}}+2F_{\mathbf{AX}}+F_{\mathbf{X}}}$$

$$\log \tilde{\mathbf{\sigma }}=\log {\mathbf{\sigma }}_{1:F_{\mathbf{A}}}$$

$$\displaystyle{}\operatorname*{\scalebox{1.0}[1.5]{$\parallel$}}\log\bm{\sigma}_{F_{\mathbf{A}}+2F_{\mathbf{AX}}+1:F_{\mathbf{A}}+2F_{\mathbf{AX}}+F_{\mathbf{X}}}$$

Full text

1 Introduction

Although it is relatively easy to obtain the proxy social network and various individual features for users of online social platforms, the combined characterisation of these types of information is still challenging our methodology. While current approaches have been able to approximate the observed marginal distributions of node and network features separately, their combined consideration was usually done via summary network statistics merged with otherwise independently built feature sets of nodes. However, the entanglement between structural patterns and feature similarities appears to be fundamental to a deeper understanding of network formation and dynamics. The value of this joint information then calls for the development of statistical tools for the learning of combined representation of network and feature information and their dependencies.

The formation of ties and mesoscopic structures in online social networks is arguably determined by several competing factors. Considering only network information, neighbour similarities between people are thought to explain network communities, where triadic closure mechanisms [1, 2] induce ties between peers with larger fractions of common friends [3]. Meanwhile, random bridges [2] are built via focal closure mechanisms, optimising the structure for global connectedness and information dissemination. At the same time, people in the network can be characterised by various individual features such as their socio-demographic background [4, 5], linguistic characters [6, 5, 7], or the distributions of their topics of interests [8, 9], to mention just a few. Such features generate homophilic tie creation preferences [10, 11], which induce links with higher probability between similar individuals, whom in turn form feature communities of shared interest, age, gender, or socio-economic status, and so on [4, 12]. Though these mechanisms are not independent and lead to correlations between feature and network communities, it is difficult to define the causal relationship between the two: first, because simultaneously characterising similarities between multiple features and a complex network structure is not an easy task; second, because it is difficult to determine, which of the two types of information, features or structure, is driving network formation to a greater extent. Indeed, we do not know what fraction of similar people initially get connected through homophilic tie creation, versus the fraction that first get connected due to structural similarities before influencing each other to become more similar [13, 14].

Over the last decade popular methods have been developed to characterise structural and feature similarities and to identify these two notions of communities. The detection of network communities has been a major challenge in network science with various concepts proposed [15, 16, 17] to solve it as an unsupervised learning task [18, 19]. Commonly, these algorithms rely solely on network information, and their output is difficult to cross-examine without additional meta-data, which is usually disregarded in their description. On the other hand, methods grouping similar people into feature communities typically ignore network information, and exclusively rely on individual features to solve the problem as a data clustering challenge [20, 21]. Some semi-supervised learning tasks, such as link prediction, may take feature and structural information simultaneously into account, but only by enriching individual feature vectors with node-level network characteristics such as degree or local clustering [22, 23, 9]. Methods that would take higher order network correlations and multivariate feature information into account at the same time are still to be defined. Their development would offer huge potential in understanding the relation between individuals’ characters, their social relationships, the content they are engaged with, and the larger communities they belong to. This would not only provide us with deeper insight about social behaviour, it would give us predictive tools for the emergence of network structure, individual interests and behavioural patterns.

In this paper we propose a contribution to solve this problem by developing a joint feature-network embedding built on multitask Graph Convolutional Networks [24, 25, 26, 27] and Variational Autoencoders (GCN-VAE) [28, 29, 30, 31, 32], which we call the Attributed Node to Vector method (AN2VEC). In our model, different dimensions of the generated embeddings can be dedicated to encode feature information, network structure, or shared feature-network information separately. Unlike previous embedding methods dealing with features [33, 34, 35, 36], this interaction model [37] allows us to explore the dependencies between the disentangled network and feature information by comparing the embedding reconstruction performance to a baseline case where no shared information is extracted. Using this method, we can identify an optimal reduced embedding, which indicates whether combined information coming from the structure and features is important, or whether their non-interacting combination is sufficient for reconstructing the featured network.

In practice, as this method solves a reconstruction problem, it may give important insights about the combination of feature- and structure-driven mechanisms which determine the formation of a given network. As an embedding, it is useful to identify people sharing similar individual and structural characters. And finally, by measuring the optimal overlap between feature– and network-associated dimensions, it can be used to verify network community detection methods to see how well they identify communities explained by feature similarities.

In what follows, after summarising the relevant literature, we introduce our method and demonstrate its performance on synthetic featured networks, for which we control the structural and feature communities as well as the correlations between the two. As a result, we will show that our embeddings, when relying on shared information, outperform the corresponding reference without shared information, and that this performance gap increases with the correlation between network and feature structure since the method can capture the increased joint information. Next, we extensively explore the behaviour of our model on link prediction and node classification on standard benchmark datasets, comparing it to well-known embedding methods. Finally, we close our paper with a short summary and a discussion about potential future directions for our method.

2 Related Work

The advent of increasing computational power coupled with the continuous release and ubiquity of large graph-structured datasets has triggered a surge of research in the field of network embeddings. The main motivation behind this trend is to be able to convert a graph into a low-dimensional space where its structural information and properties are maximally preserved [38]. The aim is to extract unseen or hard to obtain properties of the network, either directly or by feeding the learned representations to a downstream inference pipeline.

2.1 Graph embedding survey: from matrix factorisation to deep learning

In early work, low-dimensional node embeddings were learned for graphs constructed from non-relational data by relying on matrix factorisation techniques. By assuming that the input data lies on a low dimensional manifold, such methods sought to reduce the dimensionality of the data while preserving its structure, and did so by factorising graph Laplacian eigenmaps [39] or node proximity matrices [40].

More recent work has attempted to develop embedding architectures that can use deep learning techniques to compute node representations. DeepWalk [41], for instance, computes node co-occurrence statistics by sampling the input graph via truncated random walks, and adopts a SkipGram neural language model to maximise the probability of observing the neighbourhood of a node given its embedding. By doing so the learned embedding space preserves second order proximity in the original graph. However, this technique and the ones that followed [42, 43] present generalisation caveats, as unobserved nodes during training cannot be meaningfully embedded in the representation space, and the embedding space itself cannot be generalised between graphs. Instead of relying on random walk-based sampling of graphs to feed deep learning architectures, other approaches have used the whole network as input to autoencoders in order to learn, at the bottleneck layer, an efficient representation able to recover proximity information [31, 44, 34]. However, the techniques developed herein remained limited due to the fact that successful deep learning models such as convolutional neural networks require an underlying euclidean structure in order to be applicable.

2.2 Geometric deep learning survey: defining convolutional layers on non-euclidean domains

This restriction has been gradually overcome by the development of graph convolutions or Graph Convolutional Networks (GCN). By relying on the definition of convolutions in the spectral domain, Bruna et al. [25] defined spectral convolution layers based on the spectrum of the graph Laplacian. Several modifications and additions followed and were progressively added to ensure the feasibility of learning on large networks, as well as the spatial localisation of the learned filters [45, 27]. A key step is made by [46] with the use of Chebychev polynomials of the Laplacian, in order to avoid having to work in the spectral domain. These polynomials, of order up to $r$, generate localised filters that behave as a diffusion operator limited to $r$ hops around each vertex. This construction is then further simplified by Kipf and Welling by assuming among others that $r\approx 2$ [24].

Recently, these approaches have been extended into more flexible and scalable frameworks. For instance, Hamilton et al. [26] extended the original GCN framework by enabling the inductive embedding of individual nodes, training a set of functions that learn to aggregate feature information from a node’s local neighborhood. In doing so, every node defines a computational graph whose parameters are shared for all the graphs nodes.

More broadly, the combination of GCN with autoencoder architectures has proved fertile for creating new embedding methods. The introduction of probabilistic node embeddings, for instance, has appeared naturally from the application of variational autoencoders to graph data [29, 28, 30], and has since led to explorations of the uncertainty of embeddings [36, 32], of appropriate levels of disentanglement and overlap [47], and of better representation spaces for measuring pairwise embedding distances (see in particular recent applications of the Wasserstein distance between probabilistic embeddings [32, 48]). Such models consistently outperform earlier techniques on different benchmarks and have opened several interesting lines of research in fields ranging from drug design [49] to particle physics [50]. Most of the more recent approaches mentioned above can incorporate node features (either because they rely on them centrally, or as an add-on). However, with the exception of DANE [33], they mostly do so by assuming that node features are an additional source of information, which is congruent with the network structure (e.g. multi-task learning with shared weights [34], or fusing both information types together [35]). That assumption may not hold in many complex datasets, and it seems important to explore what type of embeddings can be constructed when we lift it, considering different levels of congruence between a network and the features of its nodes.

We therefore set out to make a change to the initial GCN-VAE in order to: (i) create embeddings that are explicitly trained to encode both node features and network structure; (ii) make it so that these embeddings can separate the information that is shared between network and features, from the (possibly non-congruent) information that is specific to either network or features; and (iii) be able to tune the importance that is given to each type of information in the embeddings.

3 Methods

In this section we present the architecture of the neural network model we use to generate shared feature-structure node embeddings111The implementation of our model is available online at github.com/ixxi-dante/an2vec.. We take a featured network as input, with structure represented as an adjacency matrix and node features represented as vectors (see below for a formal definition). Our starting point is a GCN-VAE, and our first goal is a multitask reconstruction of both node features and network adjacency matrix. Then, as a second goal, we tune the architecture to be able to scale the number of embedding dimensions dedicated to feature-only reconstruction, adjacency-only reconstruction, or shared feature-adjacency information, while keeping the number of trainable parameters in the model constant.

3.1 Multitask graph convolutional autoencoder

We begin with the graph-convolutional variational autoencoder developed by [30], which stacks graph-convolutional (GC) layers [24] in the encoder part of a variational autoencoder [29, 28] to obtain a lower dimensional embedding of the input structure. This embedding is then used for the reconstruction of the original graph (and in our case, also of the features) in the decoding part of the model. Similarly to [24], we use two GC layers in our encoder and generate Gaussian-distributed node embeddings at the bottleneck layer of the autoencoder. We now introduce each phase of our embedding method in formal terms.

3.1.1 Encoder

We are given an undirected unweighted feautured graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, with $N=|\mathcal{V}|$ nodes, each node having a $D$-dimensional feature vector. Loosely following the notations of [30], we note $\mathbf{A}$ the graph’s $N\times N$ adjacency matrix (diagonal elements set to [math]), $\mathbf{X}$ the $N\times D$ matrix of node features, and $\mathbf{X}_{i}$ the D-dimensional feature vector of a node $i$.

The encoder part of our model is where $F$-dimensional node embeddings are generated. It computes $\bm{\mu}$ and $\bm{\sigma}$, two $N\times F$ matrices, which parametrise a stochastic embedding of each node:

[TABLE]

Here we use two graph-convolutional layers for each parameter set, with shared weights at the first layer and parameter-specific weights at the second layer:

[TABLE]

In this equation, $W^{enc}_{0}$ and $W^{enc}_{1,\alpha}$ are the weight matrices for the linear transformations of each layer’s input; $\operatorname{ReLU}$ refers to a rectified linear unit [51]; and following the formalism introduced in [24], $\hat{\mathbf{A}}$ is the standard normalised adjacency matrix with added self-connections, defined as:

[TABLE]

where $\mathbf{I}_{N}$ is the $N\times N$ identity matrix.

3.1.2 Embedding

The parameters $\bm{\mu}$ and $\bm{\sigma}$ produced by the encoder define the distribution of an $F$-dimensional stochastic embedding $\bm{\xi}_{i}$ for each node $i$, defined as:

[TABLE]

Thus, for all the nodes we can write a probability density function over a given set of embeddings $\bm{\xi}$, in the form of an $N\times F$ matrix:

[TABLE]

3.1.3 Decoder

The decoder part of our model aims to reconstruct both the input node features and the input adjacency matrix by producing parameters of a generative model for each of the inputs. On one hand, the adjacency matrix $\mathbf{A}$ is modelled as a set of independent Bernoulli random variables, whose parameters come from a bi-linear form applied to the output of a single dense layer:

[TABLE]

Similarly to above, $W^{dec}_{\mathbf{A},0}$ is the weight matrix for the first adjacency matrix decoder layer, and $W^{dec}_{\mathbf{A},1}$ is the weight matrix for the bilinear form which follows.

On the other hand, features can be modelled in a variety of ways, depending on whether they are binary or continuous, and if their norm is constrained or not. Features in our experiments are one-hot encodings, so we model the reconstruction of the feature matrix $\mathbf{X}$ by using $N$ single-draw $D$-categories multinomial random variables. The parameters of those multinomial variables are computed from the embeddings with a two-layer perceptron:222Other types of node features are modelled according to their constraints and domain. Binary features are modelled as independent Bernoulli random variables. Continuous-range features are modelled as Gaussian random variables in a similar way to the embeddings themselves.

[TABLE]

In the above equations, $\operatorname{sigmoid}(z)=\frac{1}{1+e^{-z}}$ refers to the logistic function applied element-wise on vectors or matrices, and $\operatorname{softmax}(\mathbf{z})_{i}=\frac{e^{z_{i}}}{\sum_{j}e^{z_{j}}}$ refers to the normalised exponential function, also applied element-wise, with $j$ running along the rows of matrices (and along the indices of vectors).

Thus we can write the probability density for a given reconstruction as:

[TABLE]

3.1.4 Learning

The variational autoencoder is trained by minimising an upper bound to the marginal likelihood-based loss [29] defined as:

[TABLE]

Here $\mathcal{L}_{KL}$ is the Kullback-Leibler divergence between the distribution of the embeddings and a Gaussian Prior, and $\bm{\theta}$ is the vector of decoder parameters whose associated loss $\mathcal{L}_{\bm{\theta}}$ acts as a regulariser for the decoder layers.333Indeed, following [29] we assume $\bm{\theta}\sim\mathcal{N}(0,\kappa_{\bm{\theta}}\mathbf{I})$, such that $\mathcal{L}_{\bm{\theta}}=-\log p(\bm{\theta})=\frac{1}{2}\dim(\bm{\theta})\log(2\pi\kappa_{\bm{\theta}})+\frac{1}{2\kappa_{\bm{\theta}}}||\bm{\theta}||_{2}^{2}$. Computing the adjacency and feature reconstruction losses by using their exact formulas is computationally not tractable, and the standard practice is instead to estimate those losses by using an empirical mean. We generate $K$ samples of the embeddings by using the distribution $q(\bm{\xi}|\mathbf{A},\mathbf{X})$ given by the encoder, and average the losses of each of those samples444In practice, $K=1$ is often enough. [29, 28]:

[TABLE]

Finally, for diagonal Gaussian embeddings such as the ones we use, $\mathcal{L}_{KL}$ can be expressed directly [28]:

[TABLE]

3.1.5 Loss adjustments

In practice, to obtain useful results a few adjustments are necessary to this loss function. First, given the high sparsity of real-world graphs, the $A_{ij}$ and $1-A_{ij}$ terms in the adjacency loss must be scaled respectively up and down in order to avoid globally near-zero link reconstruction probabilities. Instead of penalising reconstruction proportionally to the overall number of errors in edge prediction, we want false negatives ($A_{ij}$ terms) and false positives ($1-A_{ij}$ terms) to contribute equally to the reconstruction loss, independent of graph sparsity. Formally, let $d=\frac{\sum_{ij}A_{ij}}{N^{2}}$ denote the density of the graph’s adjacency matrix ($d=\frac{N-1}{N}\times\operatorname{density}(\mathcal{G})$); then we replace $\mathcal{L}_{\mathbf{A}}$ by the following re-scaled estimated loss (the so-called “balanced cross-entropy”):

[TABLE]

Second, we correct each component loss for its change of scale when the shapes of the inputs and the model parameters change: $\mathcal{L}_{KL}$ is linear in $N$ and $F$, $\tilde{\mathcal{L}}_{\mathbf{A}}$ is quadratic in $N$, and $\mathcal{L}_{\mathbf{X}}$ is linear in $N$ (but not in $F$, remember that $\sum_{j}X_{ij}=1$ since each $\mathbf{X}_{i}$ is a single-draw multinomial).

Beyond dimension scaling, we also wish to keep the values of $\tilde{\mathcal{L}}_{\mathbf{A}}$ and $\mathcal{L}_{\mathbf{X}}$ comparable and, doing so, maintain a certain balance between the difficulty of each task. As a first approximation to the solution, and in order to avoid more elaborate schemes which would increase the complexity of our architecture (such as [52]), we divide both loss components by their values at maximum uncertainty555That is, $p(A_{ij}|\bm{\xi},\bm{\theta})=\frac{1}{2}\quad\forall i,j$, and $p(X_{ij}|\bm{\xi},\bm{\theta})=\frac{1}{D}\quad\forall i,j$., respectively $\log 2$ and $\log D$.

Finally, we make sure that the regulariser terms in the loss do not overpower the actual learning terms (which are now down-scaled close to 1) by adjusting $\kappa_{\bm{\theta}}$ and an additional factor, $\kappa_{KL}$, which scales the Kullback-Leibler term.666We use $\kappa_{KL}=2\kappa_{\bm{\theta}}=10^{3}$. These adjustments lead us to the final total loss the model is trained for:

[TABLE]

where we have removed constant terms with respect to trainable model parameters.

3.2 Scaling shared information allocation

The model we just presented uses all dimensions of the embeddings indiscriminately to reconstruct the adjacency matrix and the node features. While this can be useful in some cases, it cannot adapt to different interdependencies between graph structure and node features; in cases where the two are not strongly correlated, the embeddings would lose information by conflating features and graph structure. Therefore our second aim is to adjust the dimensions of the embeddings used exclusively for feature reconstruction, or for adjacency reconstruction, or used for both.

In a first step, we restrict which part of a node’s embedding is used for each task. Let $F_{\mathbf{A}}$ be the number of embedding dimensions we will allocate to adjacency matrix reconstruction only, $F_{\mathbf{X}}$ the number of dimensions allocated to feature reconstruction only, and $F_{\mathbf{AX}}$ the number of dimensions allocated to both. We have $F_{\mathbf{A}}+F_{\mathbf{AX}}+F_{\mathbf{X}}=F$. We further introduce the following notation for the restriction of the embedding of node $i$ to a set of dedicated dimensions $\{a,\dots,b\}$777Note that the order of the indices does not change the training results, as the model has no notion of ordering inside its layers. What follows is valid for any permutation of the dimensions, and the actual indices only matter to downstream interpretation of the embeddings after training.:

[TABLE]

This extends to the full matrix of embeddings similarly:

[TABLE]

Using these notations we adapt the decoder to reconstruct adjacency and features as follows:

[TABLE]

In other words, adjacency matrix reconstruction relies on $F_{\mathbf{A}}+F_{\mathbf{AX}}$ embedding dimensions, feature reconstruction relies on $F_{\mathbf{X}}+F_{\mathbf{AX}}$ dimensions, and $F_{\mathbf{AX}}$ overlapping dimensions are shared between the two. Our reasoning is that for datasets where the dependency between features and network structure is strong, shallow models with higher overlap value will perform better than models with the same total embedding dimensions $F$ and less overlap, or will perform on par with models that have more total embedding dimensions and less overlap. Indeed, the overlapping model should be able to extract the information shared between features and network structure and store it in the overlapping dimensions, while keeping the feature-specific and structure-specific information in their respective embedding dimensions. This is to compare to the non-overlapping case, where shared network-feature information is stored redundantly, both in feature- and structure-specific embeddings, at the expense of a larger number of distinct dimensions.

Therefore, to evaluate the performance gains of this architecture, one of our measures is to compare the final loss for different hyperparameter sets, keeping $F_{\mathbf{A}}+F_{\mathbf{AX}}$ and $F_{\mathbf{X}}+F_{\mathbf{AX}}$ fixed and varying the overlap size $F_{\mathbf{AX}}$. Now, to make sure the training losses for different hyperparameter sets are comparable, we must maintain the overall number of trainable parameters in the model fixed. The decoder already has a constant number of trainable parameters, since it only depends on the number of dimensions used for decoding features ($F_{\mathbf{X}}+F_{\mathbf{AX}}$) and adjacency matrix ($F_{\mathbf{A}}+F_{\mathbf{AX}}$), which are themselves fixed.

On the other hand, the encoder requires an additional change. We maintain the dimensions of the encoder-generated $\bm{\mu}$ and $\bm{\sigma}$ parameters fixed at $F_{\mathbf{A}}+2F_{\mathbf{AX}}+F_{\mathbf{X}}$ (independently from $F_{\mathbf{AX}}$, given the constraints above), and reduce those outputs to $F_{\mathbf{A}}+F_{\mathbf{AX}}+F_{\mathbf{X}}$ dimensions by averaging dimensions $\{F_{\mathbf{A}}+1,\dots,F_{\mathbf{A}}+F_{\mathbf{AX}}\}$ and $\{F_{\mathbf{A}}+F_{\mathbf{AX}}+1,\dots,F_{\mathbf{A}}+2F_{\mathbf{AX}}\}$ together.888Formally:

$\displaystyle\tilde{\bm{\mu}}=\bm{\mu}_{1:F_{\mathbf{A}}}$ \displaystyle{}\operatorname*{\scalebox{1.0}[1.5]{\parallel}}\frac{1}{2}(\bm{\mu}_{F_{\mathbf{A}}+1:F_{\mathbf{A}}+F_{\mathbf{AX}}}+\bm{\mu}_{F_{\mathbf{A}}+F_{\mathbf{AX}}+1:F_{\mathbf{A}}+2F_{\mathbf{AX}}})

\displaystyle{}\operatorname*{\scalebox{1.0}[1.5]{\parallel}}\bm{\mu}_{F_{\mathbf{A}}+2F_{\mathbf{AX}}+1:F_{\mathbf{A}}+2F_{\mathbf{AX}}+F_{\mathbf{X}}}

$\displaystyle\log\tilde{\bm{\sigma}}=\log\bm{\sigma}_{1:F_{\mathbf{A}}}$ \displaystyle{}\operatorname*{\scalebox{1.0}[1.5]{\parallel}}\frac{1}{2}(\log\bm{\sigma}_{F_{\mathbf{A}}+1:F_{\mathbf{A}}+F_{\mathbf{AX}}}+\log\bm{\sigma}_{F_{\mathbf{A}}+F_{\mathbf{AX}}+1:F_{\mathbf{A}}+2F_{\mathbf{AX}}})

\displaystyle{}\operatorname*{\scalebox{1.0}[1.5]{\parallel}}\log\bm{\sigma}_{F_{\mathbf{A}}+2F_{\mathbf{AX}}+1:F_{\mathbf{A}}+2F_{\mathbf{AX}}+F_{\mathbf{X}}}

where \operatorname*{\scalebox{1.0}[1.5]{\parallel}} denotes concatenation along the columns of the matrices. In turn, this model maintains a constant number of trainable parameters, while allowing us to adjust the number of dimensions $F_{\mathbf{AX}}$ shared by feature and adjacency reconstruction (keeping $F_{\mathbf{A}}+F_{\mathbf{AX}}$ and $F_{\mathbf{X}}+F_{\mathbf{AX}}$ constant). Figure 1 schematically represents this architecture.

4 Results

We are interested in measuring two main effects: first, the variation in model performance as we increase the overlap in the embeddings, and second, the capacity of the embeddings with overlap (versus no overlap) to capture and benefit from dependencies between graph structure and node features. To that end, we train overlapping and non-overlapping models on synthetic data with different degrees of correlation between network structure and node features.

4.1 Synthetic featured networks

We use a Stochastic Block Model [53] to generate synthetic featured networks, each with $M$ communities of $n=10$ nodes, with intra-cluster connection probabilities of $0.25$, and with inter-cluster connection probabilities of $0.01$. Each node is initially assigned a colour which encodes its feature community; we shuffle the colours of a fraction $1-\alpha$ of the nodes, randomly sampled. This procedure maintains constant the overall count of each colour, and lets us control the correlation between the graph structure and node features by moving $\alpha$ from 0 (no correlation) to 1 (full correlation).

Node features are represented by a one-hot encoding of their colour (therefore, in all our scenarios, the node features have dimension $M=N/n$). However, since in this case all the nodes inside a community have exactly the same feature value, the model can have difficulties differentiating nodes from one another. We therefore add a small Gaussian noise ($\sigma=.1$) to make sure that nodes in the same community can be distinguished from one another.

Note that the feature matrix has less degrees of freedom than the adjacency matrix in this setup, a fact that will be reflected in the plots below. However, opting for this minimal generative model lets us avoid the parameter exploration of more complex schemes for feature generation, while still demonstrating the effectiveness of our model.

4.2 Comparison setup

To evaluate the efficiency of our model in terms of capturing meaningful correlations between network and features, we compare overlapping and non-overlapping models as follows. For a given maximum number of embedding dimensions $F_{max}$, the overlapping models keep constant the number of dimensions used for adjacency matrix reconstruction and the number of dimensions used for feature reconstruction, with the same amount allocated to each task: $F^{ov}_{\mathbf{A}}+F^{ov}_{\mathbf{AX}}=F^{ov}_{\mathbf{X}}+F^{ov}_{\mathbf{AX}}=\frac{1}{2}F_{max}$. However they vary the overlap $F^{ov}_{\mathbf{AX}}$ from 0 to $\frac{1}{2}F_{max}$ by steps of 2. Thus the total number of embedding dimensions $F$ varies from $F_{max}$ to $\frac{1}{2}F_{max}$, and as $F$ decreases, $F^{ov}_{\mathbf{AX}}$ increases. We call one such model $\mathcal{M}^{ov}_{F}$.

Now for a given overlapping model $\mathcal{M}^{ov}_{F}$, we define a reference model $\mathcal{M}^{ref}_{F}$, which has the same total number of embedding dimensions, but without overlap: $F^{ref}_{\mathbf{AX}}=0$, and $F^{ref}_{\mathbf{A}}=F^{ref}_{\mathbf{X}}=\frac{1}{2}F$ (explaining why we vary $F$ with steps of 2). Note that while the reference model has the same information bottleneck as the overlapping model, it has less trainable parameters in the decoder, since $F^{ref}_{\mathbf{A}}+F^{ref}_{\mathbf{AX}}=F^{ref}_{\mathbf{X}}+F^{ref}_{\mathbf{AX}}=\frac{1}{2}F$ will decrease as $F$ decreases. Nevertheless, this will not be a problem for our measures, since we will be mainly looking at the behaviour of a given model for different values of $\alpha$ (i.e. the feature-network correlation parameter).

For our calculations (if not noted otherwise) we use synthetic networks of $N=1000$ nodes (i.e. 100 clusters), and set the maximum embedding dimensions $F_{max}$ to 20. For all models, we set the intermediate layer in the encoder and the two intermediate layers in the decoder to an output dimension of $50$, and the internal number of samples for loss estimation at $K=5$. We train our models for 1000 epochs using the Adam optimiser [54] with a learning rate of $0.01$ (following [30]), after initialising weights following [55]. For each combination of $F$ and $\alpha$, the training of the overlapping and reference models is repeated $20$ times on independent featured networks.

Since the size of our synthetic data is constant, and we average training results over independently sampled data sets, we can meaningfully compare the averaged training losses of models with different parameters. We therefore take the average best training loss of a model to be our main measure, indicating the capacity to reconstruct an input data set for a given information bottleneck and embedding overlap.

4.3 Advantages of overlap

4.3.1 Absolute loss values

Figure 2 shows the variation of the best training loss (total loss, adjacency reconstruction loss, and feature reconstruction loss) for both overlapping and reference models, with $\alpha$ ranging from 0 to 1 and $F$ decreasing from 20 to 10 by steps of 2. One curve in these plots represents the variation in losses of a model with fixed $F$ for data sets with increasing correlation between network and features; each point aggregates 20 independent trainings, used to bootstrap 95% confidence intervals.

We first see that all losses, whether for overlapping model or reference, decrease as we move from the uncorrelated scenario to the correlated scenario. This is true despite the fact that the total loss is dominated by the adjacency reconstruction loss, as feature reconstruction is an easier task overall. Second, recall that the decoder in a reference model has less parameters than its corresponding overlapping model of the same $F$ dimensions (except for zero overlap), such that the reference is less powerful and produces higher training losses. The absolute values of the losses for overlap and reference models are therefore not directly comparable. However, the changes in slopes are meaningful. Indeed, we note that the curve slopes are steeper for models with higher overlap (lower $F$) than for lower overlap (higher $F$), whereas they seem relatively independent for the reference models of different $F$. In other words, as we increase the overlap, our models seem to benefit more from an increase in network-feature correlation than what a reference model benefits.

4.3.2 Relative loss disadvantage

In order to assess this trend more reliably, we examine losses relative to the maximum embedding models. Figure 3 plots the loss disadvantage that overlap and reference models have compared to their corresponding model with $F=F_{max}$, that is, $\frac{\mathcal{L}_{\mathcal{M}_{F}}-\mathcal{L}_{\mathcal{M}_{F_{max}}}}{\mathcal{L}_{\mathcal{M}_{F_{max}}}}$. We call this the relative loss disadvantage of a model. In this plot, the height of a curve thus represents the magnitude of the decrease in performance of a model $\mathcal{M}^{ov|ref}_{F}$ relative to the model with maximum embedding size, $\mathcal{M}^{ov|ref}_{F_{max}}$. Note that for both the overlap model and the reference model, moving along one of the curves does not change the number of trainable parameters in the model.

As the correlation between network and features increases, we see that the relative loss disadvantage decreases in overlap models, and that the effect is stronger for higher overlaps. In other words, when the network and features are correlated, the overlap captures this joint information and compensates for the lower total number of dimensions (compared to $\mathcal{M}^{ov|ref}_{F_{max}}$): the model achieves a better performance than when network and features are more independent. Strikingly, for the reference model these curves are flat, thus indicating no variation in relative loss disadvantage with varying network-feature correlations in these cases. This confirms that the new measure successfully controls for the baseline decrease of absolute loss values when the network-features correlation increases, as observed in Figure 2. Our architecture is therefore capable of capturing and taking advantage of some of the correlation by leveraging the overlap dimensions of the embeddings.

Finally note that for high overlaps, the feature reconstruction loss value actually increases a little when $\alpha$ grows. The behaviour is consistent with the fact that the total loss is dominated by the adjacency matrix loss (the hardest task). In this case it seems that the total loss is improved more by exploiting the gain of optimising for adjacency matrix reconstruction, and paying the small cost of a lesser feature reconstruction, than decreasing both adjacency matrix and feature losses together. If wanted, this strategy could be controlled using a gradient normalisation scheme such as [52].

4.4 Standard benchmarks

Finally we compare the performance of our architecture to other well-known embedding methods, namely spectral clustering (SC) [56], DeepWalk (DW) [41], the vanilla non-variational and variational Graph Auto-Encoders (GAE and VGAE) [30], and GraphSAGE [26] which we look at in more detail. We do so on two tasks: (i) the link prediction task introduced by [30] and (ii) a node classification task, both on the Cora, CiteSeer and PubMed datasets, which are regularly used as citation network benchmarks in the literature [57, 58]. Note that neither SC nor DW support feature information as an input.

The Cora and CiteSeer datasets are citation networks made of respectively 2708 and 3312 machine learning articles, each assigned to a small number of document classes (7 for Cora, 6 for CiteSeer), with a bag-of-words feature vector for each article (respectively 1433 and 3703 words). The PubMed network is made of 19717 diabetes-related articles from the PubMed database, each assigned to one of three classes, with article feature vectors containing term frequency-inverse document frequency (TF/IDF) scores for 500 words.

4.4.1 Link prediction

The link prediction task consists in training a model on a version of the datasets where part of the edges has been removed, while node features are left intact. A test set is formed by randomly sampling 15% of the edges combined with the same number of random disconnected pairs (non-edges). Subsequently the model is trained on the remaining dataset where 15% of the real edges are missing.

We pick hyperparameters such that the restriction of our model to VGAE would match the hyperparameters used by [30]. That is a 32-dimensions intermediate layer in the encoder and the two intermediate layers in the decoder, and 16 embedding dimensions for each reconstruction task ($F_{\mathbf{A}}+F_{\mathbf{AX}}=F_{\mathbf{X}}+F_{\mathbf{AX}}=16$). We call the zero-overlap and the full-overlap versions of this model AN2VEC-0 and AN2VEC-16 respectively. In addition, we test a variant of these models with a shallow adjacency matrix decoder, consisting of a direct inner product between node embeddings, while keeping the two dense layers for feature decoding. Formally: $A_{ij}|\bm{\xi}_{i},\bm{\xi}_{j}\sim\operatorname{Ber}(\bm{\xi}^{T}_{i}\bm{\xi}_{j})$. This modified overlapping architecture can be seen as simply adding the feature decoding and embedding overlap mechanics to the vanilla VGAE. Consistently, we call the zero-overlap and full-overlap versions AN2VEC-S-0 and AN2VEC-S-16.

We follow the test procedure laid out by [30]: we train for $200$ epochs using the Adam optimiser [54] with a learning rate of $.01$, initialise weights following [55], and repeat each condition $10$ times. The $\bm{\mu}$ parameter of each node’s embedding is then used for link prediction (i.e. the parameter is put through the decoder directly without sampling), for which we report area under the ROC curve and average precision scores in Table 1.999Note that in [30], the training set is also 85% of the full dataset, and test and validation sets are formed with the remaining edges, respectively 10% and 5% (and the same amount of non-edges). Here, since we use the same hyperparameters as [30] we do not need a validation set. We therefore chose to use the full 15% remaining edges (with added non-edges) as a test set, as explained above.

We argue that AN2VEC-0 and AN2VEC-16 should have somewhat poorer performance than VGAE. These models are required to reconstruct an additional output, which is not directly used to the link prediction task at hand. First results confirmed our intuition. However, we found that the shallow decoder models AN2VEC-S-0 and AN2VEC-S-16 perform consistently better than the vanilla VGAE for Cora and CiteSeer while their deep counterparts (AN2VEC-0 and AN2VEC-16) outperforms VGAE for all datasets. As neither AN2VEC-0 nor AN2VEC-16 exhibited over-fitting, this behaviour is surprising and calls for further explorations which are beyond the scope of this paper (in particular, this may be specific to the link prediction task). Nonetheless, the higher performance of AN2VEC-S-0 and AN2VEC-S-16 over the vanilla VGAE on Cora and CiteSeer confirms that including feature reconstruction in the constraints of node embeddings is capable of increasing link prediction performance when feature and structure are not independent (consistent with [33, 35, 34]). An illustration of the embeddings produced by AN2VEC-S-16 on Cora is shown in Figure 4.

On the other hand, performance of AN2VEC-S-0 on PubMed is comparable with GAE and VGAE, while AN2VEC-S-16 has slightly lower performance. The fact that lower overlap models perform better on this dataset indicates that features and structure are less congruent here than in Cora or CiteSeer (again consistent with the comparisons found in [34]). Despite this, an advantage of the embeddings produced by the AN2VEC-S-16 model is that they encode both the network structure and the node features, and can therefore be used for downstream tasks involving both types of information.

We further explore the behaviour of the model for different sizes of the training set, ranging from 10% to 90% of the edges in each dataset (reducing the training set accordingly), and compare the behaviour of AN2VEC to GraphSAGE. To make the comparison meaningful we train two variants of the two-layer GraphSAGE model with mean aggregators and no bias vectors: one with an intermediate layer of 32 dimensions and an embedding layer of 16 dimensions (roughly equivalent in dimensions to the full overlap AN2VEC models), the second with an intermediate layer of 64 dimensions and an embedding layer of 32 dimensions (roughly equivalent to no overlap in AN2VEC). Both layers use neighbourhood sampling, 10 neighbours for the first layer and 5 for the second. Similarly to the shallow AN2VEC decoder, each pair of node embeddings is reduced by inner product and a sigmoid activation, yielding a scalar prediction between 0 and 1 for each possible edge. The model is trained on minibatches of 50 edges and non-edges (edges generated with random walks of length 5), learning rate 0.001, and 4 total epochs. Note that on Cora, one epoch represents about 542 minibatches,101010One epoch is 2708 nodes $\times$ 5 edges per node $\times$ 2 (for non-edges) = 27080 training edges or non-edges; divided by 50, this makes 541.6 minibatches per epoch. such that 4 epochs represent about 2166 gradient updates; thus with a learning rate of 0.001, we remain comparable to the 200 full batches with learning rate 0.01 used to train AN2VEC.

Figure 5 plots the AUC produced by AN2VEC and GraphSAGE for different training set sizes and different embedding sizes (and overlaps, for AN2VEC), for each dataset. As expected, the performance of both models decreases as the size of the test set increases, though less so for AN2VEC. For Cora and CiteSeer, similarly to Table 1, higher overlaps and a shallow decoder in AN2VEC give better performance. Notably, the shallow decoder version of AN2VEC with full overlap is still around .75 for a test size of 90%, whereas both GraphSAGE variants are well below .65. For PubMed, as in Table 1, the behaviour is different to the first two datasets, as overlaps 0 and 16 yield the best results. As for Cora and CiteSeer, the approach taken by AN2VEC gives good results: with a test size of 90%, all AN2VEC deep decoder variants are still above .75 (and shallow decoders above .70), whereas both GraphSAGE variants are below .50.

4.4.2 Node classification

Since the embeddings produced also to encode feature information, we then evaluate the model’s performance on a node classification task. Here the models are trained on a version of the dataset where a portion of the nodes (randomly selected) have been removed; next, a logistic classifier111111Using Scikit-learn’s [59] interface to the liblinear library, with one-vs-rest classes. is trained on the embeddings to classify training nodes into their classes; finally, embeddings are produced for the removed nodes, for which we show the F1 scores of the classifier.

Figure 6 shows the results for AN2VEC and GraphSAGE on all datasets. The scale of the reduction in performance as the test size increases is similar for both models (and similar to the behaviour for link prediction), though overlap and shallow versus deep decoding seem to have less effect. Still, the deep decoder is less affected by the change in test size than the shallow decoder; and contrary to the link prediction case, the 0 overlap models perform best (on all datasets). Overall, the performance levels of GraphSAGE and AN2VEC on this task are quite similar, with slightly better results of AN2VEC on Cora, slightly stronger performance for GraphSAGE on CiteSeer, and mixed behaviour on PubMed (AN2VEC is better for small test sizes and worse for large test sizes).

4.4.3 Variable embedding size

We also explore the behaviour of AN2VEC for different embedding sizes. We train models with $F_{\mathbf{A}}=F_{\mathbf{X}}\in\{8,16,24,32\}$ and overlaps 0, 8, 16, 24, 32 (whenever there are enough dimensions to do so), with variable test size. Figure 7 shows the AUC scores for link prediction, and Figure 8 shows the F1-micro scores for node classification, both on CiteSeer (the behaviour is similar on Cora, though less salient). For link prediction, beyond confirming trends already observed previously, we see that models with less total embedding dimensions perform slightly better than models with more total dimensions. More interestingly, all models seem to reach a plateau at overlap 8, and then exhibit a slightly fluctuating behaviour as overlap continues to increase (in models that have enough dimensions to do so). This is valid for all test sizes, and suggests (i) that at most 8 dimensions are necessary to capture the commonalities between network and features in CiteSeer, and (ii) that having more dimensions to capture either shared or non-shared information is not necessarily useful. In other words, 8 overlapping dimensions seem to capture most of what can be captured by AN2VEC on the CiteSeer dataset, and further increase in dimensions (either overlapping or not) would capture redundant information.

Node classification, on the other hand, does not exhibit any consistent behaviour beyond the reduction in performance as the test size increases. Models with less total dimensions seems to perform slightly better at 0 overlap (though this behaviour is reversed on Cora), but neither the ordering of models by total dimensions nor the effect of increasing overlap are consistent across all conditions. This suggests, similarly to Figure 6a, that overlap is less relevant to this particular node classification scheme than it is to link prediction.

4.4.4 Memory usage and time complexity

Finally, we evaluate the resources used by our implementation of the method in terms of training time and memory usage. We use AN2VEC with 100-dimensions intermediate layers in the encoder and the (deep) decoder with 16 embedding dimensions for each reconstruction task ($F_{\mathbf{A}}+F_{\mathbf{AX}}=F_{\mathbf{X}}+F_{\mathbf{AX}}=16$), and overlap $F_{\mathbf{AX}}\in\{0,8,16\}$. We train that model on synthetic networks generated as in section Synthetic featured networks (setting $\alpha=0.8$, and without adding any other noise on the features), with $M\in\{50,100,200,500,1000,2000,5000\}$ communities of size $n=10$ nodes.

Only CPUs were used for the computations, running on a 4 $\times$ Intel Xeon CPU E7-8890 v4 server with 1.5 TB of memory. Using 8 parallel threads for training,121212There are actually two levels of threading: the number of threads used in our code for computing losses, and the number of threads used by the BLAS routines for matrix multiplication. We set both to 8, and since both computations alternate this leads to an effective 8 compute threads, with some fluctuations at times. we record the peak memory usage,131313Using the top utility program. training time, and full job time141414As reported by our scripts and by GNU Parallel. for each network size, averaged over the three overlap levels. Results are shown in Figure 9. Note that in a production setting, multiplying the number of threads by $n$ will divide compute times by nearly $n$, since the process is aggressively parallelised. A further reduced memory footprint can also be achieved by using sparse encoding for all matrices.

5 Conclusions

In this work, we proposed an attributed network embedding method based on the combination of Graph Convolutional Networks and Variational Autoencoders. Beyond the novelty of this architecture, it is able to consider jointly network information and node attributes for the embedding of nodes. We further introduced a control parameter able to regulate the amount of information allocated to the reconstruction of the network, the features, or both. In doing so, we showed how shallow versions of the proposed model outperform the corresponding non-interacting reference embeddings on given benchmarks, and demonstrated how this overlap parameter consistently captures joint network-feature information when they are correlated.

Our method opens several new lines of research and applications in fields where attributed networks are relevant. As an example one can take a social network with the task of predicting future social ties, or reconstruct existing but invisible social ties. Solutions to this problem can rely on network similarities in terms of overlapping sets of common friends, or on feature similarities in terms of common interest, professional or cultural background, and so on. While considering these types of information separately would provide us with a clear performance gain in the prediction, these similarities are not independent. For example, common friends may belong to the same community. By exploiting these dependencies our method can provide us with an edge in terms of predictive performance and could indicate which similarities, structural, feature-related, or both, better explain why a social tie exists at all. Another setting where we believe our framework might yield noteworthy insights is when applied to the prediction of side effects of drug pairs (polypharmacy). This problem has recently been approached by Zitnik et al. [60] by extending GraphSAGE for multirelational link prediction in multimodal networks. In doing so, the authors were able to generate multiple novel candidates of drug pairs susceptible to induce side effects when taken together. Beyond using drug feature vectors to generate polypharmacy edge probabilities, our overlapping encoder units would enable a detailed view on how these side effects occur due to confounding effects of particular drug attributes. It would pinpoint the feature pairs that interacting drugs might share (or not), further informing the drug design process. Furthermore, we expect that our method will help yield a deeper understanding between node features and structure, to better predict network evolution and ongoing dynamical phenomena. In particular, it should help to identify nodes with special roles in the network by clarifying whether their importance has structural or feature origin.

In this paper our aim was to ground our method and demonstrate its usefulness on small but controllable featured networks. Its evaluation on more complex synthetic datasets, in particular with richer generation schemes, as well as its application to larger real datasets, are therefore our immediate goals in the future.

List of abbreviations

VAE

Variational Autoencoder

GCN

Graph Convolutional Network

GCN-VAE

Graph Convolutional Variational Autoencoder

GC

Graph Convolutional layer

MLP

Multi-layer perceptron

Ber

Bernoulli random variable

ReLU

Rectified Linuar Unit

KL

Kullback-Leibler divergence

AN2VEC

Attributed Node to Vector model:

• •

AN2VEC-0: zero overlap model

• •

AN2VEC-16: 16-dimensions overlap model

• •

AN2VEC-S-0: zero overlap model with shallow adjacency decoder

• •

AN2VEC-S-16: 16-dimensions overlap model with shallow adjacency decoder

SC

Spectral Clustering

DW

DeepWalk embedding model

DANE

Deep Attributed Network Embedding

GAE

Graph Autoencoder

VGAE

Variational Graph Autoencoder

TF/IDF

Term-frequency-inverse-document-frequency

AUC

Area under the ROC curve

AP

Average precision

ROC

Receiver operating characteristic

Availability of data and material

The synthetic datasets generated for this work are stochastically created by our implementation, available at github.com/ixxi-dante/an2vec.

The datasets used for standard benchmarking (Cora, CiteSeer, and PubMed) are available at linqs.soe.ucsc.edu/data.

Our implementation of AN2VEC is made using the Julia programming language, and particularly making heavy use of Flux [61]. Parallel computations were run using GNU Parallel [62]. Finally, we used StellarGraph [63] for the GraphSAGE implementation.

Competing interests

The authors declare that they have no competing interests.

Funding

This project was supported by the LIAISON Inria-PRE project, the SoSweet ANR project (ANR-15-CE38-0011), and the ACADEMICS project financed by IDEX LYON.

Author’s contributions

MK, JLA and SL participated equally in designing and developing the project, and in writing the paper. SL implemented the model and experiments. SL and JLA developed and implemented the analysis of the results.

Acknowledgements

We thank E. Fleury, J-Ph. Magué, D. Seddah, and E. De La Clergerie for constructive discussions and for their advice on data management and analysis. Some computations for this work were made using the experimental GPU platform at the Centre Blaise Pascal of ENS Lyon, relying on the SIDUS infrastructure provided by E. Quemener.