# A Hierarchical Block Distance Model for Ultra Low-Dimensional Graph
Representations

Nikolaos Nakis
Section for Cognitive Systems, Technical University of Denmark
and
Abdulkadir Çelikkanat
Section for Cognitive Systems, Technical University of Denmark
and
Sune Lehmann Jørgensen
Section for Cognitive Systems, Technical University of Denmark
and
Morten Mørup
Section for Cognitive Systems, Technical University of Denmark The authors
gratefully acknowledge the Independent Research Fund Denmark for supporting
this work [grant number: 0136-00315B]

###### Abstract

Graph Representation Learning (GRL) has become central for characterizing
structures of complex networks and performing tasks such as link prediction,
node classification, network reconstruction, and community detection. Whereas
numerous generative GRL models have been proposed, many approaches have
prohibitive computational requirements hampering large-scale network analysis,
fewer are able to explicitly account for structure emerging at multiple
scales, and only a few explicitly respect important network properties such as
homophily and transitivity. This paper proposes a novel scalable graph
representation learning method named the Hierarchical Block Distance Model
(HBDM). The HBDM imposes a multiscale block structure akin to stochastic block
modeling (SBM) and accounts for homophily and transitivity by accurately
approximating the latent distance model (LDM) throughout the inferred
hierarchy. The HBDM naturally accommodates unipartite, directed, and bipartite
networks whereas the hierarchy is designed to ensure linearithmic time and
space complexity enabling the analysis of very large-scale networks. We
evaluate the performance of the HBDM on massive networks consisting of
millions of nodes. Importantly, we find that the proposed HBDM framework
significantly outperforms recent scalable approaches in all considered
downstream tasks. Surprisingly, we observe superior performance even imposing
ultra-low two-dimensional embeddings facilitating accurate direct and
hierarchical-aware network visualization and interpretation.

Keywords: graph representation learning, link prediction, low-dimensional
embeddings, latent-space, networks, large-scale

##  1 Introduction

Networks naturally arise in a plethora of scientific areas to model
interactions between entities from physics to sociology and biology, with many
instances such as collaboration, protein-protein interaction, and brain
connectivity networks (Newman, 2003) to mention but a few. In recent years,
Graph Representation Learning (GRL) approaches have attracted great interest
due to their outstanding performance compared to classical techniques for
arduous problems such as link prediction (Liben-Nowell and Kleinberg, 2003),
node classification (Getoor and Taskar, 2007; Grover and Leskovec, 2016), and
community detection (Fortunato, 2010).

Numerous GRL methods have been proposed, see also (Zhang et al., 2020) for a
survey. The leading initial works are the random walk-based methods (Perozzi
et al., 2014; Grover and Leskovec, 2016; Tang et al., 2015; Çelikkanat and
Malliaros, 2020; Nguyen and Malliaros, 2018), leveraging the Skip-Gra m
algorithm (Mikolov et al., 2013) to learn the node representations. Matrix
factorization-based algorithms (Cao et al., 2015; Zhang et al., 2020) have
also become prominent, extracting the embedding vectors by decomposing a
designed feature matrix. Furthermore, neural network models (Zhang et al.,
2020; Hamilton et al., 2017) have been proposed for graph-structured data,
returning outstanding performance by combining node attributes and network
structure when learning embeddings. Although recent studies (Zhang et al.,
2019) aim to alleviate the computational burden of these algorithms through
matrix sparsification tools (Qiu et al., 2019), hierarchical representations
(Bhowmick et al., 2020; Chen et al., 2018), or by fast hashing schemes (Zhang
et al., 2018), the performance of the methods in downstream tasks
significantly drops, and they require larger embedding sizes to compensate for
this loss.

Latent Space Models (LSMs) for the representation of graphs have been quite
popular over the past years, especially for social networks analysis (Sosa and
Buitrago, 2020; Kim et al., 2017). LSMs utilize the generalized linear model
framework to obtain informative latent node embeddings while preserving
network characteristics. The choice of latent effects in modeling the link
probabilities between the nodes leads to different expressive capabilities
characterizing network structure. In particular, in the Latent Distance Model
(LDM) (Hoff et al., 2002) nodes are placed closer in the latent space if they
are similar or vice-versa. LDM obeys the triangle inequality and thus
naturally represents transitivity (”a friend of a friend is a friend”) and
network homophily (a tendency where similar nodes are more likely to connect
to each other than dissimilar ones). The approach has been extended to
bipartite networks in Friel et al. (2016) by introducing mode-specific
embedding vectors and community detection by endowing the LDM with a Gaussian
Mixture Model prior promoting cluster structures in the latent space forming
the latent position clustering model (LPCM)Handcock et al. (2007a); Krivitsky
et al. (2009). However, the LDM is unable to account for possible stochastic
equivalence as defined by the stochastic block models Stochastic Blockmodels
(Holland et al., 1983; Nowicki and Snijders, 2001), i.e. (”groups of nodes
defined by shared intra- and inter-group relationships”) defining non
positive-semi-definite latent representations. The LSMs were advanced to
characterize such stochastic equivalence by imposing an Eigenmodel admitting
negative eigenvalues (Hoff, 2007). These latent space methods are attractive
due to their simplicity, as they define well-structured inference problems and
are characterized by high explanatory power. The time and space complexities
are their main drawbacks as the likelihood function scales by the number of
node pairs (i.e., quadratically in the number of nodes for a unipartite graph)
typically addressed using subsampling strategies (Raftery et al., 2012).

Many real-world networks are composed of structures emerging at multiple
scales which can be expressed using hierarchical representations (Ravasz and
Barabási, 2003). Several methods have thus been advanced to such hierarchical
representations including stochastic block model approaches (Clauset et al.,
2008; Roy et al., 2007; Roy and Teh, 2008; Herlau et al., 2012, 2013; Peixoto,
2014) as well as agglomerative (Blondel et al., 2008; Ahn et al., 2010;
Blundell and Teh, 2013) and recursive partitioning (Li et al., 2020)
procedures relying on various measures of similarity. Importantly, learning
node representations characterizing structure at multiple scales of the
network can facilitate network visualization and the understanding of the
inner dynamics of networks. Hierarchical representation of bipartite networks
is of special interest due to the fact that most unipartite hierarchical
clustering algorithms do not generalize to the bipartite case beyond
clustering each mode separately or transforming the bipartite network into a
unipartite representation. In the work of Wieling and Nerbonne (2010), the
authors used the spectral partitioning algorithm of (Dhillon, 2001) and then
applied k-means on the spectral space to get initial bi-clusters which were
followed by divisive bi-splits to create a dendrogram. In this case, the
spectral embedding space was not constructed to reflect explicitly the
clustering criterion. In addition to divisive procedures, agglomerative
clustering has also been proposed for bipartite networks. In the work of Cai
and Liu (2020) a multi-objective function was designed and combined with
classical community construction algorithms. One limitation here is that the
network should be transformed to a unipartite structure.

Despite the many advantages of hierarchical structures and block models, one
major limitation remains to accurately account for homophily (Hoff, 2007)
which is a key characteristic of social networks. More specifically, block
models have been extended to explicitly impose a community structure Rosvall
and Bergstrom (2007); Mørup and Schmidt (2012) but notable this only provides
within cluster homogeneity and thus homophily-like properties for the
community relative to the other communities but not a hierarchy complying with
such a structure. Whereas LPCM accounts for homophily it does not account for
hierarchical structures and cluster structures are not strictly imposed beyond
a prior promoting the latent positions to form groups.

In this work, we propose a novel node representation learning approach, the
Hierarchical Block Distance Model (HBDM)111Source code: Hierarchical Block
Distance Model., as a reconciliation between hierarchical block structures of
different scales and network properties such as homophily and transitivity. In
particular, we propose a framework combining embedding and hierarchical
characterization for graph representation learning. Importantly, we design a
hierarchical structure that respects a linearithmic total time and space
complexity, in terms of the number of nodes (i.e., $\mathcal{O}(N\log N)$),
and at the same time provides an accurate interpretable representation of
structure at different scales. Using the HBDM, we embed moderate-sized and
large-scale networks containing more than a million nodes and establish the
performance of our model in terms of link prediction and node classification
to existing prominent graph embedding approaches. We further highlight how the
inferred hierarchical organization can facilitate accurate visualization of
network structure even when using only $D=2$ dimensional representations
providing favorable performance in all the considered GRL tasks; link-
prediction, node classification, and network reconstruction. Additionally, we
show how our proposed framework extends the hierarchical multi-resolution
structure to bipartite networks and provides the characterization of
communities at multiple scales.

##  2 The Hierarchical Block Distance Model

