Spectral Embedding of Graph Networks

We derive the proposed graph representation from the curve enclosing the feasibility domain of vertex and spectral spreads of the BCG. The definition below is given for a general weighted graphs.

Note that $\mathbb{D}_{(s_{\mathbf{L}},s_{\mathbf{W}})}$ contains all possible values $(s_{\mathbf{L}},s_{\mathbf{W}})$ that can be obtained on $\cal G$ with respect to a non zero signal $\mathbf{x}\in l^{2}(\cal G)$.

The lower boundary of the curve enclosing the feasibility domain of the spreads $\left(\mathbf{s}_{\mathbf{L}^{BE}}(\mathbf{x}),\mathbf{s}_{\mathbf{W}^{BE}}(\mathbf{x})\right)$, with respect to a unit norm vector $\mathbf{x}\in l^{2}(\cal G)$ is defined as

To solve (7), we use the Lagrangian

with $\beta\in\mathbb{R}$. Differentiating the Lagrangian and comparing the results to zero, we obtain the following generalized eigenvalue problem:

where the eigenvector $\mathbf{x}$ solving (9) is also a minimizer for (7).

Denoting

we can write the generalized eigenvalue problem (9) using the matrix pencil $\mathbf{\tilde{L}}(\beta)$

The scalar $\beta\in\mathbb{R}$ can be considered as a parameter that controls the trade-off between the total vertex and spectral spreads. An example for the trade-off in the barbell graph as captured by $\tilde{\mathbf{L}}^{BE}(\beta$) is illustrated in Fig. (a) and (b) using $\beta$ = -200, and $\beta$ = -0.2, respectively. Functions colored with green correspond to the eigenvectors of $\tilde{\mathbf{L}}^{BE}(\beta$), functions colored with blue correspond to the eigenvectors of ${\mathbf{L}}^{BE}$ . As can be seen, when $|\beta|$ is large, the eigenvectors of $\tilde{\mathbf{L}}^{BE}(\beta$) reveal structure which is similar to the eigenvectors of ${\mathbf{L}}^{BE}$, while small values of $|\beta|$ produce structure which is similar to those corresponding to ${\mathbf{W}}^{BE}$.

For analysis and practical considerations it is often useful to encode the network captured in $\mathbf{\tilde{L}}$ using a semi positive definite operator. For practical consideration, we transform $\mathbf{\tilde{L}}$ into a semi positive definite matrix $\mathbf{L}_{\Delta}$ using a simple perturbation matrix $\Delta$, where $\Delta=\mu\mathbf{I}$, $\mu=-\tilde{\lambda}_{0}$, where $\tilde{\lambda}_{0}$ is the smallest eigenvalue of $\mathbf{\tilde{L}}$, and $\mathbf{I}$ is the identity matrix. Setting

ensures that $\mathbf{L}_{\Delta}$ is semi positive definite (SPD) matrix. Letting $\left\{\tilde{\lambda}_{i}\right\}_{i=1}^{N}$ and $\left\{\lambda^{\Delta}_{i}\right\}_{i=1}^{N}$ be the eigenvalues corresponding to $\mathbf{\tilde{L}}$ and $\mathbf{{L}}_{\Delta}$, respectively, we can see that the choice made in (12) ensures that the original spacing in eigenvalues of $\mathbf{\tilde{L}}$ is preserved in $\mathbf{{L}}_{\Delta}$ $\lambda^{\Delta}_{i+1}-\lambda^{\Delta}_{i}=\tilde{\lambda}_{i+1}-\tilde{\lambda}_{i}$. The embedding method proposed using $\mathbf{L}_{\Delta}$ is coined Graph Spectral Spread Embedding (GSSE). Note that there is a geometric interpretation which is related to the way $\mathbf{\tilde{L}}(\beta)$ was obtained from the lower boundary curve enclosing the feasibility domain $\mathbb{D}_{(s_{\mathbf{L}},s_{\mathbf{W}})}$; Setting instead another perturbation matrix $\Delta=\mu\mathbf{I}$ with $\mu=-\lambda_{N}^{W}$ and defining

