{'fraction_non_alphanumeric': 0.029117709245778837, 'fraction_numerical': 0.014362521521128463, 'mean_word_length': 4.485461022284815, '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': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comparing directed networks via denoising graphlet distributions

Network comparison is a widely-used tool for analyzing complex systems, with applications in varied domains including comparison of protein interactions or highlighting changes in structure of trade networks. In recent years, a number of network comparison methodologies based on the distribution of graphlets (small connected network subgraphs) have been introduced. In particular, NetEmd has recently achieved state of the art performance in undirected networks. In this work, we propose an extension of NetEmd to directed networks and deal with the significant increase in complexity of graphlet structure in the directed case by denoising through linear projections. Simulation results show that our framework is able to improve on the performance of a simple translation of the undirected NetEmd algorithm to the directed case, especially when networks differ in size and density.


Introduction
Complex networks represent relationships between actors of systems with non-trivial structural properties and are ubiquitous in a myriad of domains. Studying the networks underlying these complex systems is then a step towards gaining an understanding about the systems themselves. Comparing different objects is a fundamental part of human cognition, therefore it is only natural that analysis of networks encompasses a comparison element. It is sometimes straightforward to show two networks are different, through finding a mismatch between the set of nodes and edges that compose the networks, which would and claim a difference if these sets do not match. On the other hand, knowing if two networks are exactly the same is known as the graph isomorphism problem, which has been shown to be in the NP complexity class [6]. Between these two extremes, network comparison has been established as an area of study within network analysis that combines network statistics and computable properties of a network to tell how similar or dissimilar two networks are. Network comparison has been the subject of an increasing amount of research, see for example [1,2,41,43,46], and its applications are curate representations of that phenomena due to unseen data, measurement errors or simply because the systems they represent are too complex to fully describe. These inaccuracies lead to noise in the statistics used to describe the network and similar topographies become harder to recognize, especially when networks differ greatly in size.
The distributions of graphlets are not immune to these errors, so our second contribution in this work is using denoising methods to reduce the impact of these errors so that NetEmd is able to more accurately distinguish networks according to those core structures that compose each network. To this end, we propose a framework that employs Principal Component Analysis (PCA) [13] as the most commonly used denoising method before comparing the graphlet distributions with EMD. PCA finds uncorrelated directions of ordered by contribution to variance, which is sufficient for independence under the assumption of normality [42]. Graphlet distributions for complex networks, however, break this assumption as, although hard to characterize exactly, they are thought to be similar to degree distributions, which, whether or not they are approximately power-law as often claimed [3], are at least often far from normal due to being heavy-tailed [4]. Therefore, we hypothesize that PCA struggles to create an appropriate model for these distributions, which can make the denoising process ineffective, so we also propose an alternative intermediate step that uses Independent Component Analysis (ICA) [5], a method that is explicitly designed to work with heavy-tailed distributions [25].
We test our contributions in clustering tasks involving synthetic and real world networks, following the same experimental procedure as the original NetEmd paper [46], which adheres to the framework of graph comparison introduced by Yaveroglu et al. [50]. We find that our proposed method of denoising using component analysis techniques is able to improve NetEmd's performance on clustering tasks, both in undirected and directed networks. This improvement is particularly significant when comparing networks of different sizes and densities, where NetEmd is already the state of the art.

Graph theoretical concepts
A graph G is composed of a set of vertices (or nodes) V (G) and a set of edges E(G) ⊆ V (G) × V (G), represented by pairs (a, b) ∈ E(G) for a, b ∈ V (G). A graph can be either directed, when the order of the vertices in the pairs expresses direction, meaning that (a, b) ∈ E(G) does not imply that (b, a) ∈ E(G), or undirected otherwise (i.e. when (a, b) ∈ E(G) if and only if (b, a) ∈ E(G)). The size of a graph is the number of vertices in the graph, written as |V (G)|, and the density is the portion of edges present in the graph over the total potential ones (|E(G)|/ |V (G)| 2 in undirected networks). For a directed network, the reciprocity ρ is the probability that, given an edge (a, b) selected uniformly at random from E(G), (b, a) is also a member of E(G). A graph is called simple if it does not contain multiple edges (two or more edges connecting the same pair of vertices) or self-loops (an edge of the form (a, a) that connects a vertex to itself). In this work, we consider only simple graphs.
The neighbourhood of a vertex u ∈ V (G) is defined as N(u) = {v : (v, u) ∈ E(G) ∨ (u, v) ∈ E(G)}. All nodes are assigned consecutive integer numbers starting from 0 and running to |V (G)| − 1. The degree of a vertex is the number of edges it participates in, which is equivalent to the size of the node's neighbourhood. In directed networks, the degree of a vertex can be split into the in-degree and outdegree, counting the number of incoming and outgoing edges, respectively, from that node. In directed networks, the degree of a vertex is the sum of its in-degree and out-degree. In undirected networks, the degree of a node u is simply |N(u)|. The distribution of degrees in a network is called the degree sequence.
FIG. 1. Undirected graphlets of size 2, 3 and 4 and directed graphlets of size 2 and 3. Nodes are numbered according to their orbit, counted by type (directed or undirected) from the size 2 graphlets, and nodes with the same number represent the same orbit.
A subgraph of size k, G k , of a graph G is a graph with k vertices, such that V (G k ) ⊆ V (G) and and is said to be connected when all pairs of vertices have a sequence of edges connecting them. Two graphs G and H are isomorphic, written as G ∼ H, if there is a bijection between V (G) and V (H) such that two vertices are adjacent in G if and only if their correspondent vertices in H are adjacent. A match of a graph H in a larger graph G is a set of nodes that induce the respective subgraph H. In other words, it is a subgraph G k of G that is isomorphic to H. The frequency of a subgraph G k is then the number of different matches of G k in G.
Orbits are unique positions of a graph, calculated by partitioning the set of vertices into equivalence classes where two vertices belong to the same class if there is an automorphism that maps one into the other [37].
Small, connected, non-isomorphic, induced subgraphs are commonly called graphlets [32]. The smallest graphlet considered is a single edge, which can be seen as a size-2 subgraph. An undirected edge has a single orbit (three in the directed case) and the frequency of a node in this orbit is equivalent to the degree of the node. The graphlet degree vector of a node is an extension of the definition of degree to size-k graphs, representing how many times the node occurs in each orbit. The graphlet degree matrix of a graph is the collection of graphlet degree vectors of each node in the graph. Figure 1 shows the graphlets of size 2, 3 and 4 in undirected networks and size 2 and 3 in directed networks, alongside the respective orbits.

NetEmd
NetEmd [46] is a network comparison measure that relies on structural features of the network, mainly the distribution of orbit frequencies. The core idea is formalizing the intuition that the shape of the degree distribution is indicative of the network's generation mechanisms, for instance, a network with a power law degree distribution is generated by a process distinct from a network with a uniform degree distribution. As the graphlet degree vector is a generalization of the degree distribution for graphlets of size k 3, the shapes of the distributions of each orbit also carry information about the topology of the network. Note that because the shape of a distribution is invariant under linear transformations such as translations, using the shape as the focus of the comparison is well-suited to comparing networks of different sizes and densities.
Wegner et al. [46] postulate that "any metric that aims to capture the similarity of shapes should be invariant under linear transformations of its inputs." Thus, they define a measure of similarity between distributions p and q, with non-zero and finite variances, using the Earth Mover's Distance (EMD) [40]: wherep andq are the distributions resulting of scaling p and q to variance 1. Any distance metric d can be used to generate d * ; EMD was used by the authors as it has been shown to be an appropriate metric to compare shapes of distribution in domains such as information retrieval and it produced better results than other distance metrics, like the Kolmogorov or L 1 distances.
Given two networks G and H and a set of m orbits O = {o 1 , o 2 , . . . , o m }, the NetEmd measure is defined as: where p o i (G) and p o i (H) are the distributions of orbit i in graphs G and H, respectively. Note that O may be replaced by any set of network features; in this work we focus on orbits only.

Motivation
In this section, we describe how we couple dimensionality reduction techniques with NetEmd, using them as noise reduction techniques. The rationale is that processes that generate networks are inherently noisy due to the complexity of the systems they represent. This is particularly evident when we consider random networks generated by the same model; we expect them to have similar structure and that similarity to be reflected in the distribution of subgraph and orbit frequencies, but the stochastic nature of the generation process introduces noise in these distributions that makes them harder to identify. Considering the graphlets and orbits shown in Figure 1, it is clear that (for example) directed graphlet G 2 cannot exist if G 0 does not, and similarly directed graphlet G 12 cannot exist if G 1 does not. More generally, it is clear that the different orbits do not represent independent degrees of freedom, and are expected to be constrained by the requirement of combinatorial consistency of the underlying full graph G as well as correlations induced by the generation process.
In practice, both combinatorial enumeration of graphlets and a priori determination of statistical relationships between them (beyond comparison to simple random graphs as in [28]) are computationally infeasible. While this does not matter in practice for the undirected case, for directed graphs as we consider here, the relatively large number and complexity of orbits makes these much more important. We therefore seek an approach that can adjust for the non-independence of orbit counts.
Proceeding somewhat formally, let F G denote the graphlet degree matrix of graph G for a set of orbits O, a n × m matrix where n = |V (G)|, m = |O|, such that [F G ] i, j represent the frequency of orbit j for node i. Further, let f i be the vector of frequencies of node i, i.e. (f i ) j = [F G ] i, j meaning f i is the i-th row of F G . Now note that we are only interested into networks up to isomorphism; in fact, our distance measures depend on empirical histograms of orbit counts, deriving the distribution p o i in (2.1) from histogram heights like which is the proportion of nodes with y counts of orbit o j , where 1 is an indicator function taking the value 1 if its argument is true and the value 0 otherwise. The right-hand side of (3.1) is a sum of random variables, one for each node, and from (2.1) we will further sum over orbits for a function of two such histogram heights. In general, we will expect a random graph model to assign a probability to each graph in some finite set, and according to this measure we expect a probability distribution to be induced on h G = [h o j ,y ]. As discussed above, we expect both combinatorial constraints and correlations between orbit counts, meaning that a sum over all elements of h G as implied by (3.1) and (2.1) will inflate, due to the additivity of variances of random variables, the amount of noise in EMD realisations compared to the amount that is absolutely necessary under the random graph model. Therefore, a natural approach is to denoise using techniques that do not require explicit solution of or simulation from random graph models. In particular, we will use the linear techniques of principal component analysis (PCA) and independent component analysis (ICA), noting their computational efficiency compared to potentially more general nonlinear methods [11]. The main idea of our methodology is therefore to project the orbit frequencies to a lower dimension, training the dimension reduction model to lose minimal information while removing noise.
Since it is not guaranteed that each graph will dimensionally reduce to the same size, the approaches we use allow for expansion of the dimensionally reduced features back to the original (high) dimension. The NetEmd comparison from Equation 2.1 is consequently applied to the reconstructed frequencies. By using the first L components of the linear methods, which contain the most signal from the dataset, to reconstruct the original counts, this can be intuitively understood as decreasing noise within the orbit frequencies when we do not have a good model for that noise.

Principal Component Analysis
Principal component analysis (PCA) is a technique for dimensionality reduction that preserves as much variability of the data as possible, by computing principal components of the data. Principal components are sequences of orthogonal unit vectors that form an orthonormal basis, in which the original dimensions of the data are linearly uncorrelated. When projecting the data onto this new basis, these vectors can be seen as the directions that maximize the variance of the projected data, with the first principal component representing the maximum variance.
Let F G denote the graphlet degree matrix of graph G for a set of orbits O, a n × m matrix where n = |V (G)|, m = |O| and [F G ] i, j represents the frequency of orbit j for node i. We assume that the frequencies have been normalized, a preprocessing step advised before applying PCA (to prevent notation overload, we use the same F G to denote the normalized version of the graphlet degree matrix). The principal components of this matrix are defined by V = F G W, where W is a m × m matrix whose columns are the eigenvectors of F G F G . To allow for only the first L components, the matrix W can be truncated to W L , with dimensions m × L, leading to a transformation V L = F G W L . This truncation is done such that the L eigenvectors chosen to correspond to the largest eigenvalues of W. The original graphlet degree matrix can be reconstructed asF G = V L W L = F G W L W L . The goal of PCA is to learn W L such that the variance of the original data preserved is maximized, while also minimizing the total squared reconstruction error We use the reconstructed graphlet degree matrix to compare networks with NetEmd, keeping Equation 2.1 virtually unchanged from the original formulation: wherep i (G) andp i (H) are obtained fromF G andF H respectively, instead of F G and F H .

3.2.1
Choosing the number of components. Choosing an appropriate number of components to project the data down affects the reconstruction error and therefore the amount of noise reduced. If we take L = m, then clearly we are able to reconstruct the original frequencies perfectly and no noise has been removed from the data. On the other hand, picking a value for L too low, leads to foregoing descriptiveness of the data. A common strategy for picking the number of components is through the amount of variance explained by each component, a strategy that allows us to adapt the number of components to suit the networks we are comparing. The sample covariance matrix of each orbit, i.e. F G i,1 , F G i,2 , . . . , F G i,m , is proportional to F G F G (the derivation for this result can be found in [18], pp. [30][31] and contains the variance of the frequency of each orbit in its diagonal. As a square matrix, according to the spectral theorem, the covariance matrix can be diagonalized by its eigenvectors and the values of the resulting diagonal matrix are the corresponding eigenvalues. These eigenvalues represent variability of each axis in the projected space. Therefore, by taking the sum of the L highest eigenvalues, the ones corresponding to the first L components, we get the variance explained by the first L principal components. The ratio ∑ L 1 λ i / ∑ m 1 λ i , where λ i is the eigenvalue corresponding to the ith eigenvector, measures the proportion of variance explained by the first L principal components. Finally, to condition the number of components on the proportion of variance explained, the smallest L is calculated such that ∑ L 1 λ i / ∑ m 1 λ i r where r × 100% is the percentage of variance explained.

Independent Component Analysis
Independent component analysis (ICA) is a statistical technique in which an observed vector of random variables is thought to be the linear combination of unknown independent components. These components are assumed to be mutually statistically independent and with zero mean. The classic example of an ICA application is a party where a microphone is picking up multiple conversations and the goal is to separate the source signal from each conversation from the mixed one picked up by the microphone.
Transposing ICA to the domain of network comparison, we can draw an analogy from the above application by considering that the voices in the conversation are the nodes in the network and the orbit frequencies are the data points captured by the microphone. Our goal is then to search for the source signals, i.e., network characteristics stemming from its generation mechanism that generate such frequency distribution. We conjecture that by using these source signals to reduce noise in the orbit frequency distributions, our network comparison measure becomes able to more accurately distinguish networks with different generation mechanisms.
As before, let F G denote the graphlet degree matrix of graph G for a set of orbits O and let f i be the vector of frequencies of node i, i.e. f i, j = F G i, j . We can write f i as f i = As i , where A is called the mixing matrix and s i are the L independent components. The goal of ICA is to estimate A and s i from f i only. In practice, algorithms to calculate independent components compute a weight matrix W, a pseudo-inverse of A and obtain the independent components through s i = Wf i . Using the FastICA algorithm [15,16], W is constructed iteratively by finding unit vectors w such that the projection w f i maximizes nongaussianity, measured by the approximation of negentropy where G(x) = log cosh x and ν is a standard Gaussian random variable with mean 0 and variance 1. These unit vectors w are combined into the weight matrix W and decorrelated to prevent convergence to the same maxima by W = W/ WW (the matrix norm is calculated using the l 2 norm); this procedure is then repeated with W = 3 2 W − 1 2 WW W until convergence. Upon obtaining the weight matrix, we calculate its pseudo-inverse to obtain the mixing matrix A and calculateF G similarly to what was done with PCA: we project F G to a lower dimensional space using W and apply the mixing matrix A to retrieve the original data. The final NetEmd formula then takes a similar form as for PCA NetEmd: withp i (G) andp i (H) obtained fromF G andF H respectively, as before.

Directed NetEmd
Directed networks pose challenges unlike the ones observed in undirected networks, due to the combinatorial explosion in number of orbits introduced by distinguishing (u, v) and (v, u) as different edges between u and v. There are 730 orbits when considering directed graphlets of size up 4, compared to 15 in the undirected case, and scaling up to size 5 in directed networks becomes unfeasible as the number of orbits rises to 45, 637. This sharp increase makes the task of counting orbit frequencies for each node even harder, with no known analytical approaches like ORCA [12] that rely on crafting sets of equations that exploit combinatorial relationships between smaller graphlets to compute orbit counts. Instead, enumeration-based approaches are required in the directed case, which are known to be at least an order of magnitude slower than analytical approaches in undirected networks. In order to adapt NetEmd to directed networks, we use the G-Trie [38] data structure and the counting algorithm proposed by Aparicio et al. [2], publicly available at [36], modified to return the graphlet degree matrix instead of the graphlet degree distribution. How to handle this increased number of orbits in the comparison is also a challenging problem, especially when comparing networks that contain only a small subset of orbits and differences within those orbits can get diluted when taking the average over the whole set. Aparicio et al. [2] argue that the comparison should only be done over the orbits present in at least one of the networks, using networks G and H in Figure 2 as an example. In this case, the authors show that using the original formulation of the Graphlet Degree Distribution Agreement (GDA) (described in Section 5.4.2) inflates the similarity score between the two networks far beyond what would be expected from two so distinct looking networks, with an agreement score of 0.92 against 0.32 when using their modified version that compares only orbits in at least one of the networks. This modification to GDA is well-founded because GDA is a measure of agreement, bounded between 0 and 1 and meant to be interpretable; a score close to 1 means that the networks are similar and a score close to 0 the opposite, so the magnitude of the value that the measure outputs is of interest. However, NetEmd is a measure of distance between networks and the relative difference between distances is more informative than the absolute value of that distance, so including information that an orbit is missing in both networks improves our ability to tell if they are less distant from each other when comparing against another pair of networks.
With the above in mind, we propose two versions of directed NetEmd. Given two networks G and all the orbits in graphlets of size up to 3 or 4 and the formulation is the same as the original NetEmd in Equation 2.1. The second is using the same idea of Aparicio et al. [2], which we refer to as weighted NetEmd, where we restrict the set of orbits to those that occur in at least one of the networks G or H.
Graphs G and H in Figure 2 show an analogous example to the one shown by Aparicio et al. [2] to demonstrate the difference between weighted and original NetEmd applied to directed networks. Although it is usually impractical to scale up to size 5 in directed networks, in a small example like this one it becomes feasible to do so. In this context, we calculate a distance of 5 × 10 −5 using all orbits and 0.12 when using the weighted version.

Experimental Setup
Experiments were performed on an AMD Opteron Processor 6380 with 1.4 GHz and 2 MB of cache memory, using Python version 3.5.4. The implementation of our methods is available online at https: //github.com/migueleps/denoise-dir-netemd.

Measure of Cluster Performance
Given a set of networks G = {G 1 , G 2 , . . . , G N }, divided in disjoint classes C = {c 1 , c 2 , . . . , c m }, we use the performance measure proposed by Wegner et al. [46] that captures the idea that networks of the same class should be nearer to each other than to networks of other classes. It is defined as an empirical probability P(G) = d(G, G 1 ) < d(G, G 2 ), where G 1 is a network selected randomly from the same class as G and G 2 is randomly selected from a different class, and d is the network comparison statistic. The performance over the whole dataset is the average P(G) over all the networks in G , and we write this as P = P(G ) = 1 |G | ∑ G∈G P(G). In Appendix D, we present results using area under precision-recall curve (AUPR) and adjusted Rand index (ARI).

Synthetic datasets
We reproduce the experimental setup of Wegner et al. [46], testing our proposed modifications to NetEmd on the same eight random network models. These models are Erdős-Rényi [7], Barabasi-Albert preferential attachment [3], configuration model, geometric random graphs [8], geometric gene duplication model [10], duplication divergence of Vazquez et al. [44], duplication divergence of Ispolatov et al. [17] and Watts-Strogatz [45]. Details of parameters for each model are described in Appendix E.
We generate 10 networks per model per combination of number of nodes and average degree. The set of number of nodes used is N ∈ {1250, 2500, 5000, 10000} and the set of average degrees is k ∈ {10, 20, 40, 80}. This leads to a total of 160 networks per model and 1280 networks in total.
For the directed datasets, we add varying levels of reciprocity ρ ∈ {0, 0.25, 0.5, 0.75, 1}. To generate a directed network with reciprocity ρ, we take the undirected version, we duplicate and invert all edges (if (u, v) is in the undirected network, we add (v, u) to the directed version) and we then take a proportion 1 − ρ of these pairs and choose randomly a direction to remove (either (u, v) or (v, u)). This is done incrementally such that if (u, v) is removed from the dataset with 75% reciprocity, then it is also removed from the datasets with 50%, 25% and 0% reciprocity.
These datasets of synthetic networks give rise to two tasks aimed at gauging how well a network comparison measure separates clusters according to their generation mechanism, with the ground truth given by the random network model used to generate the networks in each cluster. The first, more simple task is to separate networks with the same number of nodes and the same density. To this end, we create 16 groups based on the combinations of N and k, each group containing 10 realizations of each random network model, for a total of 80 networks. This task is equivalent to RG 1 in the original NetEmd paper [46]. The second task, equivalent to RG 3 in the original NetEmd paper, is to compare all 1280 networks simultaneously, finding the 8 clusters of 160 networks. This task measures how sensitive the comparison measure is to differences in an order of magnitude in number of nodes and edges, making the separation by model type more difficult. In directed networks, we repeat Task 1 and 2 for each level of reciprocity, determining how sensitive the network comparison is to this third parameter that impacts the set of orbits available for comparison (the set of orbits with 0% reciprocity is disjoint from the set of orbits with 100% reciprocity).

Real world datasets
We use the dataset of Onnela et al. [30] to validate our methodology on a mix of real world and synthetic networks. The multiple sources that make up this dataset lead to a heterogeneous set of networks; however, the ability to separate these networks according to their domain is desired from a network comparison method. We were unable to find sources for the original 746 networks, so we use the reduced set proposed by Ali et al. [1] with 151 unweighted and undirected networks. There is no ground truth for the true clusters in this dataset and the dendograms constructed through the methods proposed by Onnela et al. [30] and Ali et al. [1] disagree on the composition of each cluster. Therefore, we aim to reconstruct clusters according to the type of data, which can be visualized in the supplementary material of Ali et al. [1].
For directed networks, we download four datasets from the SNAP library [23]: Gnutella peer-to-peer file sharing network from August 2002 [22,39] (9 networks with an average reciprocity of 0%), CAIDA autonomous systems relationships datasets from January 2004 to November 2007 [21] (122 networks with an average reciprocity of 100%), ego networks of circles from Google+ [27] (132 networks with an average reciprocity of 28%) and ego networks of lists from Twitter [27] (973 networks with an average reciprocity of 54%).

5.4
Other network comparison methods 5.4.1 GCD. Yaveroglu et al. [49] note that orbit counts have dependencies between them, making them redundant since they can be expressed as a linear combination of other orbits. The authors identify 11 out of 15 orbits in graphlets of size up to 4 and 56 out of 73 orbits in graphlets of size up to 5 as non-redundant. The authors construct the graphlet degree vector of each node in the network using this reduced number of orbits (although the full set of orbits may also be used), making up a matrix where each row is the graphlet degree vector of each node, and compute the Spearman's correlation coefficient between each pair of orbits, i.e., each pair of columns in this matrix. The pairwise Spearman's correlation coefficients are aggregated in a square m × m matrix called the Graphlet Correlation Matrix (GCM), where m is the number of orbits used (e.g., 11 if using non-redundant orbits of graphlets up to size 4).
To compare two networks, the authors propose the Graphlet Correlation Degree (GCD), which is the Euclidean distance between the upper triangle of the GCM of each network.
Sarajlić et al. [41] extend GCD to directed networks, finding 13 non-redundant orbits in graphlets of size up to 3 and 129 in graphlets of size up to 4.
To compare the GDD distributions of the same orbit in different networks G and H, the authors propose the GDD-agreement (GDA) metric defined as: As a measure of agreement, the output of GDA is 1 if the distributions are identical, as opposed to NetEmd which is a measure of distance meaning it will output 0 for identical distributions. To aggregate the agreements of the multiple orbits, the authors propose either a arithmetic or a geometric mean. We compare against the implementation of Aparicio et al. [2], which uses the arithmetic mean and is available at [36]. This implementation uses the idea that if a graphlet has a frequency of 0 in both networks, then the orbits from that graphlet are excluded from the arithmetic mean, for both in the directed and undirected versions.

Undirected Results
We present the results for Tasks 1 and 2 in undirected networks and the results for the Onnela et al. [30] dataset in Table 1 We find that our proposed modifications to NetEmd achieve an increase in performance, with ICA NetEmd outperforming original NetEmd in synthetic datasets and PCA NetEmd in the Onnela et al. dataset. The performance gain in Task 1 is not significant as the results are within a standard error of each other, but on Task 2, the more difficult task, the difference is more pronounced. This

Directed Results
We show the results for Task 1 in directed networks in Table A.3 and the results for Task 2 and the dataset of real world directed networks in Table 2. The results for Task 1 show the average and standard error for the 16 values of P, one for each combination of N and k, for each level of reciprocity. The results for Task 2 and real world directed networks dataset show the single value of P after comparing the 1280, for each level of reciprocity, and 1232 networks, respectively. Similarly to the undirected case, we find that the results for each algorithm in Task 1 are similar across different levels of reciprocity, meaning that the various different versions of NetEmd and DGCD are able to distinguish the generation mechanism regardless of how many edges are reciprocated in a directed network. Results show that DGCD with 129 orbits achieves best performance for 0%, 25%, 50% and 75% and PCA NetEmd with 90% explained variance for 100% reciprocity. We observe no significant difference between using NetEmd with all orbits, with weighted orbits or coupled with dimen-Parameter  Table 2. Results for Task 2 in directed networks, for each level of reciprocity, and for the real world networks dataset. The metric used for this task is the sole value of P after comparing the 1280 and 1231 networks, respectively. The parameter column indicates: graphlet size used in directed and weighted NetEmd; percentage of variance explained to determine the number of components in PCA NetEmd (using orbits in graphlets of size up to 4); the number of components used in ICA NetEmd, using orbits in graphlets of size up to 4; the number of orbits used by DGCD; graphlet sizes used in GDA. The bolded values are the maximum for each dataset. Note that results for GDA with graphlet size 4 took longer than a week to return results for Task 2 and in the real world networks dataset, at which point we stopped the computation. sionality reduction techniques. Results with TriadEMD are non-significantly better than our NetEmd version with orbits from size 3 graphlets in this task, indicating that using graphlets of size 2 for the comparison in this task does not help differentiating between different models. This is likely due to the configuration model, which shares the same degree distribution as the duplication divergence model of Vazquez et al. [44]. In Task 2, our results show a noticeable gain in performance when using ICA NetEmd over other NetEmd versions and DGCD, in particular when using only 2 components and in reciprocity levels smaller than 100%. We find that using orbits in graphlets of only size 2 and 3 orbits achieves better results for this task than using ones with size up to 4, with the exception of the two extremities in reciprocity (0% and 100%), where having more available orbits for comparison helps performance. In this task, we also find that using Weighted NetEmd degrades performance compared to using all orbits. Results also show that our versions of NetEmd perform very similarly to TriadEMD, with a slight advantage to TriadEMD at 0, 75% and 100% reciprocity, which is perhaps expected again due to the configuration model adding a confounding factor when using graphlets G 0 and G 1 from Figure 1 in the comparison.
On the other hand, Weighted NetEmd reaches the best performance in the real directed networks dataset, where, similarly to Task 2 in synthetic networks, using smaller graphlets yields better results. Table 2 shows that ICA NetEmd with 10 and 15 components obtains better performance than NetEmd without dimensionality reduction and with size 4 graphlets. This comparison is unfair to ICA NetEmd, as results suggest that when there is a large discrepancy between the sizes of the networks, using smaller graphlet sizes improves the performance of NetEmd. Therefore, although Weighted NetEmd achieves the best performance out of the parameters we include in Table 2, the top performance in the real world directed networks dataset is achieved by ICA NetEmd with size 3 graphlets and 10 components. All four versions of NetEmd that we propose achieve better results than DGCD and GDA in this dataset. This is also true when comparing our versions with size 3 graphlets against TriadEMD, indicating that in real world networks the degree distribution is able to contribute in a positive manner towards being able to discriminate networks of different sources

Performance by number of ICA components
We examine the impact on performance of varying the number of components used for ICA NetEmd, with the caveat that increasing the number of components leads to a greater computational complexity as the FastICA algorithm takes more iterations to converge. We also explore the performance of ICA NetEmd under smaller graphlet sizes, namely graphlets sizes 3 and 4 for undirected networks and size 3 for directed. There are two different ways of using smaller graphlets on ICA NetEmd, the more natural way in which we compute the independent components using the set of orbits related to the desired graphlet size, or an alternative way in which the components are computed using a larger graphlet size and the reconstructed graphlet degree matrix is posteriorly truncated to include only the orbits of the desired graphlet size. For example, when calculating network distances using undirected graphlets of size up to 4, there are 15 orbits. In the first method, the graphlet degree matrix F G is a |V (G)| × 15 matrix and we simply calculateF G as described in Section 3, constraining the number of components used to a maximum of 14 components. In the second method, F G is a |V (G)| × 73 matrix, allowing a maximum of 72 components to computeF G , which is then truncated to 15 columns. Figures 3 to 8 show, in order, performance with varying number of components in undirected networks Task 1, undirected networks Task 2, directed networks Task 1, directed networks Task 2, the Onnela et al. dataset and directed real world networks.
We find that adding more components to ICA does not always translate to a better performance, instead we find that there is a performance maximum obtained with a low number of components that depends on the task. Adding more components after this maximum leads to a decline in performance, either an immediately visible decline (Figures 4, 6 and 7) or after plateauing at that maximum value (Figures 3, 5 and 8). This behaviour indicates that our ICA algorithm is working as expected, as we find a suitable number of components that maximizes the noise reduced in the orbit frequency distributions, so adding more components reduces the difference between F G andF G which negates the positive impact on that noise reduction. Figures 5 and 6 show a sharp decline in performance when using graphlets of size up to 3 with reciprocity 0% and 100%, after 13 and 4 components respectively. This is caused by these extreme levels of reciprocity limiting the orbits that can occur in the networks we compare. With 0% reciprocity, only graphlets with no reciprocal edges can occur in these networks, so all orbits in graphlets with reciprocal edges have a frequency of 0. As there are only 13 orbits belonging to graphlets with no reciprocal edges, when we try to compute more than 13 independent components the performance degrades as more components are added. With 100% reciprocity, we find a similar situation but there are only 4 orbits in graphlets of size 2 and 3 for which all edges are reciprocal (orbits 2, 27, 28 and 32 in Figure 1).

Conclusion
We presented two extensions of NetEmd to directed networks, one a direct extension that uses all orbits of graphlets up to size 4 and another that compares networks based only on the orbits that exist in at least one of the networks. We showed that our methodology achieves state of the art performance in large datasets of synthetic networks with heterogeneous network sizes and in datasets of real world networks.
We also proposed to add dimensionality reduction techniques, namely PCA and ICA, as a preprocessing step to the network comparison. The goal of this step is to use dimensionality reduction as a process to attenuate noise within orbit frequencies, that can be introduced by random number generation (in synthetic networks) or data collection/representation (in real world networks). Results show that this preprocessing improves the performance of NetEmd not only when comparing networks with the same number of nodes and edges, but especially in large datasets containing networks of different sizes.
From an end-user perspective, our extensive testing allows us to recommend guidelines for which version of NetEmd to use in which situation. If the set of networks is homogeneous in number of nodes and average degree, using the largest graphlet size available (size 5 in undirected, size 4 in directed) coupled with ICA (10 to 15 components) leads to the best results. In the cases where there is a large variety in number of nodes or average degree between the networks in the dataset, using smaller graphlets is generally advised. The exception is in directed networks, if the connections within the networks have no reciprocity, then using graphlets of size 4 (with ICA) is more likely to yield a more accurate comparison, due to having a larger set of orbits to with non-zero frequency. On the other hand, if the average reciprocity allows using the whole range of graphlets of size 3 and 4, then it is more likely that using graphlets of size 3 coupled with ICA will lead to a more accurate distinction between networks. In this case, 2 components might be enough if the size of the networks are all within one order of magnitude from each other but 5 to 10 components if there is a wider difference.

MEPS is funded by Engineering and Physical Sciences Research Council
Manchester Centre for Doctoral Training in Computer Science (grant number EP/I028099/1). TH is supported by the UKRI through the JUNIPER modelling consortium (grant no. MR/V038613/1) and the Engineering and Physical Sciences COVID-19 scheme (grant number EP/V027468/1), the Royal Society (grant number INF/R2/180067), and the Alan Turing Institute for Data Science and Artificial Intelligence.

A. Discussion on the time complexity of NetEmd
Calculating the NetEmd measure is a process that can be separated in two phases: acquiring the distributions of the network statistics of interest and comparing the statistics using the EMD. In the case of the algorithm we propose in this work, obtaining the distributions of network statistics involves calculating the graphlet degree matrix and performing PCA or ICA on this matrix.
Obtaining the graphlet degree matrix is the most computationally expensive step of this process, with a time complexity of O(Nd m−1 ), where N is the number of nodes in the network, d the maximum degree of any node in the network and m the size of the graphlets being enumerated. In undirected networks, we use the combinatorial algorithm ORCA [12] that relies on an analytical approach to set up a system of linear equations that relates different orbit frequencies. In directed networks, no such approaches are known [37], so we rely on G-Tries [2,38], a data structure that supports a graphlet enumeration algorithm by representing subgraphs in a prefix tree, the state of the art in enumerating directed graphlets, to our knowledge [2,37].
The complexity of PCA is split in two parts, computing the covariance matrix is O(N p 2 ) and doing the eigenvalue decomposition is O(p 3 ), where N is the number of nodes in the network and p the number of orbits being considered. This leads to an overall complexity of O (N p 2 + p 3 ).
The complexity of ICA is harder to characterize. Firstly, the FastICA algorithm assumes that the data has been centered and whitened [16]. The implementation we use, from scikit-learn [31], uses PCA to do the whitening preprocessing, so its complexity is at least as high as PCA. The iterative algorithm to find the weight matrix W is repeated until the matrix has converged, but there is no prior assurance that the algorithm converges or how many iterations it takes to do so. The scikit-learn implementation defines a maximum number of iterations to stop execution, we set this value to 1000. Each iteration involves calculating log cosh(w f i ) for each component and each node, which the scikit-learn implementation does in O(Nc 2 ), where c is the number of components to be calculated, and the decorrelation step, which is done in O(c 2 ). The overall complexity of ICA then becomes O(N p 2 + p 3 + INc 2 ), where I is the number of iterations.
Wegner et al. [46] calculate the complexity of comparing two graphlet distributions using EMD * to be O(k(N + N ) log(N + N )), where N and N are the number of nodes of each network and k is the maximum number of function calls to the optimization algorithm used to align the distributions.

C. Discussion on orbits from smaller graphlet sizes
Wegner et al. [46] have shown that reducing the graphlet size as input to NetEmd can improve the quality of the clusters under certain conditions. A more common reason to use smaller graphlets for comparison is a computational load argument, as we show in Appendix A, increasing the graphlet size leads to more computation time in each step of the NetEmd framework. In the first step, computing orbit frequencies of graphlets of size 4 in directed networks or size 5 in undirected networks can be prohibitive for very large graphs. Another consideration is the time taken to compute the EMD between orbit distributions, particularly in directed networks with graphlets of size 4 requiring 22 times more calls to the EMD * function than using size 3. Finally, computing the principal or independent components also contain a complexity component dependent on the number of orbits we consider as features to these methods, which also contribute to greater computational load when increasing the graphlet sizes as input to NetEmd. Therefore, it is important to understand in which situations decreasing the graphlet size as input to NetEmd leads to similar or better performance to avoid paying the high computational cost.
We previously presented the performance of ICA NetEmd under smaller graphlet sizes, for multiple components used, in Figures 3 to 8. The results for Task 1, both in directed and undirected networks, demonstrate that when comparing networks of similar sizes and densities using bigger graphlet sizes yields a more accurate comparison. The same holds for Task 2 in the undirected case and in the extreme cases of reciprocity (0% and 100%) in Task 2 of the directed experiments, where ICA achieves the top performance with size 5 undirected and size 4 directed graphlets, respectively. On the other hand, in the datasets of real world networks and in directed Task 2 with reciprocity of 25%, 50% and 75%, size 5 graphlets in undirected networks and size 4 in directed networks perform worse than size 4 and size 3, respectively. For the datasets of real networks in particular, this difference in performance can be explained by noticing that these datasets contain very small networks such that some graphlets have no matches in them. This issue is amplified the larger the graphlet size considered is and it adds a confounding factor to the network comparison, similar to the example of Figure 2. Incidentally, this is the only dataset where Weighted NetEmd achieves better performance than following the original formulation.   Figures 3 to 8 also show the performance results of using the full set of size k orbits to denoise the distributions of k − 1 or k − 2 orbits. In this second way of using ICA, the only time saved compared to using the full set of orbits is the time to compute the EMD between orbit distributions and we find that in no case it leads to the best clustering performance. However, Figure 4 shows a breakpoint at 30 components, where the performance using size 3 and 4 orbits after denoising with size 5 succeeds in performing better than size 5. The reason for this is that, as more components are added, the smaller the reconstruction error ||F G −F G || 2 is, so the performance with those parameters tends to approximate the performance of the original NetEmd with those graphlets sizes.
The performance of PCA NetEmd under smaller graphlet sizes follows a similar set of rules to ICA NetEmd, with better performances in Task 1 measured when using larger graphlets and in Task 2 and real world networks when using smaller graphlet sizes. Tables A.4 to A.6 detail the results for multiple combinations of graphlet size and percentage of variance explained. Unlike ICA NetEmd, our experimental setup does not allow us to claim that after a certain value of explained variance performance degrades or stabilizes. Instead, results suggest that the optimal threshold depends on both the dataset and the size of graphlet chosen. Like with ICA NetEmd, choosing a low value for percentage of explained variance for Task 2 in directed networks yields the best performance for 0, 25, 50 and 75% reciprocity, but otherwise choosing 90 or 95% explained variance seems to be a good rule of thumb for best performance with PCA NetEmd.   Table A.6. Results for PCA NetEmd with different combinations of graphlet size and percentage of variance explained in Task 2 in directed networks, for each level of reciprocity, and for the real world networks dataset. The metric used for evaluation is the mean (and standard error of the mean) of the 16 values for P, in each combination of number of nodes, network density and reciprocity. The bolded values are the maximum for each reciprocity level.

D.1 Area Under Precision-Recall (AUPR)
Sarajlic et al. [41] evaluate the performance of their graph comparison tool using the Area Under Precision -Recall curve (AUPR), after the framework proposed by Yaveroglu et al. [50]. The metric is calculated as follows: for each value of a parameter ε 0, if the distance between two networks is less than ε, then the networks belong to the same cluster. The parameter ε ranges from 0 to 1, with values incrementing by 5 × 10 −3 . For each value of ε, we compute the number of true positives, false positives, true negatives and false negatives. In this context, we define these concepts as: • True positive (TP): the distance between the networks is smaller than ε and they were generated by the same random network model.
• False positive (FP): the distance between the networks is smaller than ε but they were generated by different random network models.  Table A The bolded values are the maximum for each task. Italic values are within one standard error of the maximum performance. Note that results for GDA with graphlet size 5 took longer than a week to return results for Task 2, at which point we stopped the computation.
• True negative (TN): the distance between the networks is greater than ε and they were generated by different random network models.
• False negative (FN): the distance between the networks is greater than ε but they were generated by the same random network model.
From these quantities, we calculate precision as T P T P+FP and recall as T P T P+FN , obtaining 200 tuples of (precision, recall), from which we are able to plot the precision-recall curve. The area under this curve is a metric relevant to our problem as it puts an emphasis on the positive predictive value of the model, disregarding the true negatives that compose the majority of datasets with imbalanced labels such as this one (with 8 random network models, positive labels comprise only 12.5% of the dataset). Table A.7 shows the AUPR for Task 1 and Task 2 in the undirected case and the Onnela et al. dataset.

D.1.1 Undirected Results
In Task 1, according to AUPR, using PCA or ICA does not lead to a gain in performance, unlike what we observed when using P. This is not surprising as the improvement we displayed with P was  Table A.8. Results for Task 2 in directed networks, for each level of reciprocity, and for the real world networks dataset. The metric used for this task is the sole value of AUPR after comparing the 1280 and 1231 networks, respectively. The parameter column indicates: graphlet size used in directed and weighted NetEmd; percentage of variance explained to determine the number of components in PCA NetEmd (using orbits in graphlets of size up to 4); the number of components used in ICA NetEmd, using orbits in graphlets of size up to 4; the number of orbits used by DGCD; graphlet sizes used in GDA. The bolded values are the maximum for each dataset. Note that results for GDA with graphlet size 4 took longer than a week to return results for Task 2 and in the real world networks dataset, at which point we stopped the computation.
within a standard error and both methodologies are able to separate networks between these random models with a very high degree of accuracy. Task 1 serves as a proof of concept that our proposed measure is able to perform the most simple task. The results for Task 2 are aligned with what we observe when using P, coupling ICA with NetEmd leads to the best separation between models when different network sizes and densities are present. In the Onnela et al. dataset, the results with AUPR also agree with P, both performance metrics consider the original NetEmd with size 4 orbits to be one of the top performers with this dataset. The difference between the two metrics is that when measuring with P, we observe that using PCA with 90% explained variance leads to similar performance (Table A.4), but the top performer in this dataset according to AUPR is PCA with 99% explained variance, achieving an AUPR score of 0.791.

D.1.2 Directed Results
The results using AUPR for Task 1 in directed networks are shown in Table  A.9 and for Task 2 and the dataset of real world directed networks in Table A. 8. In Task 1, the main difference we find is that GDA with size 3 orbits shows an improvement over the other algorithms, becoming the top performer with 25% (tied with DGCD with 129 orbits) and 75%   reciprocity. DGCD with 129 orbits still shows the best performance with 0% reciprocity, but PCA with 90, 95 and 99% explained variance is within one standard error, similarly to the results measured with P. At a 100% reciprocity, NetEmd with size 4 orbits achieves the best performance, but PCA with 90, 95 and 99% explained variance and ICA with 2 components obtain similar AUPR scores. As before, we highlight that the high scores achieved by the 6 algorithms in Table A.9 serve as a proof of concept that they are able to distinguish the random network models for different levels of reciprocity, even if the AUPR score is not within a standard error of the best score. In Task 2, similarly to the undirected version, we find an agreement between the AUPR score and P, with ICA with 2 components highlighted as the highest performing algorithm for 0, 25, 50 and 75% reciprocity datasets and NetEmd with size 4 orbits for 100% reciprocity. The same agreement between AUPR and P is found on the dataset of real world directed networks, with the top performer being ICA with 6 components with an AUPR of 0.891 but not highlighted in Table A. 8.
The conclusions regarding the difference between our proposed NetEmd methods and TriadEMD also hold when using AUPR instead of P; in the synthetic datasets TriadEMD achieves better AUPR scores than both NetEmd and Weighted NetEmd with orbits from size up to 3 graphlets, but worse results than NetEmd with size 4 graphlets. In the real world directed networks dataset, both NetEmd and Weighted NetEmd with size 3 orbits obtain better AUPR scores than TriadEMD.

D.2 Adjusted Rand Index
The adjusted Rand index [14,34] is a metric to evaluate the similarity of data clusterings; a corrected for chance version of the Rand index [34] was proposed by Hubert and Arabie [14]. Let G = {G 1 , G 2 , . . . , G N } be the set of networks we wish to partition into clusters. For the synthetic datasets, the ground truth partition of this set is known, we have 8 models of random graphs and we evaluate a network comparison method by its ability to reconstruct these groups. For the Onnela et al. [30] dataset, the ground truth partition used is the one shown in the supplementary material of [1]. Let P * = {P * 1 , P * 2 , . . . , P * c } be the ground truth partitions and P = {P 1 , P 2 , . . . , P c } be the clusters produced by the network comparison algorithms. These clusters are generated by applying a hierarchical clustering algorithm to the matrix of pair-wise distances calculated by the network comparison; in our case we used the implementation by scikit-learn [31], the AgglomerativeClustering class with complete linkage.
Given the two set of clusters, we define the following four quantities: • a: the number of pairs of elements of G that are in the same cluster in P * and in the same cluster in P.
• b: the number of pairs of elements of G that are in different clusters in P * and also in different clusters in P.
• c: the number of pairs of elements of G that are in the same cluster in P * but in different clusters in P.
• d: the number of pairs of elements of G that are in different clusters in P * but in the same cluster in P.
These quantities can be interpreted as the number of true positives, true negatives, false negatives and false positives, respectively, if we consider the attribution of a pair of networks to a cluster as a decision problem. The Rand index is then defined as: The general form for correcting a index for chance is: Hubert and Arabie [14] show how to calculate E(RI), which we omit for brevity, leading to the formula for the adjusted Rand index (assuming a maximum value of 1 for the Rand index when c and d are 0): where n i j is the number of networks in common between P * i and P j , i.e., n i j = |P * i ∩ P j |, a i = ∑ c j=1 n i j and b j = ∑ c i=1 n i j . Table A.10 shows the ARI for Task 1, Task 2 and real world networks in the undirected case.

D.2.1 Undirected Results
As with AUPR, we find that using ARI as a performance metric does not change the conclusions drawn from the results with P. Recalling that the ARI metric measures the quality of hierarchical clusters generated from a distance matrix, for Task 1, the output of original NetEmd with size 5 orbits leads the best score using this metric. More evidence of the quality of the clusters produced by the original NetEmd with size 5 orbits is to notice that increasing the number of components in ICA or the percentage of variance explained in PCA (and therefore approximating the results of original NetEmd) leads to higher ARI score.
For Task 2, the gain measured by P and AUPR when using ICA NetEmd with 2 components is also reflected in the ARI score. For the Onnela et al. dataset, the best performance we measure not included in Table A.10 is ICA NetEmd with size 4 graphlets and 8 components, with the components calculated with size 5 graphlets, with a score of 0.713.

D.2.2 Directed Results
The results using ARI for Task 1 in directed networks are shown in Table  A.12 and for Task 2 and the dataset of real world directed networks in Table A. 11.
In Task 1, we find that, according to ARI, the performance of all NetEmd variants is significantly below the performance of DGCD with 129 orbits for 0, 25 and 50% reciprocity and of GDA with size 3 graphlets for 75% reciprocity. At 100% reciprocity, NetEmd, PCA NetEmd with 80% variance and GDA, all with size 4 graphlets, share the top performance, but with a smaller difference to other algorithms and parameters. These results are a stark contrast to the ones we observe using AUPR and P, where scores for this task did not differ significantly between all the algorithms.
In Task 2, unlike in the other metrics, ICA NetEmd with 2 components is no longer considered the best performer for 0% reciprocity, with TriadEMD achieving the best result instead. On the other hand, according to ARI, ICA NetEmd with 2 components is the top performer for 100% reciprocity. For 25, 50 and 75% reciprocity, ICA NetEmd with 2 components performs significantly worse when comparing to the results using AUPR and P. ARI scores indicate that the best performers in these levels of reciprocity are PCA NetEmd with 99% variance, NetEmd with size 3 graphlets and PCA NetEmd with 95% variance. As we mention previously, using a smaller graphlet size usually leads to better performance when different network sizes and densities are involved. We find partial agreement to this claim when using ARI as the performance metric. At 0% reciprocity, we find that ICA NetEmd with 2 components and size 3 orbits achieves an ARI of 0.509 (higher than the score of ICA NetEmd with 2 components and size 4 orbits in Table A.11, but lower than the score of TriadEMD); at 75% reciprocity, PCA NetEmd with 90% variance explained and size 3 orbits achieves an ARI of 0.495 (higher than the score of PCA NetEmd with 95% variance explained and size 4 orbits in Table A.11).
For the dataset of real world directed networks, ARI disagrees with AUPR and P, assigning the best performance to GDA with size 3 orbits. This task also favours using smaller graphlets for the comparison, so the PCA NetEmd and ICA NetEmd results are at a disadvantage compared to the other methods as we only report results with size 4 in Table A.11. For these two algorithms with size 3 graphlets, the best ARI scores on this task are 0.524, when using PCA NetEmd with 90% explained variance, and 0.489, when using ICA NetEmd with 19 components.
On the overall, in spite of some agreements between ARI and AUPR or P, we find ARI to be an unreliable metric to gauge the performance of network distance measures, as it relies on an intermediate Parameter  Table A.11. Results for Task 2 in directed networks, for each level of reciprocity, and for the real world networks dataset. The metric used for this task is the sole value of ARI after comparing the 1280 and 1231 networks, respectively. The parameter column indicates: graphlet size used in directed and weighted NetEmd; percentage of variance explained to determine the number of components in PCA NetEmd (using orbits in graphlets of size up to 4); the number of components used in ICA NetEmd, using orbits in graphlets of size up to 4; the number of orbits used by DGCD; graphlet sizes used in GDA. The bolded values are the maximum for each dataset. Note that results for GDA with graphlet size 4 took longer than a week to return results for Task 2 and in the real world networks dataset, at which point we stopped the computation.
step of clustering that is unrelated to the measures themselves. Using a different clustering method or even a different linkage leads to disparate results compared to the ones we present in this section, which adds a confounding factor when measuring the performance of the network comparison itself. An example of this unreliability is the ARI score of 0 given to GDA with size 3 orbits and to TriadEMD at 100% reciprocity, when the other metrics showed a performance comparable to DGCD.
• Erdős-Rényi (ER) model [7]. A random graph with n is generated by picking m unique edges chosen at random from the n 2 possible edges. m was chosen to generate networks with the  Table A.12. Results for Task 1 in directed networks, for each level of reciprocity. The metric used for evaluation is the mean (and standard error of the mean) of the 16 values for ARI, in each combination of number of nodes, network density and reciprocity. The parameter column indicates: graphlet size used in directed and weighted NetEmd; percentage of variance explained to determine the number of components in PCA NetEmd (using orbits in graphlets of size up to 4); the number of components used in ICA NetEmd, using orbits in graphlets of size up to 4; the number of orbits used by DGCD; graphlet sizes used in GDA. The bolded values are the maximum for each reciprocity level. appropriate density.
• Barabasi-Albert (BA) preferential attachment model [3]. An initial graph is created with m nodes and new nodes are added iteratively to the network, each new node is connected to m existing nodes, picked randomly with probability proportional to their degree. m was chosen to generate networks with the appropriate density.
• Geometric random graphs [8], using a 3-dimensional square (D = 3). Nodes are randomly embedded in a D-dimensional space and are connected if the Euclidean distance between them is smaller than a threshold r. This threshold parameter was determined by grid search to generate networks with the desired number of edges.
• Geometric gene duplication model [10]. Starting with an initial network of 5 nodes embedded in a 3-dimensional space, on each iteration a node is chosen randomly to be duplicated. The duplicated node is placed randomly in the 3-dimensional space at an Euclidean distance of at most 2 from the original node. This process is repeated until the desired number of nodes, and nodes at a distance of r or less at connected. This distance r is chosen to generate a network with the appropriate number of edges.
• Duplication divergence model of Vazquez et al. [44]. The network grows in two stages. In the first stage, a node v is chosen randomly to be duplicated into a node v , that keeps all the edges from v. The nodes v and v are connected with probability 0.05. In the second stage, one of the duplicated edges ((v, u) or (v , u)) is chosen randomly and deleted with probability q. This process is repeated until the network grows to the desired number of nodes and q is chosen to generate networks with the appropriate density.
• Duplication divergence model of Ispolatov et al. [17]. Starting from a network with a single edge, a node is chosen randomly and duplicated, with the new node keeping each of the original's neighbours with probability p. This parameter is chosen to generate networks with the desired number of edges.
• Configuration model, using the Duplication divergence model of Vazquez et al. [44] as the generator for the graphical degree sequence.
• Watts-Strogatz (WS) model [45]. Nodes are placed in a ring and connected to their k nearest neighbours in both sides of the ring. Each edge is then rewired with probability 0.05 to a new node selected at random. The parameter q is chosen to generate networks with the appropriate average degree.

E.2 Real world networks
The summary statistics of the Onnela et al. [30] and the real world directed networks are given in Table A.13.