Title: SPECTRAL GRAPH WAVELETS FOR STRUCTURAL ROLE SIMILARITY IN NETWORKS Anonymous

Abstract: Nodes residing in different parts of a graph can have similar structural roles within their local network topology. The identification of such roles provides key insight into the organization of networks and can also be used to inform machine learning on graphs. However, learning structural representations of nodes is a challenging unsupervised-learning task, which typically involves manually specifying and tailoring topological features for each node. Here we develop GRAPHWAVE, a method that represents each node's local network neighborhood via a low-dimensional embedding by leveraging spectral graph wavelet diffusion patterns. We prove that nodes with similar local network neighborhoods will have similar GRAPHWAVE embeddings even though these nodes may reside in very different parts of the network. Our method scales linearly with the number of edges and does not require any hand-tailoring of topological features. We evaluate performance on both synthetic and real-world datasets, obtaining improvements of up to 71% over state-of-the-art baselines.

1. INTRODUCTION:

Structural role discovery in graphs focuses on identifying nodes which have topologically similar local neighborhoods (i.e., similar local structural roles) while residing in potentially distant areas of the network (Figure 1 ). Such alternative definition of node similarity is very different than more traditional notions (Perozzi et al., 2014; Grover and Leskovec, 2016; Yang et al., 2016; Monti et al., 2016; Kipf and Welling, 2017; Hamilton et al., 2017a; b; Garcia-Duran and Niepert, 2017) , which all assume some notion of "smoothness" over the graph and thus consider nodes residing in close network proximity to be similar. Such structural role information about the nodes can be used for a variety of tasks, including as input to machine learning problems, or even to identify key nodes in a system (principal "influencers" in a social network, critical hubs in contagion graphs, etc.).

When structural roles of nodes are defined over a discrete space, they correspond to different topologies of local network neighborhoods (e.g., edge of a chain, center of a star, a bridge between two clusters). However, such discrete roles must be pre-defined, requiring domain expertise and manual inspection of the graph structure. A more powerful and robust method for identifying structural similarity involves learning a continuous vector-valued structural signature χ a of each node a in an unsupervised way. This motivates a natural definition of structural similarity in terms of closeness of topological signatures: For any > 0, nodes a and b are defined to be -structurally similar with respect to a given distance if: dist(χ a , χ b ) ≤ . Thus, a robust structural similarity metric must introduce both an appropriate signature and an adequate distance metric.

While several methods have been proposed for structural role discovery in graphs, existing approaches are extremely sensitive to small perturbations in the topology and typically lack one or more desirable properties. They often require manually hand-labeling topological features (Henderson et al., 2012) , rely on non-scalable heuristics (Ribeiro et al., 2017) , and/or return a single similarity score instead of a multidimensional structural signature (Jin et al., 2011; 2014) .

Here we address the problem of structure learning on graphs by developing GRAPHWAVE. Building upon techniques from graph signal processing (Coifman et al., 2006; Hammond et al., 2011; Shuman et al., 2013) , our approach learns a structural embedding for each node based on the diffusion of a spectral graph wavelet centered at that node. Intuitively, each node propagates a unit of energy over the graph and characterizes its neighboring topology based on the response of the network to this probe. In contrast to prior work that characterizes the wavelet diffusion as a function of the wavelet scaling parameter, we study how the wavelet diffuses through the network at a given scale as a function of the initial source node. We prove that the coefficients of this wavelet directly relate to graph topological properties. Hence, these coefficients contain all the necessary information to recover structurally similar nodes, without requiring the hand-labeling of features. However, the wavelets are, by design, localized on the graph. Therefore to compare structural signatures for nodes that are far away from each other, typical graph signal processing methods (using metrics like correlation between wavelets or 2 distance) cannot be used without specifying an exact one-to-one mapping between nodes for every pairwise comparison, a computationally intractable task.

Graph:

To overcome this challenge, we propose a novel way of treating the wavelets as probability distributions over the graph. This way the structural information is contained in how the diffusion spreads over the network rather than where it spreads. In order to provide vector-valued signatures which can then be used as input to any machine learning algorithm, we embed these wavelet distributions using the empirical characteristic function (Lukacs, 1970) . The advantage of empirical characteristic functions is that they capture all the moments of a given distribution. This allows GRAPHWAVE to be robust to small perturbations in the local edge structure, as we prove mathematically. Computational complexity of GRAPHWAVE is linear in the number of edges, thus allowing it to scale to large (sparse) networks. Finally, we compare GRAPHWAVE to several state-of-the-art baselines on both real and synthetic datasets, obtaining improvements of up to 71% and demonstrating how our approach is a useful tool for characterizing structural signatures in graphs.

Summary of contributions.:

The main contributions of our paper are as follows:

• We develop a novel use of spectral graph wavelets by treating them as probability distributions and characterizing the distributions using empirical characteristic functions. • We leverage these insights to develop a scalable method (GRAPHWAVE) for learning node embeddings based on structural similarity in graphs. • We prove that GRAPHWAVE accurately recovers structurally similar nodes.

Further related work. Prior work on discovering nodes with similar structural roles has typically relied on explicit featurization of nodes. These methods generate an exhaustive listing of each node's local topological properties (e.g., node degree, number of triangles it participates in, number of k-cliques, its PageRank score) before computing node similarities based on such heuristic representations. A notable example of such approaches is RolX (Henderson et al., 2012) , which aims to recover a soft-clustering of nodes into a predetermined number of K distinct roles using recursive feature extraction (Henderson et al., 2011) . Similarly, struc2vec (Ribeiro et al., 2017) uses a heuristic to construct a multilayered graph based on topological metrics and simulates random walks on the graph to capture structural information. In contrast, our approach does not rely on heuristics (we mathematically prove its efficacy) and does not require explicit manual feature engineering or hand-tuning of parameters.