where ${\mbox{min}}({\lambda}(\mathbf{L}_{\Delta}(\beta)))$ corresponds to the minimum eigenvalue of $\mathbf{L}_{\Delta}$ . we have that

which defines a half plane in $\mathbb{D}_{(s_{\mathbf{L}},s_{\mathbf{W}})}$.

Properties: The feasibility domain $\mathbb{D}_{(s_{\mathbf{L}^{BE}},s_{\mathbf{W}^{BE}})}$ is a bounded set since $\forall(s_{\mathbf{L}^{BE}},s_{\mathbf{W}^{BE}})\in\mathbb{D}_{(s_{\mathbf{L}^{BE}},s_{\mathbf{W}^{BE}})}$

where $\lambda_{N}^{W^{BE}}$corresponds to the largest eigenvalue of the affinity graph. We provide quantitative analysis that bound the largest eigenvalue of the betweeness graph centrality, $\lambda_{N}^{W^{BE}}$ by the largest value of the vertex betweeness centrality of the graph (see Lemma 2 in the appendix).

In this section we construct a node embedding using a linear mapping between the subspaces $\mathbf{W}^{BE}$ and $\mathbf{{L}}^{BE}$, which capture local and structural properties. We exploit the fact that the correspondence between the nodes represented by the two subspaces is known, and propose using the solution of a Sylvester equation as the node feature representation. This choice is motivated by Roth’s removal rule, which states that one can find a complementary subspace by solving the Sylvester equation. $\mathbf{W}^{BE}$and $\mathbf{{L}}^{BE}$ can be considered two different operators with known correspondence, the identity map. Thus, the resulting embedding given by the eigensystem of the Sylvester equation will be composed of two hybrid representations associated with the graph networks and local connectivity. The method is coined Graph Sylvester Embedding (GSE).

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

for $\mathbf{C}=\rho\mathbf{I}$, a scalar regularizer, $\mathbf{I}$ is the identity matrix, where $\mathbf{L}^{BE},\mathbf{W}^{BE},\mathbf{C}\in\mathbb{R}^{N\times N}$. The resulting eigensystem of $\mathbf{X}$ is described in the next Lemma below.

Applying Lemma 1 with $\mathbf{A}=\mathbf{L}^{BE}$, and $\mathbf{B}=\mathbf{W}^{BE}$ we obtain that the set of eigenvectors corresponding to the solution $\mathbf{X}$ to the Sylvester equation are $\left\{\phi_{l}^{BE}+\phi_{l}^{W^{BE}}\right\}_{l=1}^{N}$, where $\phi_{l}^{BE},\phi_{l}^{W^{BE}}$ correspond to the set of eigenvectors of $\mathbf{L}^{BE}$, $\mathbf{W}^{BE}$ , respectively.
The proposed embedding is constructed from the set of eigenvectors and associated eigenvalues corresponding to $\mathbf{L}_{\Delta}^{BE}$ and $\mathbf{X}$, respectively. The performance of these embeddings is now evaluated empirically.

Datasets. We compare our method using the Barbell graph (Figure 1) for a qualitative comparison, and the mirrored Karate network for a quantitative comparison of node classification of structural similarity. The mirrored Karate network (based on a experiment suggested in [30]) is composed of two copies corresponding to the Karate Club network [38], where each node has a mirror node in the duplicated graph. We follow the exact same protocol which was suggested in [19], where the performance is evaluated by measuring the percentage of nodes whose nearest neighbor in the embedding space corresponds to its structural equivalent role.

Barbell Graph Our graph embeddings is able to identify nodes 5 and 16 as ones that stand out and have a large distance in the embeddings space from the rest of the nodes to which they are connected to (nodes 1 to 4 and 17 to 20, that are considered locally similar to 5 and 16, respectively). For example, the heat maps shown in Fig.1(b) and (c) represent embeddings using the same number of eigenvectors corresponding to ${\mathbf{L}}$ and ${\mathbf{L}^{BE}}$ respectively. In these cases, nodes 5 and 16 do not stands out and their distances to locally similar nodes is very small. Karate network We compare to a number of state-of-the-art methods in node feature embeddings including node2vec [22] struc2vec [30], Deep Walk [29], and GraphWave [19]. The baseline which is the WKS descriptor [4] combined with the Laplacian of the unweighted graph is also reported in Table 1, where our approach shows a significant improvement over the baseline. The experimental results in Table (showing the best and average results for the varying number of edges) demonstrates that our method consistently performs well for a varying number of mirrored edges, and is improving the accuracy by more than $14\%$ comparing to the best method.