We firstly concentrate our study on undirected networks and later generalize
our approach to bipartite graphs. We will firstly provide the necessary
definitions required throughout the paper. Let $\mathcal{G}=(V,E)$ be a graph
where $N:=\left|V\right|$ is the number of nodes and $Y_{N\times
N}=\left(y_{i,j}\right)\in\\{0,1\\}^{N\times N}$ be the adjacency matrix of
the graph such that $y_{i,j}=1$ if the pair $(i,j)\in E$ otherwise it is equal
$0$, for all $1\leq i<j\leq N$. We denote the latent representations of nodes
by $\mathbf{Z}=(z_{i,d})\in\mathbb{R}^{N\times D}$ where each row vector,
$\mathbf{z}_{i}\in\mathbb{R}^{D}$, indicates the corresponding embedding of
node $i\in\mathcal{V}$ in a $D$-dimensional space.

The most natural choice for modeling homophily and transitivity can be found
in the Latent Space Model (LSM) which defines an $\mathbb{R}^{D}$-dimensional
latent space in which every node of the graph is characterized through the
unobserved but informative node-specific variables
$\\{\mathbf{z}_{i}\in\mathbb{R}^{D}\\}$. These variables are considered
sufficient to describe and explain the underlying relationships between the
nodes of the network. The probability of an edge occurring is considered
conditionally independent given the unobserved latent positions. Consequently,
the total probability distribution of the network can be written as:

$P(Y|\bm{Z},\bm{\theta})=\prod_{i<j}^{N}p(y_{i,j}|\mathbf{z}_{i},\mathbf{z}_{j},\bm{\theta}_{i,j}),$
(1)

where $\bm{\theta}$ denotes any potential additional parameters, such as
covariate regressors. A popular and convenient parameterization of Equation
(1) for binary data is through the logistic regression model (Hoff et al.,
2002; Handcock et al., 2007b; Krivitsky et al., 2009; Hoff, 2005). In
contrast, we adopt the Poisson regression model similar to Hoff (2005) under a
generalized linear model framework for the LSM. The use of a Poisson
likelihood for modelling binary relationships in a network does not decrease
the predictive performance nor the ability of the model to detect the network
structure, as shown in Wind and Mørup (2012). It also generalizes the analysis
to integer weighted graphs. In addition, the exchange of the logit to a log
link function when transitioning from a Bernoulli to a Poisson model defines
nice decoupling properties over the predictor variables in the likelihood
(Karrer and Newman, 2011; Herlau et al., 2014). Utilizing the Poisson Latent
Distance Model (LDM) of the LSM family framework, the rate of an occurring
edge depends on a distance metric between the latent positions of the two
nodes. In our formulation, we consider the LDM Poisson rate with node-specific
biases or random-effects (Hoff, 2005; Krivitsky et al., 2009) such that the
expression for the Poisson rate becomes:

$\lambda_{ij}=\exp\big{(}\gamma_{i}+\gamma_{j}-d(\mathbf{z}_{i},\mathbf{z}_{j})\big{)},$
(2)

where $\gamma_{i}\in\mathbb{R}$ denotes the node-specific random-effects and
$d_{ij}(\cdot,\cdot)$ denotes any distance metric obeying the triangle
inequality $\big{\\{}d_{ij}\leq d_{ik}+d_{kj},\>\forall(i,j,k)\in
V^{3}\big{\\}}$. Considering variables $\\{\mathbf{z}_{i}\\}_{i\in V}$ as the
latent characteristics, Equation (2) shows that similar nodes will be placed
closer in the latent space, yielding a high probability of an occurring edge
and thus modeling homophily and satisfies network transitivity and reciprocity
through the triangle inequality whereas the node specific bias can account for
degree heterogeneity. The conventional LDM rate utilizing a global bias,
$\gamma^{g}$, corresponds to the special case in which
$\gamma_{i}=\gamma_{j}=0.5\gamma^{g}$. As in Hoff et al. (2002), we presently
adopt the Euclidean distance as the choice for the distance metric
$d_{ij}(\cdot,\cdot)$.

###  2.1 Designing A Linearithmic Complexity

Our goal is to design a Hierarchical Block Model preserving homophily and
transitivity properties with a total complexity allowing for the analysis of
large-scale networks. Our HBDM, defines the rate of a link between each
network dyad $\\{i,j\\}\in V\times V$ based on the Euclidean distance, as
shown in Equation (2). Therefore, we can define a block-alike hierarchical
structure by a divisive clustering procedure over the latent variables in the
Euclidean space. The total optimization cost of such a model though should
have a linearithmic upper bound complexity to make large-scale analysis
feasible. Introducing a number of clusters $K$ equal to the number of nodes
$N$ in the HBDM, leads to the same log-likelihood as of the standard LDM,
defining a sum over each ordered pair of the network, as:

$\log
P(Y|\bm{\Lambda})=\sum_{i<j}\Big{(}y_{ij}\log(\lambda_{ij})-\lambda_{ij}\Big{)}=\sum_{i<j:y_{ij}=1}\log(\lambda_{ij})\;-\;\sum_{i<j}\lambda_{ij}\>,$
(3)

For brevity, we presently ignore the linear scaling by dimensionality $D$ of
the above log-likelihood function. Notably, the link contribution
$\sum_{y_{i,j}=1}\log(\lambda_{i,j})$ is responsible for positioning ”similar”
nodes closer in the latent space, expressing the desired homophily.

In addition, large networks are highly sparse (Barabási and Pósfai, 2016) with
the number of edges for very sparse networks being proportional to the number
of nodes in the network. As a result, the computation of the link contribution
is relatively cheap scaling linearithmic or sub-linearithmic (as shown in
supplementary). Most importantly, the link term removes rotational ambiguity
between the different blocks of the hierarchy (as discussed later). For these
three reasons, no block structure is imposed on the calculation of the link
contribution. The second term acts as the repelling force for dissimilar nodes
and requires the computation of all node pairs scaling as $\mathcal{O}(N^{2})$
making the evaluation of the above likelihood infeasible for large networks.
By enforcing a block structure, i.e., akin to stochastic block models White et
al. (1976); Holland et al. (1983); Nowicki and Snijders (2001), when grouping
the nodes into $K$ clusters we define the rate between block $k$ and
$k^{\prime}$ in terms of their distance between centroids. A simple block
structure without a hierarchy would lead to the following non-link expression:

$\sum_{i<j}\lambda_{ij}\\!\approx\\!\sum_{k=1}^{K}\\!\Bigg{(}\\!\\!\sum_{\begin{subarray}{c}i<j\\\
i,j\in
C_{k}\end{subarray}}\\!\\!\exp{\\!\big{(}\gamma_{i}\\!+\\!\gamma_{j}\\!-\\!||\mathbf{z}_{i}-\mathbf{z}_{j}||_{2}\big{)}}\\!+\\!\sum_{k^{{}^{\prime}}>k}^{K}\sum_{i\in
C_{k}}\sum_{j\in
C_{k^{{}^{\prime}}}}\exp{\big{(}\gamma_{i}+\gamma_{j}-||\bm{\mu}_{k}-\bm{\mu}_{k^{\prime}}||_{2}\big{)}}\Bigg{)},$
(4)

where $\bm{\mu}_{k}$ denotes the $k$’th cluster centroid of the set
$\bm{C}=\\{C_{1},\dots,C_{K}\\}$, and has absorbed the dependency over the
variables $\bm{Z}\in\mathbb{R}^{N\times D}$. More specifically, the cluster
centroids $\bm{\mu}_{k}$ are implicit parameters defined as a function over
the latent variables, as we will show later. Overall, considering the
principle that connected and homophilic nodes will be placed closer in the
latent space, this expression generalizes this principle by introducing a
clustering procedure that obeys ”cluster-homophily” and ”cluster-transitivity”
over the latent space. More specifically, we can assume that closely related
nodes will be positioned in the same cluster while related or interconnected
clusters will also be positioned close in the latent space, providing an
accurate block structure schema. As opposed to the LPCM where clustering
structures are imposed through a prior, the above formulation strictly defines
the clustering structure as shared overall proximity between blocks as defined
by the distances between centroids of the formed groups.

(a) (i)

(b) (ii)

Figure 1: Schematic representation of the distance matrix calculation for a
hierarchical structure of the tree of height $L=3$ and for the number of
observations $N=64$.

####  2.1.1 A Hierarchical Representation

In order to obtain the desired hierarchical representation, we define
hierarchical clustering via a divisive procedure. In detail, we organize the
embedded clusters into a hierarchy using a tree structure, defining a cluster
dendrogram. The root of the tree is a single cluster containing the total
amount of latent variable embeddings $\bm{Z}$. At every level of the tree, we
perform partitioning until we obtain leaf nodes containing equal or less than
the desired number of nodes, $N_{leaf}$. This number is chosen with respect to
our linearithmic complexity upper bound and set as $N_{\text{leaf}}=\log N$,
resulting in approximately $K=N/log(N)$ total clusters. The tree-nodes
belonging to a specific tree-level are considered the clusters for that
specific tree height. Every novel partition of a non-leaf node is performed
only on the set of points allocated to the parent tree-node (cluster). For
every level of the tree, we consider the pairwise distances of datapoints
belonging to different tree-nodes as the distance between the corresponding
cluster centroids, as illustrated by Figure 1 (ii). Based on these distances,
we calculate the likelihood contribution of the blocks and continue with
binary splits, down the tree, for the non-leaf tree-nodes. When all tree-nodes
are considered as leaves, we calculate analytically the inner cluster pairwise
distances for the corresponding likelihood contribution of analytical blocks,
as shown in the last part of Figure 1 (ii). The latter analytical calculation
comes at a linearithmic cost of
$\mathcal{O}(KN_{\text{leaf}}^{2})=\mathcal{O}(N\log N)$ while enforces the
homophily and transitivity properties of the model since for the most similar
nodes the HBDM behaves explicitly as the standard LDM.

