Graph Spectral Embedding using the Geodesic Betweeness Centrality

Analysis of networks arising in social science, material science, biology and commerce, relies on both local and global properties of the graph that represents the network. Local properties include similarity of attributes among adjacent nodes, and global properties include the structural role of nodes or edges within the graph, as well as similarity among distant nodes. Graph representations based on the eigenvectors of the standard Graph Laplacian have limited expressivity when it comes to capturing local and structural similarity, which we wish to overcome.

We propose an unsupervised approach to learn embeddings that exploit both the Laplacian and the affinity matrix combined into a Sylvester equation. However, to better capture local as well as structural similarity, we consider the Laplacian and the affinity matrix not of the original graph, but rather of a modified betweenness centrality graph derived from it. Centrality is a statistic (a function of) the topology of the graph that is not captured in the original edge weights. The node embeddings are the Spectral Kernel Descriptors of the solution of the Sylvester equation.

The key to our approach is three-fold: (i) The betweeness centrality of edges, defined as the fraction of shortest paths going through that edge, to capture semi-global structure in the original graph; (ii) The use of both the Laplacian and the affinity matrix of the centrality graph, to modulate the effect of local and global structure; (iii) The use of the Sylvester equation as a natural way of combining these two graph properties into a single basis. Any number of descriptors could be built on this basis, we just choose the Spectral Kernel Descriptor for simplicity.

Our approach is related to uncertainty principles of graph harmonic analysis, that explore the extent in which signals can be represented simultaneously in the vertex and spectral domains. We extend this approach to characterize the relation between the spectral and non-local graph domain spread of signals defined on the nodes.

We illustrate the use of our method in multiple unsupervised learning tasks including node classification in social networks, network alignment in protein-to-protein interaction (PPI) for the human-SARS-CoV-2 protein-protein interactome, and forecasting failure of edges in material science.

For a graph ${\cal G}=(\cal V,W)$ with vertices in the set $\cal V$ $=\left\{1,2,...,N\right\}$ and edge weights in $\mathbf{\cal{W}}=\left\{(w_{ij})(i,j)\in\cal{V}\right\}$, where $w_{ij}$ denotes the weight of an edge between nodes $i$ and $j$, we measure the centrality of an edge by the fraction of the number of shortest paths that pass through that edge, called edge betweenness centrality (EBC). From the betweenness centrality of the edges of ${\cal G}$, we construct a modified graph with the same connectivity of ${\cal G}$, but edges weighted by their EBC value. The result is called edge betweeness centrality graph (BCG), or centrality graph for short, and denoted by $\bar{\cal G}$, which measures the importance of edges given the adjacency matrix.

Graph centrality based measurements, such as vertex centrality, eigenvector centrality [20], and edge betweeness centrality [28] have been widely used in network analysis and its diverse applications [7, 10]. There is a large body of work that focus on graph embedding methods that, similar to our approach, aims to learn a function $\Phi:\cal V\rightarrow$ $\mathbb{R}^{m}$ from the graph nodes into a vector space $\mathbb{R}^{m}$ [23, 1, 30, 40]. Such methods construct graph representations that are based on either local or structural similarities.

Graph convolutional networks [24, 11, 39] are among the most popular methods for graph based feature learning, combining the expressiveness of neural networks with the graph structures. Typically, these methods assume that node features or attributes are available, while we focus on problems where only the graph network is available. There is also a large body of work that focus on graph embedding methods that, similar to our approach, only assumes that the affinity matrix is given as an input. Such methods construct graph representations that are based on either local or structural similarities [23, 1, 30, 40].

Uncertainty principles in graphs have been explored by [2, 37], extending traditional uncertainty principles of signal processing to more general domains such as irregular and unweighted graphs [32]. [2] provides definitions of graph and spectral spreads, where the graph spread is defined for an arbitrary vertex, and studies to what extent a signal can be localized both in the graph and spectral domains.

Other related work includes different approaches in manifold learning [31, 5, 38], manifold regularization [21, 17, 13, 16, 12, 18, 14, 15], graph diffusion kernels [9] and kernel methods widely used in computer graphics [3, 36] for shape detection. Most methods assume that signals defined over the graph nodes exhibit some level of smoothness with respect to the connectivity of the nodes in the graph and therefore are biased to capture local similarity.