Effective representation of nodes with similar network topology plays an important role in the problem of protein protein interaction, where graph networks that describe protein-protein interaction (PPI) encode valuable information that can provide insights of biological function. To understand the development of COVID-19, the analysis of PPI networks is ongoing [26]. It is important to develop robust computational methods that can reveal biological processes and host signaling pathways that may be used by SARS-CoV-2 to infect cells.

One promising approach aims to develop drugs that target human proteins that the virus requires. Recent research responding to COVID-19 aims to prioritize additional human protein by employing the known human protein interactors of SARS-CoV-2 proteins and a whole-genome protein interaction network in order to predict potential interactors with high accuracy [25].

Another important problem in studying PPI networks is the problem of network alignment or node correspondence which, by finding topological similarities between networks, can be useful to understanding biological function and evolution. We consider a problem of node classification by finding node correspondence using network alignment on a recently published dataset of human proteins that physically interact with SARS-CoV-2 [21].
Dataset: We test using the STRING network from the STRING database, a PPI network that consists of 18,886 nodes and 977,789 edges. The STRING network includes 330 human proteins that physically interact with SARS-CoV-2. It is likely to have both false positives and false negatives edges. For instance, the relative expression of the viral protein could increase the number of interacting partners detected, generating false positives [25].

We aim to find network alignment between two copies of the STRING network, given the STRING network and an additional STRING network which is created by randomly removing 10$\%$ and 20 $\%$ percentage of noisy edges. Each of the networks includes 330 nodes corresponding to the known human protein interactors of SARS-CoV-2. Given the node correspondence of the protein interactors of SARS-CoV-2, we connect each pair of nodes corresponding to the same human protein interactors of SARS-CoV-2 with an edge, which is resulted in a network composed of the two STRING networks. We then apply the proposed GSSE and GSE to generate node embeddings. For comparison, we compute the node embeddings using graph embeddings methods including Laplacian Eigenmaps (LE) [5], Locally Linear Embedding (LLE)[31], Hope [28], and Regularized Laplacian (RL)[41], which are given the same input network. Additionally, we compare to network alignment methods, such as Isorank [33] and Final [39] and iNEAT [40]. For the network alignments methods, we provide as an input the corresponding affinity graphs and a matrix with the known correspondence between the nodes that correspond to human proteins that physically interact with SARS-CoV-2.

Table 4 shows the percentage of nodes whose nearest neighbor corresponds to its true node correspondence using GSE and GSSE in the STRING network, which is compared to Graph Embeddings and network alignment methods. 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. Network alignments methods perform better, while our proposed GSE is improving the accuracy by more than $\sim 16.4\%$ comparing to the best method.