We can thereby define a Hierarchical Block Distance Model with Random Effects
(HBDM-Re) as:

$\displaystyle\log P(Y|\mathbf{Z},\bm{\gamma})$
$\displaystyle=\sum_{\begin{subarray}{c}i<j\\\
y_{i,j}=1\end{subarray}}\Bigg{(}\gamma_{i}+\gamma_{j}-||\mathbf{z}_{i}-\mathbf{z}_{j}||_{2}\Bigg{)}-\sum_{k=1}^{K_{L}}\Bigg{(}\sum_{\begin{subarray}{c}i<j\\\
i,j\in
C^{(L)}_{k}\end{subarray}}\exp{(\gamma_{i}+\gamma_{j}-||\mathbf{z}_{i}-\mathbf{z}_{j}||_{2})}\Bigg{)}$
$\displaystyle-\sum_{l=1}^{L}\sum_{k=1}^{K_{l}}\sum_{k^{{}^{\prime}}>k}^{K_{l}}\exp{(-||\bm{\mu}^{(l)}_{k}-\bm{\mu}^{(l)}_{k^{\prime}}||_{2})}\Bigg{(}\sum_{i\in
C_{k}^{(l)}}\exp{(\gamma_{i})}\Bigg{)}\Bigg{(}\sum_{j\in
C_{k^{{}^{\prime}}}^{(l)}}\exp{(\gamma_{j})}\Bigg{)},$ (5)

where $l\in\\{1,\ldots,L\\}$ denotes the $l$’th dendrogram level, $k_{l}$ is
the index representing the cluster id for the different tree levels and
$\bm{\mu}_{k}^{(l)}$ the corresponding centroid. We also consider a
Hierarchical Block Distance Model (HBDM) without the random effects setting
$\gamma_{i}=0.5\gamma^{g}$. For a multifurcating tree splitting in $K$
clusters and having $N/log(N)$ terminal nodes (clusters), the number of
internal nodes are $\mathcal{O}(N/(K\log{N}))$ and each node needs to evaluate
$\mathcal{O}(K^{2})$ pairs providing an overall complexity of
$\mathcal{O}(NK/\log{N})$, thus $K\leq\log{N}^{2}$ to achieve
$\mathcal{O}(N\log{N})$ scaling (Epp, 2010).

####  2.1.2 Divisive partitioning using k-means with a Euclidean distance
metric

Whereas the likelihood in Equation (5) can be directly minimized by assigning
nodes to the clusters given by the tree structure, this evaluation for all $N$
nodes scales prohibitively as $\mathcal{O}(N^{2}/\log{N})$. To reduce this
scaling, we use a more efficient divisive partitioning procedure. For
efficiency instead of considering minimizing the full expression in Equation
(5) we consider minimizing the Euclidean norm
$||\bm{\mu}_{k_{l}}-\bm{\mu}_{k_{l}^{\prime}}||_{2}$. The divisive clustering
procedure thus relies on the following Euclidean norm objective

$J(\mathbf{r},\bm{\mu})=\sum_{i=1}^{N}\sum_{k=1}^{K}r_{ik}||\mathbf{z}_{i}-\bm{\mu}_{k}||_{2},$
(6)

where $k$ denotes the cluster id, $\mathbf{z}_{i}$ is the i’th data
observation, $r_{ik}$ the cluster responsibility/assignment, and
$\bm{\mu}_{k}$ the cluster centroid.

This objective function is unfortunately not accounted for by existing k-means
clustering algorithms relying on the squared Euclidean norm. We therefore
presently derive an optimization procedure for k-means clustering with
Euclidean norm utilizing the auxiliary function framework of Tsutsu and
Morikawa (2012) developed in the context of compressed sensing. We define an
auxiliary function for (6) as:

$J^{+}(\bm{\phi},\mathbf{r},\bm{\mu})=\sum_{i=1}^{N}\sum_{k=1}^{K}r_{ik}\Bigg{(}\frac{||\mathbf{z}_{i}-\bm{\mu}_{k}||_{2}^{2}}{2\phi_{ik}}+\frac{1}{2}\phi_{ik}\Bigg{)},$
(7)

where $\bm{\phi}$ are the auxiliary variables. Thereby, minimizing Equation
(7) with respect to $\bm{\phi}_{nk}$ yields
$\bm{\phi}_{ik}^{*}=||\mathbf{z}_{i}-\bm{\mu}_{k}||_{2}$ and by plugging
$\bm{\phi}_{ik}^{*}$ back to (6) we obtain
$J^{+}(\bm{\phi}^{*},\mathbf{r},\bm{\mu})=J(\mathbf{r},\bm{\mu})$ verifying
that (7) is indeed a valid auxiliary function for (6). The algorithm proceeds
by optimizing cluster centroids as $\bm{\mu}_{k}=\Big{(}\sum_{i\in
k}\frac{\mathbf{z_{i}}}{\phi_{ik}}/\sum_{i\in k}\frac{1}{\phi_{ik}}\Big{)}$
and assigning points to centroids as
$\operatorname*{arg\,min}_{\bm{C}}=\sum_{k=1}^{K}\sum_{\bm{z}\in
C_{k}}\Big{(}\frac{||\bm{z}-\bm{\mu}_{k}||_{2}^{2}}{2\bm{\phi}_{k}}+\frac{1}{2}\bm{\phi}_{k}\Big{)}$
upon which $\bm{\phi}_{k}$ is updated. The overall complexity of this
procedure is $\mathcal{O}(TKND)$ (Hartigan and Wong, 1979) where $T$ is the
number of iterations required to converge.

A simple approach to construct the tree structure would be to use the above
Euclidean k-means procedure to split the nodes into $K=N/\log(N)$ clusters and
construct the tree according to agglomeration as in hierarchical clustering.
Unfortunately, such a strategy is prohibitive. For that, in the coarser level
(first layer of the tree), we choose to split to the maximally allowed
clusters of $K=\log{N}$ allowing scaling of $\mathcal{O}(N\log{N})$. It would
be tempting to continue splitting into $\log{N}$ clusters, however, for a
balanced multifurcating tree with $N/\log{N}$ leaf clusters, it will result in
a height scaling as $\mathcal{O}(\log{N}/\log{\log{N}})$ and thus an overall
complexity of $\mathcal{O}(N\log^{2}(N)/\log{\log{N}})$ (Epp, 2010). Whereas a
balanced binary tree at all levels below the root results in a height scaling
as $\mathcal{O}(\log{N})$ providing an overall complexity when including the
linear scaling by dimensionality $D$ of $\mathcal{O}(DN\log{N})$ (as each
level of the tree defines $\mathcal{O}(DN)$ operations). Figure 1 (i),
illustrates the resulting tree222 For visualization purposes only, we show
equally sized clusters. for a small problem of $N=64$ nodes in which we first
split into 4 ($\approx\log(64)$) clusters and subsequently create binary
splits until each leaf cluster contains 4 ($\approx\log(64)$) nodes.

####  2.1.3 Expressing Homophily and Transitivity

(a) (i) Non-optimal rotation over leaf clusters. (b) (ii) Optimal rotation
over leaf clusters.

(c) (iii) Three latent block structures.

Figure 2: The clusters within the dashed circles denote the leaf block
structures. The red circles and blue rhombuses indicate the node embeddings
and the centroids, respectively. Gray lines represent the links and the dashed
lines the distance between the cluster centers.

A central component to preserve homophily and transitiviy of the HBDM is not
to approximate the link terms at the level of block as in (hierarchical) SBMs
but calculating analytically the link contribution of the log-likelihood
across the total hierarchy beyond the leaf/analytical blocks. In Figures 2 (i)
and (ii), two leaf clusters are illustrated connected with a link. Assume that
we only calculate the distance inside the blocks analytically and that both
the link and non-link contributions of pairs belonging to different clusters
are approximated based on their centroids’ distance. This essentially would
allow for any rotation of each cluster for all clusters in the hierarchy since
the inner-block distances (analytical), as well as the centroid distances,
would not change by such rotations, yielding exactly the same likelihood
(block-level rotational invariance). In that case, homophily would be violated
as, e.g., the distance between nodes $c$ and $d$ would not necessarily be
smaller than the distance with other inter-cluster pairs (ex: Figure 2 (i)),
showing that the rotation of the blocks substantially impacts the homophily
properties of the HBDM. Calculating the link contributions between the
different clusters analytically solves this ambiguity since the likelihood is
penalized higher when nodes $c$ and $d$ are positioned in a non-rotational-
aware way. The computational cost imposed by accouting for all the link terms
analytically is that the model complexity depends on the number of edges of
the network (a total block structure would strictly be linearithmic in
complexity). Nevertheless, we show empirically in the supplementary that the
number of network edges $E$ scales linearly with $N\log N$ and thus this
analytical term respects our complexity bound. In Figure 2 (iii), we present
clusters defining cases of block inter-connections of sparsely connected
blocks ({$C_{1},C_{3}$}, {$C_{1},C_{2}$}) and densely connected blocks
{$C_{2},C_{3}$}. Whereas the analytical inter-cluster links (blue lines) are
responsible for fixing the block rotation the inter cluster links also drives
the cluster level proximities of centroids ensuring cluster-homophily and
transitivity.