The Sylvester equation was previously employed in a variety of graph mining applications such as network alignment and subgraph matching [41, 19]. Most previous works used the continuous Sylvester equation with two input affinity matrices corresponding to the input affinity matrices. We propose a different approach, where we solve the discrete-time Sylvester equation using the affinity matrix and its associated Laplacian.

Consider an undirected, weighted graph $\cal{G}=(V,W)$. The degree $d(i)$ of a node is the sum of the weights of edges connected to it. The combinatorial graph Laplacian is denoted by $\mathbf{L}$, and defined as $\mathbf{L}=\mathbf{D}-\mathbf{W}$, with $\mathbf{D}$ the diagonal degree matrix with entries $d_{ii}=d(i)$. The eigenvalues and eigenvectors of $\mathbf{L}$ are denoted by $\lambda_{1},\ldots,\lambda_{N}$ and $\mathbf{\phi}_{1},\ldots,\phi_{N}$, respectively. The normalized Laplacian is defined as ${\mathbf{L}}_{N}=\mathbf{D}^{-\frac{1}{2}}\mathbf{L}\mathbf{D}^{-\frac{1}{2}}=\mathbf{I}-\mathbf{D}^{-\frac{1}{2}}\mathbf{W}\mathbf{D}^{-\frac{1}{2}}$ and its real eigenvalues are in the interval $[0,2]$. In this work we use the normalized Laplacian $\mathbf{L}_{N}$, which from now on we refer to as $\mathbf{L}$ for simplicity.

We suggest an embedding that goes beyond the scalar modulation of the generalized Laplacian Eq.(17) by using a linear mapping between the subspaces $\mathbf{W}^{BE}$ and $\mathbf{{L}}^{BE}$. We exploit the fact that the correspondence between the nodes represented by the two basis functions of $\mathbf{W}^{BE}$and $\mathbf{{L}}^{BE}$ is known, i.e., it is the identity map, $\pi:{\mathbf{W}}^{BE}\rightarrow{\mathbf{L}}^{BE}$ with $\pi({\mathbf{W}}^{BE}(:,i))={\mathbf{L}}^{BE}(:,i)$ for each $i\in\cal V$ (i.e., the column index corresponding to each node $i\in\cal V$ in ${\mathbf{W}}^{BE}$ is the same in ${\mathbf{L}}^{BE}$). We propose using the solution of a Sylvester equation as the node feature representation, which is also the mapping between the nodes’ representations. The resulting embedding is given by the solution to the discrete Sylvester equation (or the singular value decomposition associated with the solution) will be composed of two hybrid representations associated with the graph networks and local connectivity. The proposed method is coined Graph Sylvester Embedding (GSE).

The discrete-time Sylvester operator $\mathbf{S}(\mathbf{X})=\mathbf{A}\mathbf{X}\mathbf{B}-\mathbf{X}$ is used to express the eigenvalues and eigenvectors of $\mathbf{A}$ and $\mathbf{B}$ using a single operator $\mathbf{S}$, where we solve

using $\mathbf{C}=\mathbf{I}$, ($\mathbf{I}\in\mathbb{R}^{N\times N}$ is the identity matrix), $\mathbf{A}=\mathbf{W}^{BE},\mathbf{B}=\mathbf{L}^{BE}$.

Let $\left\{\phi_{i}^{\mathbf{W}^{BE}}\right\}_{l=1}^{N}$, $\left\{\lambda_{l}^{\mathbf{W}^{BE}}\right\}_{l=1}^{N}$, and $\left\{\phi_{i}^{\mathbf{L}^{BE}}\right\}_{l=1}^{N}$, $\left\{\lambda_{l}^{\mathbf{L}^{BE}}\right\}_{l=1}^{N}$ be the corresponding eigenvectors and eigenvalues of $\mathbf{A}$ and $\mathbf{B}$, respectively.