Forecasting fracture locations in a progressively failing disordered structure is a crucial task in structural materials [7]. Recently, networks were used to represent 2 dimensional (2D) disordered lattices and have been shown a promising ability to detect failures locations inside 2D disordered lattices. A promising aspect of this line of work is the possibility to assess failure locations of a structure without needing to study its detailed energetic states [7]. In [7], it has been shown that the EBC can serve as a meaningful measurement to forecast the failure locations observed on a contact networks in 2D granular network. As failure locations were shown to occur predominately at locations that have large EBC values, [7] proposed a thresholding scheme to compare the EBC of each edge to the mean EBC in order to forecast edges corresponding to failed beams.
Dataset: The set of disordered structures was derived from experimentally determined force networks in granular materials [6]. The network data is available in the Dryad repository [8]. We tested 6 different initial networks, with mean degrees $z={2.40,2.55,2.60,3.35}$, 3.0, 3.6, following the same datasets corresponding to different initial granular configurations. For GSSE and GSE we used $\beta=-1000$ and $\rho=1$, respectively, in all experiments. Table 3 shows the success rates in forecasting failed edges for the 6 networks. In comparison to the forecasting approach proposed in [8], our framework performs consistently better. In particular, the proposed GSE which is based on Sylvester embedding provides the best results for all networks. We employ the tools developed in Sect. 4, with simple adaptation to create edge embedding. First, by constructing node embedding, and then using the concatenated nodes features to construct edge embeddings. While our method is guided by the finding that failed edges occur predominantly in edges above the mean EBC, it doesn’t rely on a single thersholding criteria to asses failed edges. Instead, we defer to make decisions about failing edges in the edge embedding space, which provides an effective tool to process, analyze and visualize the network. Fig. 3 (b) and (c) are showing visualization using t-SNE of our proposed methods for edge embedding. The points colored in red correspond to the edges whose values are above the mean EBC, while the points colored in turquoise correspond to edges bellow the mean EBC. The blue enlarged dots correspond to the failed edges in the system which were successfully detected by each method. In the examples illustrated, our proposed embeddeding successfully forecast all failed edges, which were mapped to the same cluster including edges whose EBC value was below the mean value. Table 3 shows a quantitative results: the forecasting approach proposed [7] will miss failed edges whose EBC is below the EBC mean, while embeddings using standard spectral embedding was found even less successful.

Computational Complexity: 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 a 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.
The local clustering coefficient: Another limitation is concerned with applying our method to graphs in which the average local clustering coefficient is small, in particular when the graph contains a large number of vertices whose clustering coefficient is zero. If the degree of a node is 0 or 1, the corresponding clustering coefficient of that node is 0. In clustering tasks, our approach is shown to improves the baseline when the average local clustering coefficient is not too small, and is less effective when the average local clustering coefficient is very small.

Future work include expanding our approach to address semi-supervised applications, for example by incorporating generalized graph Laplacian or BCG in training using Graph convolutional networks. Another direction is to explore alternative approaches for finding a trade-off between local and structural similarity, tailored to a specific domain application using a functional applied to the adjacency graph in a similar spirit to the approach we propose with the BCG.

We report additional experimental results on the SARS- CoV-2 PPI networks for the task of node classification using network alignment. We follow a similar protocol as in the previous experiments that aim to find network alignment between two copies of the STRING network. An additional STRING network is created by randomly removing 20$\%$ of noisy edges from the STRING network. In this experiment, we randomly subsample the set of the 332 nodes corresponding to the known SARS-CoV-2 interactors and use only 163 nodes as the given known node correspondence. Since Graph Embeddings are found less effective for this task, we only compare to network alignments methods, including Isorank (Singh, Xu, and Berger 2008), Final (Zhang and Tong 2016), and iNEAT (Zhang et al. 2020). Table 1 shows the percentage of nodes whose nearest neighbor corresponds to its true node correspondence using GSE in the STRING network. The comparison results demonstrate that our method remains robust even when using a much smaller number of correspondence nodes.

Datasets. The European air-traffic network [30] is a data collected from the statistical office of the European union. The network has 399 nodes and 5,995 edges. Airport activity is measured by the total number of landings plus takeoffs in the corresponding period [30]. The Cora dataset [Cora] is a citation network where nodes are documents and edges are citation links. The network has 2,708 nodes and 5,429 edges. Evaluation: We evaluate our approach in clustering accuracy comparing to the ground truth using purity measure. We compute node feature representation for each node and then construct a new graph which is based on the $k=15$ nearest neighbor graph based on the new features. We compare our method to a number of clustering methods, including Spectral Clustering [1] and Affinity Propagation.
Summary of the experimental results: on the European air-traffic network The experimental results, compare to spectral clustering and affinity propagation and measured in clustering accuracy by purity evaluation, results in an improvement of 24.1$\%$ and 28.2 $\%$ using GSSE, and 20$\%$ and 20.3 $\%$ using GSE, respectively, while using BCG directly we obtain an improvement of 6 $\%$ and 10$\%$.

Citation network For the Cora dataset, the performance comparing to spectral clustering and affinity propagation using the GSSE and GSE shows an improvements of less than $1\%$, while using the BCG shows an improvement of 29.35$\%$ comparing to spectral clustering and affinity propagation.