Pairwise distances in the HBDM stays invariant to rotation, reflection, and
translation of the latent space due to its LDM inheritance (Hoff et al., 2002)
these isometries can be resolved via a Singular-Value-Decomposition procedure
as provided in the supplementary. Whereas the analytical link term
calculations provide rotational awareness to the HBDM clusters, we continue by
investigating the conditions in which a continuous operations defining
infinitesimal rotations (with respect to the cluster centroid) is admissible
leaving the loss function of Equation (5) invariant to continuous rotations.
In Lemma 2.1 (proof given in the supplementary material), we start our
investigation of this problem by showing that blocks with a unique inter-
cluster link connection reduce the clusters’ degree of rotational freedom by
one.

######  Lemma 2.1.

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a graph and let $\mathcal{C}$
be a cluster with its centroid located at $\bm{\mu}\in\mathbb{R}^{D}$ having
an edge $(i,j)\in\mathcal{E}$ for some $i\in\mathcal{C}$ and
$j\in\mathcal{V}\backslash\mathcal{C}$ such that $\mathbf{z}_{i}\neq\bm{\mu}$.
If
$\mathbf{\tilde{z}}_{i}=\bm{\mu}+\mathbf{R}(\bm{\theta})(\mathbf{z}_{i}-\bm{\mu})$
such that $\mathbf{R}(\bm{\theta})$ is a rotation matrix acting on the
embeddings of nodes in cluster $\mathcal{C}$, then the maximum degree of
freedom of any infinitesimal $\lambda_{ij}$-invariant rotation is defined by
$\bm{\theta}\in\mathbb{R}^{D-2}$.

A direct consequence of Lemma 2.1 is that for a two-dimensional embedding,
there is no possible continuous rotation of a cluster having only one external
edge. Since there is a path from one node to all others in a connected graph,
every cluster must have at least one external link. For the general case of
blocks having multiple inter-cluster edges, rotations preserving the total sum
of pairwise distances among node embeddings are highly unlikely, as discussed
in the supplementary. Consequently, we can for connected networks expect
uniqueness of a (local) minima solutions with no continuous admissible
rotations leaving the HBDM loss function of Equation (5) invariant.

####  2.1.4 Extension to Bipartite Networks

Our proposed frameworks, HBDM and HBDM-Re generalize to both directed and
bipartite graphs. In the following, we provide the mathematical extension for
the bipartite case (the directed network formulation of our proposed model can
be considered a special case of the bipartite framework in which self-links
are removed and thus omitted from the below log-likelihood). For a bipartite
network with adjacency matrix $Y^{N_{1}\times N_{2}}$ we can formulate the
log-likelihood as:

$\displaystyle\log P(Y|\bm{\Lambda})$
$\displaystyle=\sum_{\begin{subarray}{c}i<j\\\
y_{i,j}=1\end{subarray}}\Bigg{(}\psi_{i}+\omega_{j}-||\mathbf{w}_{i}-\mathbf{v}_{j}||_{2}\Bigg{)}-\sum_{k_{L}=1}^{K_{L}}\Bigg{(}\sum_{i,j\in
C_{k_{L}}}\exp{\Big{(}\psi_{i}+\omega_{j}-||\mathbf{w}_{i}-\mathbf{v}_{j}||_{2}\Big{)}}\Bigg{)}$
$\displaystyle-\sum_{l=1}^{L}\sum_{k=1}^{K_{l}}\sum_{k^{{}^{\prime}}>k}^{K_{l}}\exp{\Big{(}-||\bm{\mu}_{k}^{(l)}-\bm{\mu}_{k^{\prime}}^{(l)}||_{2}\Big{)}}\Bigg{(}\sum_{i\in
C_{k}^{(l)}}\exp{(\psi_{i})}\Bigg{)}\Bigg{(}\sum_{j\in
C_{k^{{}^{\prime}}}^{(l)}}\exp{(\omega_{j})}\Bigg{)},$ (8)

where $\\{\bm{\mu}_{k}^{(l)}\\}_{k=1}^{K_{L}}$ are the latent centroids which
have absorbed the dependency of both sets of latent variables
$\\{\mathbf{w}_{i},\mathbf{v}_{j}\\}$ while we define the Poisson rate as:

$\lambda_{ij}=\exp\big{(}\psi_{i}+\omega_{j}-d(\mathbf{w}_{i},\mathbf{v}_{j})\big{)},$
(9)

where $\psi_{i}$ and $\omega_{j}$ are the corresponding random effects and
$\\{\mathbf{w}_{i}$, $\mathbf{v}_{j}\\}$ are the latent variables of the two
disjoint sets of the vertex set of sizes $N_{1}$ and $N_{2}$, respectively. In
this setting, we use our divisive Euclidean distance hierarchical clustering
procedure over the concatenation $\mathbf{z}=[\mathbf{w};\mathbf{v}]$ of the
two sets of latent variables. Therefore, we define an accurate hierarchical
block structure for bipartite networks, with each block including nodes from
both of the two disjoint modes. Here, a centroid is considered a leaf if the
corresponding tree-cluster contains less than $\log(N_{1})$ of the latent
variables $\\{\mathbf{w}_{i}\\}_{i=1}^{N_{1}}$ or less than $\log(N_{2})$ of
$\\{\mathbf{v}\\}_{j=1}^{N_{2}}$.

##  3 Experiments

We extensively evaluate the performance of our method compared to baseline
graph representation learning approaches on networks of various sizes and
structures. We have conducted all the experiments regarding the HBDM and HBDM-
Re on a $32$ GB Tesla V100 GPU machine. For the HBDM and HBDM-Re models, we
optimize the negative log-likelihood via the Adam (Kingma and Ba, 2017)
optimizer with learning rate $lr\in[0.01,0.1]$.

Datasets: We have performed the experiments on ten undirected networks of
various sizes and structures: a citation network (Cora (Sen et al., 2008)),
social interaction graphs (Facebook (Leskovec and Mcauley, 2012), YouTube
(Yang and Leskovec, 2015; Mislove et al., 2007), Flickr (Mislove et al.,
2007), Flixster (Zafarani and Liu, 2009)), product-label network (Amazon (Yang
and Leskovec, 2015)) and collaboration networks (Dblp (Perozzi et al., 2017),
AstroPh (Leskovec et al., 2007), GrQc (Leskovec et al., 2007), HepTh (Leskovec
et al., 2007)). Each network is considered as unweighted for the consistency
of the experiments. The detailed statistics of the networks are provided by
Table 1.

Table 1: Statistics of undirected networks. $N$: number of nodes, $|E|$:
number of edges.

| Cora | Dblp | AstroPh | GrQc | Facebook | HepTh | Amazon | YouTube | Flickr | Flixster
---|---|---|---|---|---|---|---|---|---|---
$N$ | 2,708 | 27,199 | 17,903 | 5,242 | 4,039 | 8,638 | 334,868 | 1,138,499 | 1,715,255 | 2,523,386
$|E|$ | 5,278 | 66,832 | 197,031 | 14,496 | 88,234 | 24,827 | 925,876 | 2,990,443 | 15,555,042 | 7,918,801

Baseline Methods: In our experiments, we have run various graph representation
learning methods in order to evaluate the performance of our approach. The
prominent GRL frameworks used in this study are: (i) DeepWalk (Perozzi et al.,
2014), (ii) Node2Vec (Grover and Leskovec, 2016), (iii) LINE (Tang et al.,
2015), (iv) NetMF (Qiu et al., 2018), (v) NetSMF (Qiu et al., 2019), (vi)
RandNE (Zhang et al., 2018), (vii) ProNE (Zhang et al., 2019), (viii)
LouvainNE (Bhowmick et al., 2020). For more details see the supplementary
material.

###  3.1 Link Prediction

We report results for the area under the curve of the receiver operator
characteristic (AUC). For the experimental setup, we follow the commonly
applied strategy (Perozzi et al., 2014; Grover and Leskovec, 2016) and remove
half of the edges. This strategy is not feasible for large-scale networks, so
we hide 30% of the edges for these networks. Links are removed in a way so
that the residual network stays connected. Extensive details for the link
prediction experiments, as well as, Precision-Recall AUC scores are given in
supplementary. Error bars for the following AUC scores were found to be on the
scale of $10^{-3}$ and thus negligible.