We will assume that the eigenvalues of $\mathbf{A}$ and $\mathbf{B}$ satisfy that $\lambda_{i}^{\mathbf{L}^{BE}}\lambda_{j}^{\mathbf{W}^{BE}}\neq-1$, $\forall i,j$, hence the operator $\mathbf{S}(\mathbf{X}):\mathbb{R}^{N\times N}\rightarrow\mathbb{R}^{N\times N}$ is non-singular and has $N^{2}$ matrix eigenvalues and eigenvectors. Note that the Sylvester equation has a unique solution for any $\mathbf{C}$ if and only if $\mathbf{S}$ is non-singular, which occurs if and only if $\lambda_{i}^{\mathbf{L}^{BE}}\lambda_{j}^{\mathbf{W}^{BE}}\neq-1$ $\forall i,j$. By choosing $\mathbf{C}=\mathbf{I}$ and ordering the columns in $\mathbf{A}$ and $\mathbf{B}$ based on the identity map, $\pi:\mathbf{A}\rightarrow\mathbf{B}$ we obtain that each column in $\mathbf{X}$ provides a node representation to its corresponding node $i\in\cal V$.

Note that we can also express the solution $\mathbf{X}$ to the Sylvester equation Eq.(9) using the Kronecker product $\otimes$:

where we vectorize $\mathbf{X}$ and $\mathbf{C}$ to obtain its equivalent vector representation $x$ = vec $(\mathbf{X})$ and $c$ = vec$(\mathbf{C})$. Since the Sylvester equation can be written using the Kronecker product, then all known properties of the Kronecker product are carried on to the proposed graph Sylvester embedding. One important property is invariance to permutations.

Proof: See Appendix.

The resulting solution of $\mathbf{X}$ can be described using the following conditions. Assume that $\lambda_{i}^{\mathbf{W}^{BE}}\lambda_{j}^{\mathbf{L}^{BE}}\neq 1\,\forall i,j$, and that the graph representation encoded in $\mathbf{A}$ and $\mathbf{B}$ has the same order based on the identity map, and $\mathbf{C}=\mathbf{I}$. Then the solution $\mathbf{X}$ to the discrete-time Sylvester system Eq.(9) is unique where

where $\mathbf{\tilde{C}}_{i,j}=\frac{(\phi_{i}^{\mathbf{W}^{BE}})^{t}\phi_{j}^{\mathbf{L}^{BE}}}{\lambda_{i}^{\mathbf{W}^{BE}}\lambda_{j}^{\mathbf{L}^{BE}}-1}$. See proof in the appendix.

Remark 1: Eq. (12) shows the analytical solution to the Sylvester equation. In practice, one can use fast iterative methods [26]. Also note that under the assumptions that $\lambda_{i}^{\mathbf{W}^{BE}}\lambda_{j}^{\mathbf{L}^{BE}}\neq 1\,\forall i,j$, the solution $\mathbf{X}$ is unique [26]. We note that Sylvester’s equation can be generalized to include multiple terms, thus allowing one to incorporate multiple basis functions. We illustrate the concept with two bases, to capture local and global statistics of the topology of the network, and leave extension to additional bases, that can be specific to particular domains or tasks, for future work.

Remark 2: One can view the solution $\mathbf{X}$ to the Sylvester equation as a polynomial of the matrices $\mathbf{A}$ and $\mathbf{B}$ ([22]), where using higher-order polynomials in $\mathbf{W}^{BE}$, $\mathbf{L}^{BE}$ can be interpreted as imposing smoothness in the embedding space.

