{'fraction_non_alphanumeric': 0.030826970363659893, 'fraction_numerical': 0.011310987446615763, 'mean_word_length': 4.997516299285936, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Hierarchical Block Distance Model for Ultra Low-Dimensional Graph Representations

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.


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 & Kleinberg 2003), node classification (Getoor & Taskar 2007, Grover & 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 & Leskovec 2016, Tang et al. 2015, Ç elikkanat & Malliaros 2020, Nguyen & 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 & 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 viceversa. 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 & Snijders 2001, i.e. ("groups of nodes defined by shared intra-and inter-group relationships") defining non positive-semidefinite 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 & 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 & Teh 2008, Herlau et al. 2012, Peixoto 2014) as well as agglomerative (Blondel et al. 2008, Ahn et al. 2010, Blundell & 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 & 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 & 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 & Bergstrom (2007),  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 hierar-chical 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) 1 , 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., 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.

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 G = (V, E) be a graph where N := |V | is the number of nodes and Y N ×N = (y i,j ) ∈ {0, 1} N ×N be the adjacency matrix of the graph such that y i,j = 1 if the pair (i, j) ∈ E otherwise it is equal 0, for all 1 ≤ i < j ≤ N . We denote the latent representations of nodes by Z = (z i,d ) ∈ R N ×D where each row vector, z i ∈ R D , indicates the corresponding embedding of node i ∈ 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 R D -dimensional latent space in which every node of the graph is characterized through the unobserved but informative node-specific variables {z i ∈ 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: where θ 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 & 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 & 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: where γ i ∈ R denotes the node-specific random-effects and d ij (·, ·) denotes any distance metric obeying the triangle inequality variables {z i } i∈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, γ g , corresponds to the special case in which γ i = γ j = 0.5γ g . As in Hoff et al. (2002), we presently adopt the Euclidean distance as the choice for the distance metric d ij (·, ·).

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} ∈ V × V based on the Euclidean distance, as shown in Equation (2). Therefore, we can define a blockalike 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 loglikelihood as of the standard LDM, defining a sum over each ordered pair of the network, as: For brevity, we presently ignore the linear scaling by dimensionality D of the above loglikelihood function. Notably, the link contribution y i,j =1 log(λ 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 & 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 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 & Snijders (2001), when grouping the nodes into K clusters we define the rate between block k and k in terms of their distance between centroids. A simple block structure without a hierarchy would lead to the following non-link expression: where µ k denotes the k'th cluster centroid of the set C = {C 1 , . . . , C K }, and has absorbed the dependency over the variables Z ∈ R N ×D . More specifically, the cluster centroids µ 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 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 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 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  We can thereby define a Hierarchical Block Distance Model with Random Effects (HBDM-Re) as: where l ∈ {1, . . . , L} denotes the l'th dendrogram level, k l is the index representing the cluster id for the different tree levels and µ (l) k the corresponding centroid. We also consider a Hierarchical Block Distance Model (HBDM) without the random effects setting γ i = 0.5γ g . For a multifurcating tree splitting in K clusters and having N/log(N ) terminal nodes (clusters), the number of internal nodes are O(N/(K log N )) and each node needs to evaluate O(K 2 ) pairs providing an overall complexity of O(N K/ log N ), thus K ≤ log N 2 to achieve O(N log N ) scaling (Epp 2010).

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 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 ||µ k l − µ k l || 2 . The divisive clustering procedure thus relies on the following Euclidean norm objective where k denotes the cluster id, z i is the i'th data observation, r ik the cluster responsibility/assignment, and µ 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 & Morikawa (2012) developed in the context of compressed sensing. We define an auxiliary function for (6) as: where φ are the auxiliary variables. Thereby, minimizing Equation (7) with respect to φ nk yields φ * ik = ||z i − µ k || 2 and by plugging φ * ik back to (6) we obtain J + (φ * , r, µ) = J(r, µ) verifying that (7) is indeed a valid auxiliary function for (6). The algorithm proceeds by optimizing cluster centroids as and assigning points to centroids as arg min The overall complexity of this procedure is O(T KN D) (Hartigan & 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 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 O(log N/ log log N ) and  Figure 1 (i), illustrates the resulting tree 2 for a small problem of N = 64 nodes in which we first split into 4 (≈ log(64)) clusters and subsequently create binary splits until each leaf cluster contains 4 (≈ log(64)) nodes.

Expressing Homophily and Transitivity
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 ap-2 For visualization purposes only, we show equally sized clusters.
(i) Non-optimal rotation over leaf clusters.
(ii) Optimal rotation over leaf clusters.
(iii) Three latent block structures. Gray lines represent the links and the dashed lines the distance between the cluster centers. proximated 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)  Lemma 2.1. Let G = (V, E) be a graph and let C be a cluster with its centroid located at µ ∈ R D having an edge (i, j) ∈ E for some i ∈ C and j ∈ V\C such that z i = µ. If is a rotation matrix acting on the embeddings of nodes in cluster C, then the maximum degree of freedom of any infinitesimal λ ij -invariant rotation is defined by θ ∈ 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.

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 ×N 2 we can formulate the log-likelihood as: where {µ (l) k } K L k=1 are the latent centroids which have absorbed the dependency of both sets of latent variables {w i , v j } while we define the Poisson rate as: where ψ i and ω j are the corresponding random effects and {w i , 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 z = [w; 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 {w i } N 1 i=1 or less than log(N 2 ) of {v} N 2 j=1 .

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 con-

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 & 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.
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. 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.

Node classification
We assess the performance of the proposed framework in the uni/multi-label classification task and provide the Micro-F 1 scores in Table 4 (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 nonlinear (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  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.

Network visualization
The graph representation learning literature mainly focuses on embeddings with dimensionality greater than D = 2 and 3. As a direct consequence, network visualizations rely on dimensionality reduction frameworks, typically using the t-distributed Stochastic Neighbor  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 as a linkage function to form a postprocessing 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, Log2-SED = log 2 i∈C (l) k |z i − µ (l) k | 2 . 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 2dimensional 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.

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 opti- ing 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 lowdimensional 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. 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 & 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 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.