Moderate-Sized Networks: Results for the moderate-sized networks are given in
Table 2. The symbol ”-” indicates that the running time of the corresponding
model takes more than 20 hours and ”x” shows that the method is not able to
run due to insufficient memory space. We observe that the HBDM and HBDM-Re
perform significantly better or on par with the performance of the considered
baseline approaches. In particular, the HBDM and HBDM-Re perform better than
all the non-LDM baselines when $D=2$. It highlights the superiority of LDMs in
learning very low-dimensional network representations that accurately account
for the network structure. We further observe that representing degree
heterogeneity with random effects provides extended representational power as
the HBDM-Re consistently outperforms the HBDM. Comparing our framework with
the classic LDM-Re and LDM, we mostly see on par results experimentally which
we attribute to the hierarchical structure well preserving properties of
homophily and transitivity. In Figure 3a, we provide the analysis of the HBDM
likelihood when contrasted with its full likelihood estimation evaluated on
the moderate-sized network of Facebook (results for more networks are provided
in the supplementary material). We here observe that the HBDM likelihood
essentially approximates the true full likelihood providing systematically
slightly lower likelihood estimates which we attribute to the small structural
differences between calculating the distances analytically versus in a
hierarchical block manner. In addition, in Figure 3a we see fluctuations in
the likelihood which is an immediate result of building the network hierarchy
from scratch every $25$ iterations. Importantly, despite the fact that k-means
is notoriously known to be an NP-hard problem (Dasgupta, 2008; Mahajan et al.,
2012), we observe that rebuilding the hierarchy has a minimum effect on the
value of the likelihood, highlighting the stability of the inferred hierarchy
in the HBDM. In Figure 3b, we view the predictive performance as a function of
latent dimension $D$ and here in general attain modest improvements in the
predictive performance when increasing the embedding dimensions from $D=2$ to
$D=8$ with no further improvements increasing to $D=128$, highlighting the
efficiency in which HBDM and HBDM-Re utilize very low-dimensional
representations.

Table 2: AUC scores for representation sizes of $2$ and $8$ over moderate-
sized networks.

| AstroPh | GrQc | Facebook | HepTh | Cora | DBLP
---|---|---|---|---|---|---
Dimension ($D$) | $2$ | $8$ | $2$ | $8$ | $2$ | $8$ | $2$ | $8$ | $2$ | $8$ | $2$ | $8$
DeepWalk | .831 | .945 | .845 | .919 | .958 |  .986 | .773 | .874 | .684 | .782 | .803 | .939
Node2Vec | .825 | .950 | .809 | .884 | .914 | .988 | .780 | .881 | .640 | .776 | .803 | .945
LINE | .632 | .910 | .688 | .920 | .751 | .980 | .659 | .874 | .634 | .521 | .625 | .503
NetMF | .800 | .814 | .830 | .860 | .872 | .935 | .757 | .792 | .629 | .739 | .838 | .858
NetSMF | .828 | .891 | .756 | .805 | .907 | .976 | .705 | .810 | .605 | .737 | .766 | .857
RandNE | .524 | .554 | .534 | .560 | .614 | .657 | .519 | .509 | .508 | .556 | .508 | .517
LouvainNE | .798 | .813 | .861 | .868 | .957 | .958 | .774 | .874 | .767 | .747 | .900 | .904
ProNE | .768 | .907 | .818 | .883 | .900 | .971 | .678 | .823 | .675 | .764 | .813 | .924
LDM |  .925 | x | .915 | .943 |  .989 |  .992 | .855 | .919 |  .780 | .786 | .918 | x
LDM-Re | .943 | x |  .925 |  .944 | .990 | .992 |  .869 |  .917 | .770 | .787 | .926 | x
HBDM | .920 |  .960 | .917 |  .944 | .980 | .986 | .853 | .915 | .786 |  .792 |  .919 |  .956
HBDM-Re | .939 | .964 | .926 | .953 | .986 | .988 | .871 | .924 | .774 | .795 | .930 | .963

Table 3: AUC for varying representation sizes over the large-scale networks.

| YouTube | Flickr | Flixster | Amazon
---|---|---|---|---
Dimension ($D$) | $2$ | $3$ | $8$ | $2$ | $3$ | $8$ | $2$ | $3$ | $8$ | $2$ | $3$ | $8$
DeepWalk | .822 | .891 | .921 | .889 | .937 | .972 | .820 | .866 | .921 | .839 | .932 | .972
Node2Vec | - | - | - | - | - | - | - | - | - | .813 | .880 | .968
LINE | .660 | .832 | .878 | .685 | .889 | .812 | .523 | .868 | .936 | .626 | .501 | .500
NetMF | x | x | x | x | x | x | x | x | x | .829 | .831 | .858
NetSMF |  .939 |  .940 |  .949 |  .974 | .977 | .980 | .987 | .987 | .987 | .768 | .786 | .835
RandNE | .672 | .700 | .762 | .833 | .869 | .903 | .700 | .739 | .835 | .507 | .511 | .514
LouvainNE | .820 | .819 | .815 | .898 | .899 | .909 | .735 | .718 | .746 | .955 | .954 | .954
ProNE | .691 | .761 | .861 | .623 | .819 | .908 | .756 | .803 | .846 | .847 | .901 | .944
HBDM | .899 | .920 | .935 | .972 |  .979 |  .986 | .897 | .916 | .932 |  .974 |  .980 |  .988
HBDM-Re | .940 | .947 | .957 | .980 | .985 | .988 |  .962 |  .969 |  .971 | .976 | .981 | .988

(a) (a)

(b) (b)

(c) (c)

Figure 3: (a) NLL comparison between HBDM and LDM for Facebook with $D=2$. (b)
AUC-ROC performance over various networks for varying embedding sizes. (c) AUC
scores of HBDM-Re in terms of iterations for large-scale networks.

Large-Scale Networks: Results for the large-scale networks are given in Table
3. Again, we observe that HBDM-Re outperforms the baselines and it is only
surpassed in performance by NetSMF for the Flixster dataset. We also here find
that the inclusion of random effects in the LDMs improves performance
highlighting the importance of explicitly accounting for degree heterogeneity
also for large networks. Notably, we again detect very good performance for
the HBDM-Re, but also NetSMF when utilizing the very low embedding dimension
of $D=2$. In Figure 3c, we investigate the convergence of the best performing
HBDM-Re for the large networks, and we witness that the model rapidly
converges. After a couple of thousand iterations (very-scalable regime), we
already obtain competitive performance for link prediction, which then gently
increases until convergence.

Bipartite Networks: We validate the performance of our proposed framework for
bipartite structures by reporting the AUC score. We perform our experiments on
three bipartite networks: (1) Drug-Gene (Wishart et al., 2017) ($N_{1}=5,017$,
$N_{2}=2,324$, $|E|=15,138$), (2) GitHub (Chacon, 2009) ($N_{1}=56,519$,
$N_{2}=120,867$, $|E|=440,237$), and (3) Gottron-Reuters (Lewis et al., 2004)
($N_{1}=21,557$, $N_{2}=38,677$, $|E|=1,464,182$) following the same
experimental setting as in the undirected case of the moderate-sized networks
(network details are given in the supplementary). We provide the results in
Table 5 where we witness how the random-effects formulation of HBDM-Re
outperforms all the baselines and in most cases significantly. Another
interesting observation is that the random-effects model has a notably higher
performance than the corresponding global bias model, as the three studied
networks have high degree heterogeneity.

Table 4: $\text{Micro-F}_{1}$ scores varying dimension sizes for two moderate
and large-scale networks.

| Cora | DBLP | Amazon | YouTube
---|---|---|---|---
Dimension ($D$) | $2$ | $3$ | $128$ | $2$ | $3$ | $128$ | $2$ | $3$ | $8$ | $2$ | $3$ | $8$
DeepWalk | .502 | .712 |  .838 | .519 | .605 | .822 | .231 | .596 | .929 | .293 | .351 | .413
Node2Vec | .419 | .658 | .835 | .448 | .540 | .815 | .096 | .305 | .895 | - | - | -
LINE | .197 | .191 | .794 | .328 | .294 | .771 | .005 | .003 | .003 | .185 | .134 | .177
NetMF | .389 | .653 | .835 | .654 | .707 | .742 | x | x | x | x | x | x
NetSMF | .554 | .705 | .842 | .622 | .732 | .829 | .387 | .649 | .845 | .317 | .361 | .397
RandNE | .271 | .337 | .731 | .406 | .473 | .718 | .223 | .411 | .787 | .211 | .226 | .277
LouvainNE | .804 |  .811 | .801 | .780 | .812 |  .825 | .970 | .971 | .974 | .362 | .360 | .359
ProNE | .450 | .611 | .830 | .574 | .634 | .825 | .420 | .750 |  .933 | .218 | .274 | .379
LDM | .810 | .802 | .774 | x | x | x | x | x | x | x | x | x
LDM-Re | .802 | .803 | .796 | x | x | x | x | x | x | x | x | x
HBDM | .789 | .807 | .816 |  .812 |  .814 | .772 | .970 | .971 | .931 | .320 |  .366 |  .414
HBDM-Re |  .805 | .813 | .818 | .805 | .822 | .808 |  .956 |  .955 | .931 |  .326 | .367 | .414