We evaluate our method on real-world networks in several applications including material science and network alignment of Protein-Protein interactions (PPIs) networks using a recent dataset that has been used to study Covid-19, the SARS-CoV-2. We compare GSE to known and existing state of the art methods in the respective applications. There are a number of benchmarks to evaluate network embeddings built using at least partial supervision, for instance for node classification relative to a given ground truth. Our method does not require any, and therefore we focus our attention on unsupervised benchmarks. Extension of our method to graph classification, or other graph analysis tasks that can exploit supervision is beyond the scope of this paper and will be considered in future work.
Evaluation metrics: In the problem of network alignment we use the percentage of correct node correspondence found to evaluate our method, based on the known ground truth correspondence between the two networks. In detecting failure edges in material science we measure the success rate using sensitivity. To detect failured edges, we use the graph embedding obtained by each method to cluster the data into two approximately equal size clusters using spectral clustering . We then choose the cluster with the largest betweenness centrality mean value. For further validation, we test the statistical significance of each method by computing the $p$ values using hyper-geometric distribution (see the Appendix).
Comparison with other methods: Our baseline for comparison with the proposed Graph Sylvester Embedding is using the concatenated spectral embedding corresponding to both the Laplacian and the affinity matrix. We coin the concatenated bases as St.($\mathbf{W}^{BE},\mathbf{L}^{BE}$). We compare to the spectral Wave Kernel Descriptor signature (WKS) [3] (coined $\mathbf{L}$ desc.). $(-)$ in the Tables indicates methods that failed to converge or provided meaningless results.
We compare with a representative of graph embeddings methods including Laplacian Eigenmaps [5], node2vec [1],  and NetMF [23]. We also test against methods that were specifically tailored to this applications we explored including methods for graph alignments [41, 42] and methods based on betweenness centrality for detecting failured edges [7].

We have presented an embedding to represent local and global structure of a graph network, constructed without the need for any supervision nor explicit annotation. Besides the cost of time and effort, annotating data can create privacy and security risks. Our focus is on developing expressive and flexible representations that can be used in a variety of downstream tasks without human intervention. Flexibility is exercised through the choice of bases, that are combined through Sylvester’s equation. In the specific cases we have experimented, the bases are chosen to capture local and structural similarity across the graph network.

The implementation of our approach is computationally intensive, with most of the burden falling on the computation of the BCG, which is $O(mN)$ where $m$ is the total number of edges and $N$ is the total number of nodes, thus approximately $O(N^{2})$ for sparse graphs. The execution time of our method with Python code implementation using Intel Core i7 7700 Quad-Core 3.6GHz with 64B Memory on the Cora dataset with approximately 2700 nodes takes $\approx$ 24.9 seconds.

Computing BC measurements on large graphs is computationally heavy ($O(N^{2})$ for sparse graphs), and developing fast methods for BC measurements is an ongoing research area. With significant progress that has been made in recent years (e.g: [27]) we can foresee extensions to scale our approach to networks with millions of nodes. In addition, the effectiveness of BC measurements degrades on larger graphs due to noise (for example when the graph includes a large number of nodes with clustering coefficient close to zero). Applying our approach to small sub-graphs or “Network motifs” before aggregating it to to the entire graph may be one foreseeable solution.

In our study, we found that spectral embedding by solving the Sylvester equation of the edge betweeness centrality graph reveals the node’s network structure, which is not necessarily possible with the original graph weights. It is also possible to compute edge betweenness centrality in a way that takes into account the original weight information. For example, by looking into the path evaluation function that assesses a path between two nodes that is combining both the sum of edge weights and the number of shortest distance path as was proposed in [34]. Another possible direction is to solve a generalized Sylvester equation which would incorporate both matrices corresponding to the original weight information and the edge betweenness centrality. Other future work includes expanding our approach to address applications which include dynamic graphs and multi-layer graphs, which would include solving Sylvester equation with time-varying coefficients.

We report additional results and details in the problem of network alignment using the STRING network from the STRING database. In the experimental results we used all 332 human proteins known to interact with the Sars-Cov-2 as the available node correspondence between the two networks we tested on network alignment. From the 18,886 node of the String network we extract the 1000 nodes which corresponds to the nodes with the largest diffusion scores and the 332 nodes corresponding the human proteins known to interact with the Sars-Cov-2, based on the method suggested in [25]. Alignment is performed between two copies of the STRING network, where 10$\%$ and 20$\%$ edges were removed from the sub-network consisting of a total of 1332 nodes. It is evident that graph embedding methods such as LLE and LE, which are rooted in manifold learning that is biased to local similarity and heavily relies on the graph smoothness are not effective for this task.