Another line of related work are graph diffusion kernels (Coifman et al., 2006) which have been utilized for various graph modeling purposes (Kondor and Lafferty, 2002; Chung, 2007; Rustamov and Guibas, 2013; Tremblay et al., 2014) . However, to the best of our knowledge, our paper is the first to apply graph diffusion kernels for determining structural roles in graphs. Kernels have been shown to efficiently capture geometrical properties and have been successfully used for shape detection in the image processing community (Sun et al., 2009; Ovsjanikov et al., 2010; Aubry et al., 2011) . However, in contrast to shape-matching problems, GRAPHWAVE considers these kernels as probability distributions over real-world graphs. This is because the graphs that we consider are highly irregular (as opposed to the Euclidean and manifold graphs). Therefore, traditional wavelet methods, which typically analyze node diffusions across specific nodes that occur in regular and predictable patterns, do not apply. Instead, by treating wavelets as distributions, GRAPHWAVE characterizes the shape of the diffusion, rather than the specific nodes where the diffusion occurs. This key insight allows us to uncover structural signatures and to discover structurally similar nodes.

2. LEARNING STRUCTURAL SIGNATURES:

Given an undirected connected graph G = (V, E) with N nodes V = {a 1 , . . . , a N }, edges E, an adjacency matrix A (binary or weighted), and a degree matrix D ii = j A ij , we consider the problem of learning, for every node a i , a structural signature representing a i 's position in a continuous multidimensional space of structural roles.

We frame this as an unsupervised learning problem based on spectral graph wavelets (Hammond et al., 2011) and develop an approach called GRAPHWAVE that provides mathematical guarantees on the optimality of learned structural signatures.

2.1. SPECTRAL GRAPH WAVELETS:

In this section, we provide background on a spectral graph wavelet-based model (Hammond et al., 2011; Shuman et al., 2013 ) that we will use in the rest of the paper.

Let U be the eigenvector decomposition of the unnormalized graph Laplacian L = D -A and let

λ 1 < λ 2 ≤ • • • ≤ λ N (Λ = Diag(λ 1 , . . . , λ N ))

denote the eigenvalues of L. Let g s be a filter kernel with scaling parameter s. For simplicity, we use the heat kernel g s (λ) = e -λs throughout this paper, but our results apply to any low-pass filter kernel (Shuman et al., 2016) . For now, we assume that s is given; we develop a method for selecting an appropriate value of s in Appendix C.

Graph signal processing (Hammond et al., 2011; Shuman et al., 2013) defines the spectral graph wavelet associated with g s as the signal resulting from the modulation in the spectral domain of a Dirac signal centered around node a. The spectral graph wavelet Ψ a is given by an N -dimensional vector:

EQUATION

where δ a = 1(a) is the one-hot vector for node a. For notational simplicity, we drop the explicit dependency of spectral graph wavelet Ψ a on s. The m-th wavelet coefficient of this column vector is thus given by Ψ ma = N l=1 g s (λ l )U ml U al . In spectral graph wavelets, the kernel g s modulates the eigenspectrum such that the resulting signal is typically localized on the graph and in the spectral domain (Shuman et al., 2013) . Spectral graph wavelets are based on an analogy between temporal frequencies of a signal and the Laplacian's eigenvalues. Eigenvectors associated with smaller eigenvalues carry slow varying signal, encouraging nodes that are geographically close in the graph to share similar values. In contrast, eigenvectors associated with larger eigenvalues carry faster-varying signal across edges. The low-pass filter kernel g s can thus be seen as a modulation operator that discounts higher eigenvalues and enforces smoothness in the signal variation on the graph.

2.2. GRAPHWAVE ALGORITHM:

First we describe the GRAPHWAVE algorithm (Alg. 1) and then analyze it in the next section. For every node a, GRAPHWAVE returns a 2d-dimensional vector χ a representing its structural signature, where nodes with structurally similar local network neighborhoods will have similar signatures.

We first apply spectral graph wavelets to obtain a diffusion pattern for every node (Line 3), which we gather in a matrix Ψ. Here, Ψ is a N × N matrix, where a-th column vector is the spectral graph wavelet for a heat kernel centered at node a. In contrast to prior work that studies wavelet coefficients as a function of the scaling parameter s, we study them as a function of the network (i.e., how the coefficients vary across the local network neighborhood around the node a). In particular, coefficients in each wavelet are identified with the nodes and Ψ ma represents the amount of energy that node a Algorithm 1 GRAPHWAVE algorithm for learning structural signatures.

1: Input: Graph G = (V, E), scale s, evenly spaced sampling points {t 1 , t 2 , . . . , t d }. 2: Output: Structural signature χ a ∈ R 2d for every node a ∈ V 3: Compute Ψ = U g s (Λ)U T (Eq. (1)) 4: for t ∈ {t 1 , t 2 , . . . , t d } do 5: Compute φ(t) = column-wise mean(e itΨ ) ∈ R N 6:

for a ∈ V do 7:

Append Re(φ a (t)) and Im(φ a (t)) to χ a has received from node m. As we will later show nodes a and b with similar network neighborhoods have similar spectral wavelet coefficients Ψ (assuming that we know how to solve the "isomorphism" problem and find the explicit one-to-one mapping of the nodes from a's neighborhood to the nodes of the b's neighborhood). To resolve the node mapping problem GRAPHWAVE treats the wavelet coefficients as a probability distribution and characterizes the distribution via empirical characteristic functions. This is the key insight that makes it possible for GRAPHWAVE to learn nodes' structural signatures via spectral graph wavelets.