###  3.2 Node classification

We assess the performance of the proposed framework in the uni/multi-label
classification task and provide the $\text{Micro-F}_{1}$ scores in Table 4
($\text{Macro-F}_{1}$ scores are reported in the supplementary). For the
experimental setup, we randomly pick $50\%$ of nodes as the training set and
use the rest as the testing set. For an accurate comparison across different
methods, we used two simple classifiers, a linear (logistic/multinomial
regression classifier) and a non-linear (linearithmic k-nearest neighbors
(kNN) classifier) and report the highest scores. We found that all methods
benefit from using kNN. The number of neighbors was set to $k=10$ (similar
results were obtained with higher choices for $k$ as well). Lastly, we report
the average $\text{Micro-F}_{1}$ scores across $10$ repeated trials. Results
for the uni-labeled Cora and DBLP networks are reported in the two leftmost
columns of Table 4. We observe that HBDM-Re and HBDM significantly outperform
the baselines for the regimes of $D=2,3$ with only LouvainNE being
competitive. Results for large-scale and multi-labeled networks Amazon and
YouTube are provided by the two rightmost columns in Table 4. Again, the
proposed framework outperforms the baselines for the low-dimensional regime
with LouvainNE being on par. Comparing our framework with the classic LDM-Re
and LDM, we again see more or less on par performance which we attribute to
the HBDM well preserving the intrinsic properties of homophily and
transitivity.

Table 5: AUC scores for varying representation sizes over three bipartite
networks.

| Drug-Gene | GitHub | Gottron-Reuters
---|---|---|---
Dimension ($D$) | $2$ | $3$ | $8$ | $2$ | $3$ | $8$ | $2$ | $3$ | $8$
DeepWalk | .673 |  .843 |  .878 | .762 | .853 | .902 | .673 | .769 | .905
Node2Vec | .758 | .814 | .793 | .724 | .823 | .876 | .694 | .766 | .830
LINE | .798 | .836 | .867 | .805 | .766 | .902 | .715 | .696 | .850
NetMF | .576 | .598 | .742 | .711 | .711 | .708 | .747 | .747 | .730
NetSMF |  .839 | .838 | .796 |  .846 | .847 | .857 | .874 |  .934 |  .941
RandNE | .536 | .551 | .613 | .615 | .651 | .707 | .769 | .808 | .872
LouvainNE | .760 | .767 | .779 | .694 | .702 | .735 | .654 | .648 | .679
ProNE | .667 | .765 | .831 | .676 | .771 | .840 | .606 | .725 | .909
HBDM | .726 | .779 | .872 | .791 |  .855 |  .905 |  .894 | .916 | .936
HBDM-Re | .846 | .888 | .917 | .917 | .918 | .924 | .964 | .967 | .969

(a)  (b) Cora

(c) DBLP (d) $D$

(e) $2$ (TES)

(f) $128$ (t-SNES)

(g) $2$ (TES)

(h) $128$ (t-SNES)

Node2Vec

(i)

(j) ($\text{NR}=0.863$)

(k)

(l) ($\text{NR}=0.920$)

(m)

NetSMF

(n)

(o) ($\text{NR}=0.860$)

(p)

(q) ($\text{NR}=0.940$)

(r)

HBDM-Re

(s)

(t) ($\text{NR}=\mathbf{0.975}$)

(u)

(v) ($\text{NR}=\mathbf{0.985}$)

(w)

Figure 4: 2-D True Embedding Space versus 128-D t-SNE constructed space. For
TES we provide the AUC-ROC for the network reconstruction (NR) task.

(a) (i) Dendrogram (b) (ii) Embedding Space

(c) (iii) L=1

(d) (iv) L=3

(e) (v) L=5

(f) (vi) L=8

Figure 5: Amazon network dendrogram, embedding space and ordered adjacency
matrices for the learned $D=2$ embeddings of HBDM-Re and various levels $(L)$
of the hierarchy.

###  3.3 Network visualization

The graph representation learning literature mainly focuses on embeddings with
dimensionality greater than $D=2\text{ and }3$. As a direct consequence,
network visualizations rely on dimensionality reduction frameworks, typically
using the t-distributed Stochastic Neighbor Embedding (t-SNE) (van der Maaten
and Hinton, 2008). In order to verify the quality of the t-SNE constructed
Space (t-SNES), we provide the labeled-colored True Embedding Space (TES) in
Figure 4 for $D=2$, as well as for $D=128$ mapped to $D=2$ via the use of
t-SNE for Cora and DBLP. We see that the HBDM-Re frameworks provide highly
informative embeddings with no need for dimensionality reduction unlike the
rest of the baselines. This is also verified from the optimal performance in
network reconstruction, HBDM-Re can successfully express the network structure
using just $D=2$. In Figure 5 we provide the hierarchical block structure
constructed by the HBDM-Re for the Amazon network. For visualization, we used
the average within-cluster Euclidean distance to the centroid
$\Big{(}\Delta\\{A,B\\}=\frac{1}{N_{A}+N_{B}}\sum_{i\in
C_{A},C_{B}}|\bm{z}_{i}-\bm{\mu}_{A\bigcup B}|_{2}\Big{)}$, as a linkage
function to form a post-processing agglomerative clustering, for ordering the
initial $\log N$ centroids. In Figure 5 (i), we provide the dendrogram which
denotes the agglomeration result in the top-level with red lines. The
dendrogram continues with the hierarchical splits of our HBDM-Re where each
color indicates the initial $\log N$ blocks. The y-axis of the dendrogram
represents the binary logarithm of the Sum of Euclidean Distances,
$\textsl{Log2-SED}=\log_{2}\Big{(}\sum_{i\in
C_{k}^{(l)}}|\bm{z}_{i}-\bm{\mu}_{k}^{(l)}|_{2}\Big{)}$. Moreover, Figure 5
(ii) conveys the corresponding latent space, colored based on the coarse $\log
N$ split, revealing directly interpretable and accurate network
representations. In Figure 5 (iii), (iv), (v) and (vi) we showcase the
organized adjacency matrices, based on the 2-dimensional HBDM-Re learned
hierarchy for various levels $L$ of the tree. We here, observe the
representation power of the extracted hierarchy from just a $2$-dimensional
HBDM-Re defining communities and their sub-communities at finer and finer
details.

For the bipartite case, we show how HBDM-Re can enhance our understanding of
the bipartite structure at multiple scales and levels. Similar to the
undirected case, Figure 6 (i), indicates the dendrogram of the imposed
hierarchy, enriched with agglomeration for a coarse level block ordering and
proximity for the GitHub network. In addition, Figure 6 (ii), provides the
corresponding latent space, colored based on the coarse $\log N$ split.
Notably, no dimensionality-reduction is necessary to define accurate network
representation in the latent space of the two disjoint populations and
visually access and express node similarity. In Figure 6 (iii), (iv) and (v),
we exhibit how the multi-scale structure evolves through different levels of
the hierarchy defined by HBDM-Re, showcasing how a joint bi-clustering for
complex network embeddings naturally can be obtained, with no need for post-
processing steps. Our HBDM, can thus accurately characterize bipartite
networks and successfully uncover their hierarchical block structure
efficiently.

(a) (i) Dendrogram (b) (ii) Embedding Space

(c) (iii) L=1

(d) (iv) L=3

(e) (v) L=6

Figure 6: GitHub network dendrogram, embedding space and ordered adjacency
matrices for the learned $D=2$ embeddings of HBDM-Re and various levels $(L)$
of the hierarchy.

##  4 Discussion

We developed the HBDM, a scalable reconciliation of latent distance models and
their ability to account for homophily and transitivity with hierarchical
representations of network structures. We demonstrated how the proposed HBDM
provides favorable network representations by: (1) Operating with a Euclidean
distance metric providing an intuitive human perception of node similarity.
(2) Naturally representing multiscale hierarchical structure based on its
block structure and carefully designed clustering procedure optimized in terms
of Euclidean distances. (3) Directly and consistently operating in $D=2,3$
with high performance. (4) Performing well on all considered downstream tasks
highlighting its ability to account for the underlying network structure.
Importantly, the inferred hierarchical structure admits community discovery at
multiple scales as highlighted from the inferred dendrograms and ordered
adjacency matrices, and naturally extends to the characterization of
communities of bipartite networks.

Our finding of ultra-low dimensional accurate characterizations of network
structures supports the findings in Chanpuriya et al. (2020) in which a
logistic PCA model was found to enable exact low-dimensional recovery of
multiple real-world networks. Whereas the work of Chanpuriya et al. (2020)
focuses on exact network reconstruction we find that generalizable patterns
can be well extracted in ultra-low dimensional representations with
performance saturating after just $D=8$ dimensions for all networks
considered. Whereas Chanpuriya et al. (2020) found that their low-dimensional
space did not perform well in classification tasks we observed strong node
classification performance by the low-dimensional representations provided by
HBDM. Importantly, for node classification, we observed better performance
using KNN as opposed to simple linear classification based on
logistic/multinomial regression typically used for node classification. This
highlights that whereas most GRL works use linear classifiers there is no
guarantee that the embedding space will be linearly separable and performance
should therefore be compared to non-linear classifiers as they may provide
more favorable performance as observed in this study.