The discrete-time Sylvester operator $\mathbf{S}(\mathbf{X})=\mathbf{A}\mathbf{X}\mathbf{B}-\mathbf{X}$ is used to express the eigenvalues and eigenvectors of $\mathbf{A}$ and $\mathbf{B}$ using a single operator $\mathbf{S}$, where we solve

using $\mathbf{C}=\mathbf{I}$, ($\mathbf{I}\in\mathbb{R}^{N\times N}$ is the identity matrix), $\mathbf{A}=\mathbf{W}^{BE},\mathbf{B}=\mathbf{L}^{BE}$.

Proof: For a permutation matrix $\mathbf{P}\in\mathbb{R}^{N\times N}$, let $\mathbf{\tilde{A}}=\mathbf{P}^{T}\mathbf{A}\mathbf{P}$, and $\mathbf{\tilde{B}}=\mathbf{P}^{T}\mathbf{B}\mathbf{P}$. Let $\mathbf{X},\mathbf{\tilde{X}}$ are two solutions to the discrete-time Sylvester equation Eq.(24) with the associate matrices $\mathbf{{A}},\mathbf{{B}}$ and $\mathbf{\tilde{A}},\mathbf{\tilde{B}}$, respectively, $\mathbf{A}\mathbf{X}\mathbf{B}-\mathbf{X}=\mathbf{I}$ and $\mathbf{\tilde{A}}\mathbf{\tilde{X}}\mathbf{\tilde{B}}-\mathbf{\tilde{X}}=\mathbf{I}$. Next note that $\mathbf{P}^{T}\mathbf{A}\mathbf{P}\mathbf{\tilde{X}}\mathbf{P}^{T}\mathbf{B}\mathbf{P}-\mathbf{\tilde{X}}=\mathbf{I}$. Since $\mathbf{P}$ is a permutation matrix, we have that $\mathbf{P}^{-1}=\mathbf{P}^{T}$ and by multiplying $\mathbf{P}$ and $\mathbf{P}^{T}$ from the left hand and right hand side, respectively, we obtain $\mathbf{A}\mathbf{P}\mathbf{\tilde{X}}\mathbf{P}^{T}\mathbf{B}-\mathbf{\tilde{X}}=\mathbf{I}$. Denoting $\mathbf{Y}=\mathbf{P}\mathbf{\tilde{X}}\mathbf{P}^{T}$, under the assumptions that $\mathbf{{X}}$ is a unique solution, then $\mathbf{{X}}=\mathbf{Y}=\mathbf{P}\mathbf{\tilde{X}}\mathbf{P}^{T}$.$\square$

The resulting solution of $\mathbf{X}$ is described in the next Lemma below. Note that since $\mathbf{A}=\mathbf{W}^{BE},\mathbf{B}=\mathbf{L}^{BE}$ we have $\mathbf{A}$ and $\mathbf{B}$ that are diagonalizable, (with the matrices $\Phi^{\mathbf{W}^{BE}},\Phi^{\mathbf{L}^{BE}}$ corresponding to the associated eigenvectors of $\mathbf{W}^{BE},\mathbf{L}^{BE}$, respectively).

Proof: Using $\mathbf{A}=\mathbf{W}^{BE},\mathbf{B}=\mathbf{L}^{BE}$ we have that $\mathbf{A}$ and $\mathbf{B}$ are diagonalizable, $\mathbf{A}=\Phi^{\mathbf{W}^{BE}}\Lambda(\Phi^{\mathbf{W}^{BE}})^{-1}$ and $\mathbf{B}=\Phi^{\mathbf{L}^{BE}}\Lambda(\Phi^{\mathbf{W}^{BE}})^{-1}$ where $\Lambda^{\mathbf{W}^{BE}}=\mbox{diag}(\lambda_{1}^{\mathbf{W}^{BE}},...,\lambda_{N}^{\mathbf{W}^{BE}})$ and $\Lambda^{\mathbf{L}^{BE}}=\mbox{diag}(\lambda_{1}^{\mathbf{L}^{BE}},...,\lambda_{N}^{\mathbf{L}^{BE}})$. Since the matrices $\Phi^{\mathbf{W}^{BE}},\Phi^{\mathbf{L}^{BE}}$ are orthogonal, we have