More precisely, we embed spectral graph wavelet coefficient distributions into 2d-dimensional space (Line 4-7) by calculating the characteristic function for each node's coefficients Ψ a and sample it at d evenly spaced points. The characteristic function of a probability distribution X is defined as: Lukacs, 1970) . The function φ X (t) fully characterizes the distribution of X because it captures information about all the moments of probability distribution X (Lukacs, 1970) . For a given node a and scale s, the empirical characteristic function of Ψ a is defined as:

φ X (t) = E[e itX ], t ∈ R (

EQUATION

Finally, structural signature χ a of node a is obtained by sampling the 2-dimensional parametric function (Eq. ( 2)) at d evenly spaced points t 1 , . . . t d and concatenating the values:

EQUATION

Note that we sample/evaluate the empirical characteristic function φ a (t) at d points and this creates a structural signature of size 2d. This means that the dimensionality of the structural signature is independent of the graph size. Furthermore, nodes from different graphs can be embedded into the same space and their structural roles can be compared across different graphs.

Distance between structural signatures. The final output of GRAPHWAVE is a structural signature χ a for each node a in the graph. We can explore distances between the signatures through the use of the 2 distance on χ a . The structural distance between nodes a and b is then defined as:

dist(a, b) = χ a -χ b 2

. By definition of the characteristic function, this technique amounts to comparing moments of different orders defined on wavelet coefficient distributions.

Scaling parameter. The scaling parameter s determines the radius of network neighborhood around each node a (Tremblay et al. (2014) ; Hammond et al. (2011)) . A small value of s determines node signatures based on similarity of nodes' immediate neighborhoods. In contrast, a larger value of s allows the diffusion process to spread farther in the network, resulting in signatures based on neighborhoods with greater radii.

GRAPHWAVE can also integrate information across different radii of neighborhoods by jointly considering many different values of s. This is achieved by concatenating J representations χ

(sj )

a , each associated with a scale s j , where s j ∈ [s min , s max ]. We provide a theoretically justified method for finding an appropriate range s min and s max in Appendix C. In this multiscale version of GRAPHWAVE, the final aggregated structural signature for node a is a vector χ a ∈ R 2dJ with the following form:

χ a = [Re(φ (sj ) a (t i )), Im(φ (sj ) a (t i )] ti,sj .

Computational complexity. We use Chebyshev polynomials (Shuman et al., 2011) to compute Line

3. ANALYSIS OF GRAPHWAVE:

In this section, we provide theoretical motivation for our spectral graph wavelet-based model (Shuman et al., 2013) . First we analytically show that spectral graph wavelet coefficients characterize the topological structure of local network neighborhoods (Section 3.1). Then we show that structurally equivalent/similar nodes have near-identical/similar signatures (Sections 3.2 and 3.3), thereby providing a mathematical guarantee on the optimality of GRAPHWAVE.

3.1. SPECTRAL GRAPH WAVELETS AS A MEASURE OF NETWORK STRUCTURE:

We start by establishing the relationship between the spectral graph wavelet of a given node a and the topological properties of local network neighborhood centered at a. In particular, we prove that a wavelet coefficient Ψ ma provides a measure of network connectivity between nodes a and m.

We use the fact that the spectrum of the graph Laplacian is discrete and contained in the compact set [0, λ N ]. It then follows from the Stone-Weierstrass theorem that the restriction of kernel g s to the interval [0, λ N ] can be approximated by a polynomial. This polynomial approximation, denoted as P , is tight and its error can be uniformly bounded. Formally, this means:

EQUATION

where K is the order of polynomial approximation, α k are coefficients of the polynomial, and r(λ) = g s (λ) -P (λ) is the residual. We can now express the spectral graph wavelet for node a in terms of the polynomial approximation as:

EQUATION

We note that Ψ a is a function of L k = (D -A) k and thus can be interpreted using graph theory.

In particular, it contains terms of the form D k (capturing the degree), A k (capturing the number of k-length paths that node a participates in), and terms containing both A and D, which denote paths of length up to k going from node a to every other node m.

Using the Cauchy-Schwartz's inequality and the facts that U is unitary and r(λ) is uniformly bounded (Eq. ( 4)), we can bound the second term on the right-hand side of Eq. ( 5) by:

EQUATION

)

As a consequence, each wavelet Ψ a can be approximated by a K-th order polynomial that captures information about the K-hop neighborhood of node a. The analysis of Eq. ( 5), where we show that the second term is limited by , indicates that spectral graph wavelets are predominately governed by topological features (specifically, degrees, cycles and paths) according to the specified heat kernel.

The wavelets thus contain the information necessary to generate structural signatures of nodes.

3.2. SIGNATURES OF STRUCTURALLY EQUIVALENT NODES:

Let us consider nodes a and b whose K-hop neighborhoods are identical (where K is an integer less than the diameter of the graph), meaning that nodes a and b are structurally equivalent. We now show that a and b have -structurally similar signatures in GRAPHWAVE.

First, we use the Taylor expansion to obtain an explicit K-th order polynomial approximation of g s as:

P (λ, s) = K k=0 (-1) k (sλ) k /k!.

Then, for each eigenvalue λ, we use the Taylor-Lagrange equality to ensure the existence of c λ ∈ [0, s] such that:

EQUATION

)

If we take any s such that it satisfies: s ≤ ((K + 1)! )) 1/(K+1) /λ 2 , then the absolute residual |r(λ)| in Eq. ( 7) can be bounded by for each eigenvalue λ. Here, is a parameter that we can specify depending on how close we want the signatures of structurally equivalent nodes to be (note that smaller values of the scale s lead to smaller values of and thus tighter bounds).

Because a and b are structurally equivalent, there exists a one-to-one mapping π from the Khop neighborhood of a (i.e., N K (a)) to the K-hop neighborhood of b (i.e., N K (b)), such that:

N K (b) = π(N K (a)

). We extend the mapping π to the whole graph G by randomly mapping the remaining nodes. Following Eq. ( 5), we write the difference between each pair of mapped coefficients Ψ ma and Ψ π(m)b in terms of the K-th order approximation of the graph Laplacian:

|Ψ ma -Ψ π(m)b | = δ m U (P (Λ) + r(Λ) U T δ a -δ π(m) U (P (Λ) + r(Λ) U T δ b ≤ (U P (Λ)U T ) ma -(U P (Λ)U T ) π(m)a + (U r(Λ)U T ) ma + (U r(Λ)U T ) π(m)b . (8)

Here, we analyze the first term on the second line in Eq. ( 8). Since the K-hop neighborhoods around a and b are identical and by the localization properties of the k-th power of the Laplacian (k-length paths, Section 3.1), the following holds:

∀m ∈ N K (a), ( K k=0 α k L k ) ma = ( K k=0 α k L k ) π(m)b , ∀m / ∈ N K (a), ( K k=0 α k L k ) ma = ( K k=0 α k L k ) π(m)a = 0,

meaning that this term cancels out in Eq. ( 8). To analyze the second and third terms on the second line of Eq. ( 8), we use bound for the residual term in the spectral graph wavelet (Eq. ( 6)) to uniformly bound entries in matrix U r(Λ)U T by . Therefore, each wavelet coefficient in Ψ a is within 2 of its corresponding wavelet coefficient in Ψ b , i.e., |Ψ ma -Ψ π(m)b | ≤ 2 . As a result, because similarity in distributions translates to similarity in the resulting characteristic functions (Lévy's continuity theorem), then assuming the appropriate selection of scale, structurally equivalent nodes have -structurally similar signatures.

3.3. SIGNATURES OF STRUCTURALLY SIMILAR NODES:

We now analyze structurally similar nodes, or nodes whose K-hop neighborhoods are identical up to a small perturbation of the edges. We show that such nodes have similar GRAPHWAVE signatures.

Let N K (a) denote a perturbed K-hop neighborhood of node a obtained by rewiring edges in the original K-hop neighborhood N K (a). We denote by L the graph Laplacian associated with that perturbation. We next show that when perturbation of a node neighborhood is small, the changes in the wavelet coefficients for that node are small as well.

Formally, assuming a small perturbation of the graph structure (i.e., sup ||L k -Lk || F ≤ , for all k ≤ K), we use K-th order Taylor expansion of kernel g s to express the wavelet coefficients in the perturbed graph as:

EQUATION

We then use the Weyl's theorem (Coburn et al., 1966) to relate perturbations in the graph structure to the change in the eigenvalues of the graph Laplacian. In particular, a small perturbation of the graph yields small perturbations of the eigenvalues. That is, for each λ, r( λ) is close its original value r(λ):

r( λ) = r(λ) + o( ) ≤ C

, where C is a constant. Taking everything together, we get:

|Ψ ma -Ψ ma | ≤ | K k=0 α k (L k -Lk ) ma | + | Ũ r( Λ) Ũ T | ma + |U r(Λ)U T | ma = ( K k=0 |α k | + 1 + C) ,

indicating that structurally similar nodes have similar signatures in GRAPHWAVE.

4. EXPERIMENTS ON SYNTHETIC GRAPHS:

Baselines. We compare our GRAPHWAVE method against two state-of-the-art baselines, struc2vec (Ribeiro et al., 2017) and RolX (Henderson et al., 2012) . We note that RolX requires the number of desired structural classes as input, whereas the two other methods learn embeddings that capture a continuous spectrum of roles rather than discrete classes. We thus use RolX as an oracle estimator, providing it with the correct number of classes1 . We also note that homophily-based methods (Kipf and Welling (2017) ; Hamilton et al. (2017a) , etc.) are unable to recover structural similarities. 

4.1. BARBELL GRAPH:

We consider a barbell graph consisting of two dense cliques connected by a long chain (Figure 2A ). We run GRAPHWAVE, RolX, and struc2vec and plot a 2D PCA representation of learned structural signatures in Figure 2B-D .

GRAPHWAVE correctly learns identical representations for structurally equivalent nodes, providing empirical evidence for our theoretical result in Section 3.2. This can be seen by structurally equivalent nodes in Figure 2A (nodes of the same color) having identical projections in the PCA plot (Figure 2D ). In contrast, both RolX and struc2vec fail to recover the exact structural equivalences.

All three methods correctly group the clique nodes (purple) together. However, only GRAPHWAVE correctly differentiates between nodes connecting the two dense cliques in the barbell graph, providing empirical evidence for our theoretical result in Section 3.3. GRAPHWAVE represents those nodes in a gradient-like pattern that captures the spectrum of structural roles of those nodes (Figure 2D ).

4.2. COMPLEX AND NOISY GRAPH STRUCTURES:

Graphs. We next consider four types of synthetic graphs where the structural role of each node is known and used as ground truth information to evaluate performance. The graphs are given by basic shapes of one of different types ("house", "fan", "star") that are regularly placed along a cycle (Table 1 and Figure 3A ). In the "varied" setup, we mix the three basic shapes when placing them along a cycle, thus generating synthetic graphs with richer and more complex structural role patterns. Additional graphs are generated by placing these shapes irregularly along the cycle followed by adding a number of edges uniformly at random. In our experiments, we set this number to be around 5% of the edges in the original structure. This setup is designed to assess the robustness of the methods to data perturbations ("house perturbed", "varied perturbed").

Experimental setup. For each graph, we run RolX, struc2vec, and GRAPHWAVE to learn the signatures. We choose to use a multiscale version of GRAPHWAVE where the scale was set as explained in Appendix C. We then use k-means to cluster the learned signatures and use three standard metrics to evaluate the clustering quality. (1) Cluster homogeneity is the conditional entropy of the ground-truth structural roles given the proposed clustering (Rosenberg and Hirschberg, 2007) .

(2) Cluster completeness (Rosenberg and Hirschberg, 2007) evaluates whether nodes with the same structural role are in the same cluster.

(3) Silhouette score compares the mean intra-cluster distance to the mean between-cluster distance, assessing the density of the recovered clusters. This score takes a value in [-1,1] (higher is better).

Results. GRAPHWAVE consistently outperforms struc2vec, yielding improvements for the homogeneity of up to 50%, and completeness up to 69% in the "varied" setting (Table 1 ). Both GRAPHWAVE and RolX achieved perfect performance in the noise-free "house" setting, however, GRAPHWAVE outperformed RolX by up to 4% (completeness) in the more complex "varied" setting. We evaluated methods on graphs in the presence of noise ("perturbed" in Table 1 ): GRAPHWAVE outperformed RolX and struc2vec by 4% and 50% (homogeneity), respectively, providing empirical evidence for our analytical result that GRAPHWAVE is robust to noise in the edge structure. The silhouette scores also show that the clusters recovered by GRAPHWAVE are denser and better separated than for the other methods.

As an example, we show a cycle graph with attached "house" shapes (Figure 3A ). We plot 2D PCA projections of GRAPHWAVE's signatures in Figure 3B , confirming that GRAPHWAVE accurately distinguishes between nodes with distinct structural roles. We also visualize the resulting characteristic functions (Eq. ( 2)) in Figure 3C . In general, their interpretation is as follows (Appendix D):

• Nodes located in the periphery of the graph struggle to diffuse the signal over the graph, and thus span wavelets that are characterized by a smaller number of non-zero coefficients. Characteristic functions of such nodes thus span a small loop-like 2D curve. • Nodes located in the core (dense region) of the graph tend to diffuse the signal farther away and reach farther nodes for the same value of t. Characteristic functions of such nodes thus have a farther projection on the x and y axis.

In Figure 3C , different shapes of the characteristic functions capture different structural roles. We note the visual proximity between the roles of the blue, light green and red nodes that these curves carry, as well as their clear difference with the core dark green and purple nodes.

5. EXPERIMENTS ON REAL-WORLD GRAPHS:



5.1. THE ENRON EMAIL GRAPH:

Data and setup. Nodes represent Enron employees and edges correspond to email communication between the employees (Klimt and Yang, 2004) . An employee has one of seven functions in the company (e.g., CEO, president, manager) . These functions provide ground-truth information about roles of the corresponding nodes in the network. We use GRAPHWAVE to learn a structural signature for every Enron employee. We then use these signatures to compute the average 2 2 distance between every two categories of employees. Results. GRAPHWAVE captures intricate organizational structure of Enron (Figure 4 ). For example, CEOs and presidents are structurally distant from all other job titles. This indicates their unique position in the email exchange graph, which can be explained by their local graph connectivity patterns standing out from the others. Traders, on the other hand, appear very far from presidents and are closer to directors. In contrast, struc2vec is less successful at revealing intricate relationships between the job titles, yielding an almost uniform distribution of distances between every class. We assess the separation between "top" job titles (CEO and President) and lower levels in the job title hierarchy. GRAPHWAVE achieves 28% better homogeneity and 139% better completeness than RolX. We also note that the variability within each cluster of struc2vec is higher than the average distance between clusters (dark green colors on the diagonal, and lighter colors on the off-diagonal).

5.2. EUROPEAN AIRLINE GRAPHS:

Data and setup. The airline graphs are taken from the list of airlines operating flights between European airports (Cardillo et al., 2013) . Each airline is represented with a graph, where nodes represent airports and links stand for direct flights between two airports. Given an airline graph, we use GRAPHWAVE to learn structural signature of every airport. We then create a visualization using PCA on the learned signatures to layout the graph on a two-dimensional structural space.

Results. Figure 5 shows graph visualizations of two very different airlines. Air France is a national French airline, whose graph has the so-called hub and spoke structure, because the airline is designed to provide an almost complete coverage of the airports in France. We note that the signature of CDG (Charles De Gaulle, which is Air France's central airport) clearly stands out in this 2D projection, indicating its unique role in the network. In contrast to Air France, Ryanair is a low-cost airline whose graph avoids the overly centralized structure. Ryanair's network has near-continuous spectrum of structural roles that range from regional French airports (Lille (LIL), Brest Bretagne (BES)) all the way to London Stansted (STN) and Dublin (DUB). These airlines have thus developed according to different structural and commercial constraints, which is clearly reflected in their visualizations.

6. CONCLUSION:

We have developed a new method for learning structural signatures in graphs. Our approach, GRAPHWAVE, uses spectral graph wavelets to generate a structural embedding for each node, which we accomplish by treating the wavelets as a distributions and evaluating the resulting characteristic functions. Considering the wavelets as distributions instead of vectors is a key insight needed to capture structural similarity in graphs.

Our method provides mathematical guarantees on the optimality of learned structural signatures. Using spectral graph theory, we prove that structurally equivalent/similar nodes have nearidentical/similar structural signatures in GRAPHWAVE. Experiments on real and synthetic networks provide empirical evidence for our analytical results and yield large gains in performance over stateof-the-art baselines. For future work, these signatures could be used for transfer learning, leveraging data from a well-explored region of the graph to infer knowledge about less-explored regions.

NOTATION:

In the appendix, we use the same notation as in the main paper. However, in the main paper, we dropped the explicit dependency of the heat kernel wavelet Ψ on the scaling parameter s. We note that we explicitly keep this dependency throughout the appendix as we study the relationship between heat kernel wavelets and the scaling parameter.

A KNOWN PROPERTIES OF HEAT KERNEL WAVELETS:

We here list five known properties of heat kernel wavelets that we later use to derive a method for automatic selection of the scaling parameter in Appendix C.

P1. The heat kernel wavelet matrix Ψ (s) is symmetric:

Ψ (s) ma = δ T m U e -sΛ U T δ a = δ T a U e -sΛ U T δ m = Ψ (s) am .

P2. First eigenvector: By definition of the graph Laplacian L = D -A, the vector 1 is an eigenvector of L, corresponding to the smallest eigenvalue λ 1 = 0. This means that:

EQUATION

*) where the last equality follows from the fact that eigenvectors in U are orthogonal.

P3. Scaling inequality: Using the Cauchy Schwartz inequality, we get:

(Ψ (s) ma ) 2 = ( N j=1 e -λj s U aj U mj ) 2 ≤ ( N j=1 (e -λ j s 2 U aj ) 2 )( N j=1 (e -λ j s 2 U jm ) 2 ) = Ψ (s) aa Ψ (s) mm .

It thus follows that: 

EQUATION

|Ψ (s+1) aa - 1 N | = | N j=2 e -λj (s+1) U 2 ja | ≤ | N j=2 e -λj s U 2 ja | ≤ |Ψ (s) aa - 1 N |,

where:

|Ψ (s) aa - 1 N | = e -λ2s | N j=2 e -(λj -λ2)s U 2 ja |.

decreasing function of s:

Taking everything together, we get:

e -λ2s U 2 a2 ≤ |Ψ (s) aa -1 N | ≤ e -λ2s (1 -1 N ). It thus follows that: |Ψ (s+1) aa -1 N |/|Ψ (s) aa -1 N | ≤ 1.

B PROPERTIES OF HEAT KERNEL WAVELET DISTRIBUTIONS:

We prove three propositions about distributions generated by heat kernel wavelets that will be used in Appendix C.

Proposition 1. The mean of heat kernel wavelet Ψ (s)

a is equal to µ (s) a = 1 N N m=1 Ψ (s) ma = 1 N . The mean µ (s)

a is thus independent of the value of scaling parameter s and node a.

Proof. We expand the mean of heat kernel wavelet Ψ (s) a using spectral graph wavelet definition in Eq. ( 1) and Property 2 (P2) from Appendix A:

µ (s) a = 1 N N m=1 Ψ (s) ma = 1 N 1 T U e -Λs U δ a = 1 √ N δ T 1 e -Λs U δ a = 1 √ N δ T 1 (e -λ1s U 1a , e -λ2s U 2a , . . . , e -λ N s U N a ) T = 1 √ N e -λ1s U 1a = 1 N .

Proposition 2. The heat kernel wavelet coefficient Ψ (s)

aa at the initial source node a is a monotonically decreasing function of scaling parameter s. Its value is bounded by: 1

N ≤ Ψ (s) aa ≤ 1.

Proof. This follows directly using definition Ψ (s) aa = m j=1 e -λj s U 2 ja and the scaling inequality in Property 3 (P3). This way, we get:

Ψ (∞) aa = 1 N ≤ Ψ (s) aa ≤ 1 = Ψ (0) aa .

Additionally, for any m = a, the wavelet coefficient

Ψ (s)

ma is non-negative and bounded. Specifically, the wavelet coefficient can be written as: Ψ Proof. We use the definition of variance to get:

Var[Ψ (s) a ] = 1 N N m=1 (Ψ (s) ma -µ (s) a ) 2 = 1 N N m=1 (Ψ (s) ma ) 2 - 1 N 2 .

We rewrite the sum on the far right hand-side using the symmetry of wavelet matrix Ψ (Property 1):

N m=1 (Ψ (s) ma ) 2 = N m=1 Ψ (s) am Ψ (s) ma = (Ψ 2 ) aa = δ T a U e -2sΛ U T δ a = ||e -sΛ U T δ a || 2 2 = N j=1 e -2sλj U 2 aj , concluding that variance Var[Ψ (s)

a ] is decreasing, since it is sum of functions, all of which decrease as s gets larger.

C SCALING PARAMETER OF HEAT KERNEL:

We here develop a method that automatically finds an appropriate range of values for the scaling parameter s in heat kernel g s , which we use in the multiscale version of GRAPHWAVE (Section 2.2)

We find the appropriate range of values for s by specifying an interval bounded by s min and s max through the analysis of variance in heat kernel wavelets. Intuitively, whether or not a given value for s is appropriate for structural signature learning depends on the relationship between the scaling parameter and the temporal aspects of heat equation. In particular, small values of s allow little time for the heat to propagate, yielding diffusion distributions (i.e., heat kernel wavelet distributions) that are trivial in the sense that only a few coefficients have non-zero values and are thus unfit for comparison. For larger values of s the network converges to a state in which all nodes have an identical temperature equal to 1/N (Property 4 (P4)), meaning that diffusion distributions are data-independent, hence non-informative.

Next we prove Propositions 4 and 5 to provide new insights into the variance and convergence rate of heat kernel wavelets. We then use these results to select s min and s max .

Proposition 4. Given the scaling parameter s, the variance of off-diagonal coefficients in heat kernel wavelet Ψ a is proportional to:

Var[{Ψ (s) am

; m = a}] ∝ ∆ (0) a ∆ (2s) a -(∆ (s) a ) 2 , where ∆ (s) a = |Ψ (s) aa -1

N | is a monotonically decreasing function of s.

Proof. Let us denote the mean of off-diagonal coefficients in wavelet Ψ a by: μ(s)

a = m =a Ψ (s)

ma /(N -1). We use the fact that m =a Ψ (s)

ma = 1 -Ψ (s)

aa , along with the definition of the variance, to obtain:

Var[{Ψ (s) am ; m

= a}] = 1 N -1 m =a (Ψ (s) ma -μ(s) a ) 2 = 1 N -1 m =a (Ψ (s) ma ) 2 -(μ (s) a ) 2 = 1 N -1 (Ψ (2s) aa -(Ψ (s) aa ) 2 ) - 1 (N -1) 2 (1 -Ψ (s) aa ) 2 = N (N -1) 2 (Ψ (2s) aa N -1 N -(Ψ (s) aa ) 2 + 2Ψ (s) aa N - 1 N ) = N (N -1) 2 ((Ψ (2s) aa - 1 N ) N -1 N -(Ψ (s) aa - 1 N ) 2 ) = N (N -1) 2 (∆ (0) a ∆ (2s) a -(∆ (s) a ) 2 ).

Proposition 4 proves that the variance is a function of ∆

a . Therefore, to maximize the variance, we must analyze the behavior of ∆ (s) a . To ensure sufficient variability in the distribution of wavelet coefficients, we need to select a range [s min , s max ] that bounds the ∆ (s) a . Our goal thus becomes establishing that ∆ (s) a is large enough that the diffusion has had time to spread, while remaining sufficiently small to ensure that the diffusion is far from its converged state.

Proposition 5. The convergence of heat kernel wavelet coefficient Ψ (s) am is bounded by: e -λ N s ∆ (0) a ≤ ∆ (s) a ≤ e -λ2 s ∆ (0) a .

Proof. We use Property 5 from Appendix A and induction over s to complete this proof. For a given s ≥ 0 we analyze:

|Ψ (s+1) aa - 1 N | = | N j=2 e -λj (s+1) U 2 ja | ≤ e -λ2 |Ψ (s) aa - 1 N |,

and conclude that: a is smooth increasing function of s, we can take the floor/ceiling of any non-integer s ≥ 0 and this proposition must hold.

e -λ N ≤ |Ψ (s+1) aa -1 N |/|Ψ (s) aa -1 N | ≤ e -λ2 .

C.1 SELECTION OF s max:

We select s max such that wavelet coefficients are localized in the network. To do so, we use Proposition 5 and bound ∆ a is closer to e -λ N (i.e., lower bound in Proposition 5). When the bulk of the eigenvalues is closer to λ 2 , ∆ (s) a will lean towards e -λ2 (i.e., upper bound in Proposition 5). In each case, the diffusion is localized if ∆ (s) a is above a given threshold η < 1. Indeed, this ensures that ∆ (s) a has shrunk to at most η * 100 % of its initial value at s = 0, and yields a bound of the form: ∆ (s) a /∆ (0) a ≥ η. The bound implies that: e -λs ≥ η, or s ≤ -log(η)/λ. To find a middle ground between the two convergence scenarios, we take λ to be the geometric mean of λ 2 and λ N . Indeed, as opposed to the arithmetic mean, the geometric mean maintains an equal weighting across the range [λ 2 , λ N ], and a change of % in λ 2 has the same effect as a change in % of λ N . We thus select s max as: s max = -log(η)/ √ λ 2 λ N .

C.2 SELECTION OF s min:

We select s min to ensure the adequate diffusion resolution. In particular, we select a minimum value s min such that each wavelet has sufficient time to spread. That is, ∆

a /∆ (0) a ≤ γ. As in the case of s max above, we obtain a bound of s ≥ -log(γ)/λ. Hence, we set s min to: s min = -log(γ)/ √ λ 2 λ N .

To cover an appropriate range of scales, we suggest setting η = 0.90 and γ = 0.99.

D VISUALIZATION OF CHARACTERISTIC FUNCTIONS:

Here, we study the properties of characteristic functions in GRAPHWAVE. Our goal is to provide intuition to understand the behavior of these functions and how their resulting 2D parametric curves reflect nodes' local topology (see Figure 3C ).

We begin by reviewing the definition of the characteristic function.

Definition 1. The empirical characteristic function (Lukacs, 1970) of wavelet Ψ (s)

a is function φ (s) a defined as:

φ (s) a (t) = 1 N N m=1

e itΨ (s) ma , t ∈ R.

In the phase plot, a given value of t thus yields the following set of coordinates: .

By varying t, we get a characteric curve in this 2D plane. Here, we note several properties of this curve:

1. Value at t = 0: ∀s, φ 3. Behavior as s → ∞: From Property 4 in Appendix A, the coefficients all converge to the same limit: ∀m, a, lim s→∞ Ψ (s) ma = 1 N . This in turn implies that: ∀t, lim s→∞ φ (s) a (t) = e i t N . Hence, the curve converges to a circle centered at 0, with radius 1 and period 2πN . .

Since all the wavelet coefficients are non-negative, the curve is thus directed counterclockwise.

:

in Algorithm 1. As inDefferrard et al. (2016), each power of the Laplacian has a computational cost of O(|E|), yielding an overall complexity of O(K|E|), where K denotes the order Chebyshev polynomial approximation. The overall complexity of GRAPHWAVE is linear in the number of edges, which allows GRAPHWAVE to scale to large sparse networks.

All code and data for experiments can be downloaded at: https://goo.gl/ndBkwb.

Figures and Tables:

Figure 1: Nodes a and b have similar local structural roles even though they are distant in the graph. While raw spectral graph wavelet signatures/coefficients Ψ of a and b might be very different, we treat them as probability distributions and show that the coefficient distributions are indeed similar.

Figure 2: Barbell graph. The graph has 8 distinct classes of structurally equivalent nodes as indicated by color (A). 2D PCA projection of structural signatures as learned by RolX (B), struc2vec (C) and GRAPHWAVE (D). Projections in B-D contain the same number of points as there are nodes in the barbell graph. Identical signatures have identical projections, resulting in overlapping points in B-D.

Figure 4: Heat maps indicate average distance between roles in the Enron email graph, as determined by struc2vec (A) and GRAPHWAVE (B). Table of comparative statistics for the three algorithms.B

mm ). P4. Convergence: We use the fact that when s → ∞, the value of e -sλ converges to 0 for any non-zero eigenvalue λ of the graph Laplacian. Combining this fact with the spectral graph wavelet definition in Eq. (1), we get: U aj U mj = U a1 U m1 = 1/N, since it follows from Eq. (**) that U a1 = U 1a = 1/ √ N for every a. The diffusion over the network thus converges to a state where all nodes have an identical temperature 1/N . P5. Convergence rate: To analyze the convergence rate of heat kernel wavelets as a function of s, we expand |Ψ (s+1) aa -1 N | using the spectral graph wavelet definition in Eq. (1):

λj s U ja U jm = 1 N + N j=2 e -λj s U ja U jm .It can be bounded by: 0 ≤ Ψ The variance of heat kernel wavelet Var[Ψ (s) a ] is a strictly decreasing function of scaling parameter s.

Given any s ∈ N, we use the induction principle to get:e -λ N s |Ψ (0)aa -

by the graph Laplacian's eigenvalues. When the bulk of the eigenvalues leans towards λ N , ∆ (s)

the scale or initial source node.2.Behavior for s = 0: Only the coefficient corresponding to the initial source node has non-zero value (i.e., Ψ(0)aa = 1, and ∀m = a, Ψ(0)ma = 0). Hence, ∀t, φ

1 N e it + N -1 N . The phase plot of the characteristic curve is a circle of radius 1 N with period 2π centered at N -1 N .

Gradient of the curve: ∇ t φ

Structural role discovery results for different synthetic graphs. (Best seen in color.) Results averaged over 20 synthetically generated graphs. Dashed lines denote perturbed versions of the basic shapes (obtained by randomly adding and removing edges), node colors indicate structural roles.

Table content: <table><tr><td>A</td><td>B</td><td>GraphWave</td><td>C</td><td/></tr><tr><td/><td/><td/><td>Imaginary</td><td/></tr><tr><td/><td/><td/><td/><td>Real</td></tr><tr><td colspan="6">Figure 3: A cycle graph with attached "house" shapes (A). 2D PCA projection of GRAPHWAVE's</td></tr><tr><td colspan="6">structural signatures. The representation of structurally equivalent nodes overlap, and GRAPHWAVE</td></tr><tr><td colspan="6">perfectly recovers the 6 different node types (B). Empirical characteristic function for the distribution</td></tr><tr><td colspan="6">of the wavelet coefficients (C). Color of a node/curve indicates structural role. (Best seen in color.)</td></tr><tr><td colspan="2">Shapes placed along a cycle graph</td><td>Method</td><td colspan="3">Homogeneity Completeness Silhouette</td></tr><tr><td/><td/><td>RolX</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td>House</td><td/><td>struc2vec</td><td>0.845</td><td>0.885</td><td>0.584</td></tr><tr><td/><td/><td>GRAPHWAVE</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td/><td/><td>RolX</td><td>0.769</td><td>0.783</td><td>0.496</td></tr><tr><td>House</td><td/><td>struc2vec</td><td>0.537</td><td>0.554</td><td>0.511</td></tr><tr><td>perturbed</td><td/><td>GRAPHWAVE</td><td>0.802</td><td>0.807</td><td>0.566</td></tr><tr><td/><td/><td>RolX</td><td>0.927</td><td>0.915</td><td>0.746</td></tr><tr><td>Varied</td><td/><td>struc2vec</td><td>0.629</td><td>0.561</td><td>0.588</td></tr><tr><td/><td/><td>GRAPHWAVE</td><td>0.943</td><td>0.948</td><td>0.867</td></tr><tr><td/><td/><td>RolX</td><td>0.776</td><td>0.728</td><td>0.525</td></tr><tr><td>Varied</td><td/><td>struc2vec</td><td>0.537</td><td>0.482</td><td>0.540</td></tr><tr><td>perturbed</td><td/><td>GRAPHWAVE</td><td>0.807</td><td>0.806</td><td>0.589</td></tr></table>

Peer review (iclr2018):
 The paper proposes a method for quantifying the similarity between the local neighborhoods of nodes in a graph/network. There are many ways in which such a distance/similarity metric between nodes could be defined. For example, once could look at the induced subgraph G_i formed by the k-neighborhood of node i, and the induced subgraph G_j of the k-neighborhood of node j, and define the similarity as k(G_i,G_j) where k is any established graph kernel. Moreover, the task is unsupervised, which makes it hard to compare the performance of different methods. Most of the experiments in the paper seem a bit contrived. Regarding the algorithm, the question is: “sure, but why this way?”. The authors take the heat kernel matrix on the graph, treat each column as a probability distribution, compute its characteristic function, and define a distance between characteristics functions. This seems pretty arbitrary and heuristic. I also find it confusing that they refer to the heat kernel as wavelets. The spectral graph wavelets of Hammond et al is a beautiful construction, but, as far as I remember, it is explicitly emphasized that the wavelet generating function g must be continuous and satisfy g(0)=0. By setting g(.lambda)=e^{-s .lambda}, the authors just recover the diffusion/heat kernel of the graph. That’s not a wavelet. Why call this a “spectral graph wavelet” approach then? The heat kernel is much simpler. I find this misleading. I also feel that the mathematical results in the paper have little depth. Diffusion is an inherently local process. It is natural then that the diffusion matrix can be approximated by a polynomial in the Laplacian (in fact, it is sufficient to look at the power series of the matrix exponential). It is not surprising that the diffusion function captures some local properties of the graph (there are papers by Reid Andersen/ Fan Chung/ Kevin Lang, as well as by Mahoney, I believe on localized PCA in graphs following similar ideas). Again, there are many ways that this could be done. The particular way it is done in the paper is heuristic and not supported by either math or strong experiments.