Recent pioneering works (Ben Chamberlain and Clough, 2017; Nickel and Kiela,
2017) have drawn significant attention of the research community by
questioning the conventional embedding space preference. It is well known that
many real-world networks show power-law degree distribution, or they can
consist of latent hierarchical inner structures. Therefore, the Euclidean
space might not always be appropriate to represent such complex network
architectures. It might also require higher-dimensional spaces to show
comparable performance in the GRL tasks. The works of Ben Chamberlain and
Clough (2017); Nickel and Kiela (2017) demonstrated that hyperbolic spaces,
such as the Poincare disk model, can provide substantial benefits over the
Euclidean space. The presented HBDM model naturally extends to other distance
measures and future studies should explore how the HBDM can be extended to
hierarchical representations beyond Euclidean geometry.

Covariate information plays an important role in the outstanding performance
of GRL methods and especially GNNs. In the current LSM literature, side
information is accounted for by extra regressors in the logit/log link
functions expressing the likelihood of a dyad being connected. Using the
Mahalanobis distance imposing a block-diagonal covariance matrix (see
supplementary), the proposed HBDM can naturally incorporate covariate
information directly to the latent space and notably construct multi-scale
structures via the enriched and concatenated embedding of the latent variables
and the covariate information.

Our discoveries highlight the existence and importance of hierarchical multi-
scale structures in complex networks. The across hierarchy re-ordered
adjacency matrices given by HBDM, manifest sub-communities inside of what
already appears as a strongly connected community. This points to how delicate
the task of defining communities is and the importance of accounting for
communities at multiple scales, as enabled by the HBDM. Importantly, these
results generalize for bipartite networks where multi-scale geometric
representations, joint hierarchical structures, and community discovery are
arduous tasks.

The HBDM uses the LDM and thus is good at characterizing transitivity and
homophily at a node and cluster level, whereas the random effects enable
accounting for degree heterogeneity. Notably, the HBDM suffers from the
limitations of the LDM and is thus unable to model stochastic equivalence.
Future work should therefore investigate hierarchical structures imposed on
more flexible GRL procedures enabling stochastic equivalence and contrast the
performance when accounting for stochastic equivalence to the existing
hierarchical methods based on the SBM (Clauset et al., 2008; Roy et al., 2007;
Herlau et al., 2012; Blundell and Teh, 2013; Herlau et al., 2013; Peixoto,
2014).

In conclusion, we proposed the Hierarchical Block Distance Model (HBDM), a
scalable reconciliation of network embeddings using the latent distance model
(LDM) and hierarchical characterizations of structure at multiple scales via a
novel clustering framework. Notably, the model mimics the behavior of the LDM
where the use of homophily and transitivity is most important while scaling in
complexity by $\mathcal{O}(DN\log{N})$. We analyzed thirteen networks from
moderate sizes to large-scale with the HBDM having favorable performance when
compared to existing scalable embedding procedures. In particular, we observed
that the HBDM well predicts links and node classes providing accurate network
visualizations and characterization of structure at multiple scales. Our
results demonstrate that favorable performance can be achieved using ultra-low
(i.e. $D=2$) embedding dimensions and a scalable hierarchical representation
that accounts for homophily and transitivity.

## References

  * (1)
  * Ahn et al. (2010) Ahn, Y.-Y., Bagrow, J. P. and Lehmann, S. (2010), ‘Link communities reveal multiscale complexity in networks’, nature 466(7307), 761–764.
  * Barabási and Pósfai (2016) Barabási, A.-L. and Pósfai, M. (2016), Network science, Cambridge University Press.
  * Ben Chamberlain and Clough (2017) Ben Chamberlain, M. D. and Clough, J. (2017), Neural embeddings of graphs in hyperbolic space, in ‘MLG Workshop’.
  * Bhowmick et al. (2020) Bhowmick, A. K., Meneni, K., Danisch, M., Guillaume, J.-L. and Mitra, B. (2020), LouvainNE: Hierarchical louvain method for high quality and scalable network embedding, in ‘WSDM’, pp. 43–51.
  * Blondel et al. (2008) Blondel, V. D., Guillaume, J.-L., Lambiotte, R. and Lefebvre, E. (2008), ‘Fast unfolding of communities in large networks’, J. Stat. Mech. Theory Exp. 2008(10), P10008.
  * Blundell and Teh (2013) Blundell, C. and Teh, Y. W. (2013), Bayesian hierarchical community discovery, in ‘NIPS’, Vol. 26.
  * Cai and Liu (2020) Cai, Q. and Liu, J. (2020), ‘Hierarchical clustering of bipartite networks based on multiobjective optimization’, IEEE Trans. Netw. Sci. Eng. 7(1), 421–434.
  * Cao et al. (2015) Cao, S., Lu, W. and Xu, Q. (2015), GraRep: Learning graph representations with global structural information, in ‘CIKM’, pp. 891–900.
  * Çelikkanat and Malliaros (2020) Çelikkanat, A. and Malliaros, F. D. (2020), Exponential family graph embeddings, in ‘AAAI’, pp. 3357–3364.
  * Chacon (2009) Chacon, S. (2009), ‘2009 github challenge’.