Multiplying by $(\Phi^{\mathbf{W}^{BE}})^{T}$ and $\Phi^{\mathbf{L}^{BE}}$ we obtain

Setting $\mathbf{\tilde{C}}=(\Phi^{\mathbf{W}^{BE}})^{T}\mathbf{X}\Phi^{\mathbf{L}^{BE}}$ and using $\mathbf{C}=\mathbf{I}$, we have

Since $\Lambda^{\mathbf{W}^{BE}}$ is a diagonal matrix, we have that the $i$th row of $\Lambda^{\mathbf{W}^{BE}}\mathbf{\tilde{C}}\Lambda^{\mathbf{L}^{BE}}$ is $\lambda_{i}$ times the $i$th row of $\mathbf{\tilde{C}}\Lambda^{\mathbf{L}^{BE}}$and since $\Lambda^{\mathbf{L}^{BE}}$ is a diagonal matrix, then the $i$th column of $\Lambda^{\mathbf{W}^{BE}}\mathbf{\tilde{C}}\Lambda^{\mathbf{L}^{BE}}$ is $\lambda_{i}^{\mathbf{L}^{BE}}$ times the $i$th column of $\Lambda^{\mathbf{W}^{BE}}\mathbf{\tilde{C}}$. Combining the two properties, we obtain that for the $(i,j)$ entry

hence: $\mathbf{\tilde{C}}_{i,j}=\frac{(\phi_{i}^{\mathbf{W}^{BE}})^{T}\phi_{i}^{\mathbf{L}^{BE}}}{\lambda_{i}^{\mathbf{W}^{BE}}\lambda_{j}^{\mathbf{L}^{BE}}-1}$, with $\mathbf{X}=\Phi^{\mathbf{W}^{BE}}\mathbf{\tilde{C}}(\Phi^{\mathbf{L}^{BE}})^{T}$.$\square$
The next Proposition is a result of Lemma 2.1 in [22], showing the solution $\mathbf{X}$ to the Sylvester equation as a polynomial of the matrices $\mathbf{A}$ and  $\mathbf{B}$ [22], where in our case using higher-order polynomials in $\mathbf{W}^{BE}$, $\mathbf{L}^{BE}$ can be interpreted as imposing smoothness in the embedding space.

Proof: Let $\mathbf{X}$ be a solution to Eq (24), we assume that $\mathbf{A}\in\mathbb{R}^{N\times N}$and $\mathbf{B}\in\mathbb{R}^{N\times N}$are symmetric matrices with the same dimensionality, and that $\sigma(\mathbf{A})\cap\sigma(\mathbf{B})=\emptyset$. By Lemma 2.1 in [22] $\mathbf{X}$ is unique and

for any $k\geq 1$ . Taking $\mathbf{C}=\mathbf{I}$, we have that ${\mathbf{A}}^{k}\mathbf{X}-\mathbf{X}{\mathbf{B}}^{k}=\sum_{i=0}^{k-1}\mathbf{A}^{k-1-i}\mathbf{C}{\mathbf{B}}^{i}=\sum_{i=0}^{k-1}\mathbf{A}^{k-1-i}{\mathbf{B}}^{i}$. Since $\sigma(\mathbf{A})\cap\sigma(\mathbf{B})=\emptyset$ then $\sigma(\mathbf{A}^{d})\cap\sigma(\mathbf{B}^{d})=\emptyset$. Substituting $\mathbf{A}$ and $\mathbf{B}$ with $\mathbf{A}^{d}$and $\mathbf{B}^{d}$, respectively, then by similar arguments to Proposition 2.3. in [22] $\mathbf{X}$ can be represented as a polynomial of $\mathbf{A}^{d}$ and $\mathbf{B}^{d}$ using the Cayley–Hamilton Theorem with $\mathbf{X}=q(\mathbf{A}^{d})^{-1}\eta(k,\mathbf{A}^{d},\mathbf{B}^{d})$, where $q(\mathbf{A}^{d})$ is the characteristic polynomial of $\mathbf{A}^{d}$ . $\square$