https://github.blog/2009-07-29-the-2009-github-contest/

  * Chanpuriya et al. (2020) Chanpuriya, S., Musco, C., Sotiropoulos, K. and Tsourakakis, C. E. (2020), ‘Node embeddings and exact low-rank representations of complex networks’, CoRR abs/2006.05592.
  * Chen et al. (2018) Chen, H., Perozzi, B., Hu, Y. and Skiena, S. (2018), HARP: hierarchical representation learning for networks, in ‘AAAI’, pp. 2127–2134.
  * Clauset et al. (2008) Clauset, A., Moore, C. and Newman, M. E. (2008), ‘Hierarchical structure and the prediction of missing links in networks’, Nature 453(7191), 98–101.
  * Dasgupta (2008) Dasgupta, S. (2008), The hardness of k-means clustering.
  * Dhillon (2001) Dhillon, I. S. (2001), Co-clustering documents and words using bipartite spectral graph partitioning, in ‘KDD’, p. 269–274.
  * Epp (2010) Epp, S. S. (2010), Discrete Mathematics with Applications, 4th edn, Brooks/Cole, USA.
  * Fortunato (2010) Fortunato, S. (2010), ‘Community detection in graphs’, Physics Reports 486(3), 75–174.
  * Friel et al. (2016) Friel, N., Rastelli, R., Wyse, J. and Raftery, A. E. (2016), ‘Interlocking directorates in irish companies using a latent space model for bipartite networks’, PNAS 113(24), 6629–6634.
  * Getoor and Taskar (2007) Getoor, L. and Taskar, B. (2007), Introduction to Statistical Relational Learning (Adaptive Computation and Machine Learning), The MIT Press.
  * Grover and Leskovec (2016) Grover, A. and Leskovec, J. (2016), Node2Vec: Scalable feature learning for networks, in ‘KDD’, pp. 855–864.
  * Hamilton et al. (2017) Hamilton, W. L., Ying, R. and Leskovec, J. (2017), Inductive representation learning on large graphs, in ‘NIPS’.
  * Handcock et al. (2007a) Handcock, M., Raftery, A. and Tantrum, J. (2007a), ‘Model-based clustering for social networks’, J R Stat Soc Ser A Stat Soc 170, 301 – 354.
  * Handcock et al. (2007b) Handcock, M., Raftery, A. and Tantrum, J. (2007b), ‘Model-based clustering for social networks’, J. R. Stat. Soc. 170, 301 – 354.
  * Hartigan and Wong (1979) Hartigan, J. A. and Wong, M. A. (1979), ‘Algorithm AS 136: A K-Means clustering algorithm’, Applied Statistics 28(1), 100–108.
  * Herlau et al. (2013) Herlau, T., Mørup, M. and Schmidt, M. (2013), Modeling temporal evolution and multiscale structure in networks, in ‘ICML’, pp. 960–968.
  * Herlau et al. (2012) Herlau, T., Mørup, M., Schmidt, M. N. and Hansen, L. K. (2012), Detecting hierarchical structure in networks, in ‘CIP’, IEEE, pp. 1–6.
  * Herlau et al. (2014) Herlau, T., Schmidt, M. N. and Mørup, M. (2014), ‘Infinite-degree-corrected stochastic block model’, Physical review E 90(3), 032819.
  * Hoff (2005) Hoff, P. D. (2005), ‘Bilinear mixed-effects models for dyadic data’, JASA 100(469), 286–295.
  * Hoff (2007) Hoff, P. D. (2007), Modeling homophily and stochastic equivalence in symmetric relational data, in ‘NIPS’, p. 657–664.
  * Hoff et al. (2002) Hoff, P. D., Raftery, A. E. and Handcock, M. S. (2002), ‘Latent space approaches to social network analysis’, JASA 97(460), 1090–1098.
  * Holland et al. (1983) Holland, P. W., Laskey, K. B. and Leinhardt, S. (1983), ‘Stochastic blockmodels: First steps’, Social networks 5(2), 109–137.
  * Karrer and Newman (2011) Karrer, B. and Newman, M. E. (2011), ‘Stochastic blockmodels and community structure in networks’, Physical review E 83(1), 016107.
  * Kim et al. (2017) Kim, B., Lee, K., Xue, L. and Niu, X. (2017), ‘A review of dynamic network models with latent variables’.
  * Kingma and Ba (2017) Kingma, D. P. and Ba, J. (2017), ‘Adam: A method for stochastic optimization’.
  * Krivitsky et al. (2009) Krivitsky, P. N., Handcock, M. S., Raftery, A. E. and Hoff, P. D. (2009), ‘Representing degree distributions, clustering, and homophily in social networks with latent cluster random effects models’, Social Networks 31(3), 204 – 213.
  * Leskovec et al. (2007) Leskovec, J., Kleinberg, J. and Faloutsos, C. (2007), ‘Graph evolution: Densification and shrinking diameters’, ACM Trans. Knowl. Discov. Data 1(1).
  * Leskovec and Mcauley (2012) Leskovec, J. and Mcauley, J. J. (2012), Learning to discover social circles in ego networks, in ‘NIPS’, pp. 539–547.
  * Lewis et al. (2004) Lewis, D. D., Yang, Y., Rose, T. G. and Li, F. (2004), ‘Rcv1: A new benchmark collection for text categorization research’, J. Mach. Learn. Res. 5, 361–397.
  * Li et al. (2020) Li, T., Lei, L., Bhattacharyya, S., Van den Berge, K., Sarkar, P., Bickel, P. J. and Levina, E. (2020), ‘Hierarchical community detection by recursive partitioning’, JASA pp. 1–18.
  * Liben-Nowell and Kleinberg (2003) Liben-Nowell, D. and Kleinberg, J. (2003), The link prediction problem for social networks, in ‘CIKM’, p. 556–559.
  * Mahajan et al. (2012) Mahajan, M., Nimbhorkar, P. and Varadarajan, K. (2012), ‘The planar k-means problem is np-hard’, Theoretical Computer Science 442, 13–21.
  * Mikolov et al. (2013) Mikolov, T., Sutskever, I., Chen, K., Corrado, G. and Dean, J. (2013), Distributed representations of words and phrases and their compositionality, in ‘NIPS’, pp. 3111–3119.
  * Mislove et al. (2007) Mislove, A., Marcon, M., Gummadi, K. P., Druschel, P. and Bhattacharjee, B. (2007), Measurement and analysis of online social networks, in ‘IMC’.
  * Mørup and Schmidt (2012) Mørup, M. and Schmidt, M. N. (2012), ‘Bayesian community detection’, Neural computation 24(9), 2434–2456.
  * Newman (2003) Newman, M. E. J. (2003), ‘The structure and function of complex networks’, SIAM Review 45(2), 167–256.
  * Nguyen and Malliaros (2018) Nguyen, D. and Malliaros, F. D. (2018), BiasedWalk: Biased sampling for representation learning on graphs, in ‘Big Data’, pp. 4045–4053.
  * Nickel and Kiela (2017) Nickel, M. and Kiela, D. (2017), Poincaré embeddings for learning hierarchical representations, in ‘NIPS’, Vol. 30.
  * Nowicki and Snijders (2001) Nowicki, K. and Snijders, T. A. B. (2001), ‘Estimation and prediction for stochastic blockstructures’, JASA 96(455), 1077–1087.
  * Peixoto (2014) Peixoto, T. P. (2014), ‘Hierarchical block structures and high-resolution model selection in large networks’, Physical Review X 4(1).
  * Perozzi et al. (2014) Perozzi, B., Al-Rfou, R. and Skiena, S. (2014), ‘DeepWalk: Online learning of social representations’, CoRR abs/1403.6652.
  * Perozzi et al. (2017) Perozzi, B., Kulkarni, V., Chen, H. and Skiena, S. (2017), Don’t walk, skip! online learning of multi-scale network embeddings, in ‘ASONAM’, pp. 258–265.
  * Qiu et al. (2019) Qiu, J., Dong, Y., Ma, H., Li, J., Wang, C., Wang, K. and Tang, J. (2019), NetSMF: Large-scale network embedding as sparse matrix factorization, in ‘WWW’, pp. 1509–1520.
  * Qiu et al. (2018) Qiu, J., Dong, Y., Ma, H., Li, J., Wang, K. and Tang, J. (2018), Network embedding as matrix factorization: Unifying DeepWalk, LINE, PTE, and Node2Vec, in ‘WSDM’.
  * Raftery et al. (2012) Raftery, A., Niu, X., Hoff, P. and Yeung, K. Y. (2012), ‘Fast inference for the latent space network model using a case-control approximate likelihood’, J Comput Graph Stat 21.
  * Ravasz and Barabási (2003) Ravasz, E. and Barabási, A.-L. (2003), ‘Hierarchical organization in complex networks’, Phys. Rev. E 67, 026112.
  * Rosvall and Bergstrom (2007) Rosvall, M. and Bergstrom, C. T. (2007), ‘An information-theoretic framework for resolving community structure in complex networks’, PNAS 104(18), 7327–7331.
  * Roy et al. (2007) Roy, D. M., Kemp, C., Mansinghka, V. and Tenenbaum, J. (2007), Learning annotated hierarchies from relational data, in ‘NIPS’, Vol. 19.
  * Roy and Teh (2008) Roy, D. M. and Teh, Y. (2008), The mondrian process, in D. Koller, D. Schuurmans, Y. Bengio and L. Bottou, eds, ‘NIPS’, Vol. 21, Curran Associates, Inc.
  * Sen et al. (2008) Sen, P., Namata, G., Bilgic, M., Getoor, L., Gallagher, B. and Eliassi-Rad, T. (2008), ‘Collective classification in network data’, AI magazine .
  * Sosa and Buitrago (2020) Sosa, J. and Buitrago, L. (2020), ‘A review of latent space models for social networks’, CoRR abs/2012.02307.
  * Tang et al. (2015) Tang, J., Qu, M., Wang, M., Zhang, M., Yan, J. and Mei, Q. (2015), LINE: Large-scale information network embedding, in ‘WWW’, pp. 1067–1077.
  * Tsutsu and Morikawa (2012) Tsutsu, H. and Morikawa, Y. (2012), An lp norm minimization using auxiliary function for compressed sensing, in ‘Proc. Int. Multiconf. Comp. Sci. Inf. Technol.’.
  * van der Maaten and Hinton (2008) van der Maaten, L. and Hinton, G. (2008), ‘Visualizing data using t-sne’, Journal of Machine Learning Research 9(86), 2579–2605.
  * White et al. (1976) White, H. C., Boorman, S. A. and Breiger, R. L. (1976), ‘Social structure from multiple networks. i. blockmodels of roles and positions’, American journal of sociology 81(4).
  * Wieling and Nerbonne (2010) Wieling, M. and Nerbonne, J. (2010), Hierarchical spectral partitioning of bipartite graphs to cluster dialects and identify distinguishing features, in ‘TextGraphs-5’, pp. 33 – 41.
  * Wind and Mørup (2012) Wind, D. K. and Mørup, M. (2012), Link prediction in weighted networks, in ‘2012 IEEE International Workshop on Machine Learning for Signal Processing’, pp. 1–6.
  * Wishart et al. (2017) Wishart, D. S., Feunang, Y. D., Guo, A. C., Lo, E. J., Marcu, A., Grant, J. R., Sajed, T., Johnson, D., Li, C., Sayeeda, Z., Assempour, N., Iynkkaran, I., Liu, Y., Maciejewski, A., Gale, N., Wilson, A., Chin, L., Cummings, R., Le, D., Pon, A., Knox, C. and Wilson, M. (2017), ‘DrugBank 5.0: a major update to the DrugBank database for 2018’, Nucleic Acids Research 46(D1), D1074–D1082.
  * Yang and Leskovec (2015) Yang, J. and Leskovec, J. (2015), ‘Defining and evaluating network communities based on ground-truth’, Knowledge and Information Systems 42(1), 181–213.
  * Zafarani and Liu (2009) Zafarani, R. and Liu, H. (2009), ‘Social computing data repository at ASU’.
  * Zhang et al. (2020) Zhang, D., Yin, J., Zhu, X. and Zhang, C. (2020), ‘Network representation learning: A survey’, IEEE Transactions on Big Data 6(1), 3–28.
  * Zhang et al. (2019) Zhang, J., Dong, Y., Wang, Y., Tang, J. and Ding, M. (2019), Prone: Fast and scalable network representation learning, in ‘IJCAI’, pp. 4278–4284.
  * Zhang et al. (2018) Zhang, Z., Cui, P., Li, H., Wang, X. and Zhu, W. (2018), Billion-scale network embedding with iterative random projection, in ‘ICDM’, pp. 787–796.