Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.
10.1109/ACCESS.2017.DOI

This work is supported by National Natural Science Foundation of China (No.
61673086) and the Science Promotion Programme of UESTC, China (No.
Y03111023901014006).

Corresponding author: Shi-Min Cai (e-mail: shimin.cai81@gmail.com).

# SSNE: Effective Node Representation for Link Prediction in Sparse Networks

MING-REN CHEN1  PING HUANG1  YU LIN4  SHI-MIN CAI1,2,3  School of Computer
Science and Engineering, University of Electronic Science and Technology of
China, Chengdu 611731, China  Institute of Fundamental and Frontier Sciences,
University of Electronic Science and Technology of China, Chengdu 611731,
China  Big data research center, University of Electronic Science and
Technology of China, Chengdu 611731, China  Research School of Computer
Science, Australian National University, Canberra, Australian

###### Abstract

Graph embedding is gaining its popularity for link prediction in complex
networks and achieving excellent performance. However, limited work has been
done in sparse networks that represent most of real networks. In this paper,
we propose a model, Sparse Structural Network Embedding (SSNE), to obtain node
representation for link predication in sparse networks. The SSNE first
transforms the adjacency matrix into the Sum of Normalized $H$-order Adjacency
Matrix (SNHAM), and then maps the SNHAM matrix into a $d$-dimensional feature
matrix for node representation via a neural network model. The mapping
operation is proved to be an equivalent variation of singular value
decomposition. Finally, we calculate nodal similarities for link prediction
based on such feature matrix. By extensive testing experiments bases on
synthetic and real sparse network, we show that the proposed method presents
better link prediction performance in comparison of those of structural
similarity indexes, matrix optimization and other graph embedding models.

###### Index Terms:

link prediction, graph embedding, node representation, sparse network

=-15pt

##  I Introduction

IN real complex systems, there are many entities, which interact with each
other in a complicated way. By treating these entities as nodes and the
corresponding interactive relationships as edges, we can abstract such systems
into the network (or graph) model. Naturally, diverse types of complex
networks are available to represent real complex systems, such like social
networks, traffic networks, brain and biological networks, infrastructure
networks, etc [1, 2, 3]. Complex networks are constantly evolving over time
and new connections between entities may occur in the future. Therefore, link
prediction becomes an important task to study the dynamic evolution of network
structure[4, 5, 6, 7, 8, 9].

In previous researches, a relatively simple link prediction framework is
proposed based on the assumption that the greater the similarity between two
nodes in the network, the greater the possibility of connection between them
[5]. Then, many similarity measurements of nodes have been proposed to compute
similarity-based indexes for link prediction. A network contains a huge amount
of structural information, which has been modeled as many similarity-based
indexes including the common neighbor (CN) index [10], the Adamic-Adar (AA)
index [11], the resource allocation (RA) index [12], the Katz index [13], the
restarted random walk (RWR) index [14], and the SimRank index [15], etc. These
indexes can mainly be divided into two categories, local and global structural
similarity indexes. The local structural similarity indexes (e.g., CN, AA, RA)
only use the local topological information of nodes, which benefit from low
computational complexity and become suitable for large-scale networks.
However, their accuracy is slightly lower compare to that of the global
structural similarity indexes (e.g., Katz, RWR, and SimRank) which considers
the global topological information at a higher computational cost.

Link predication has been studied by Markov chain and machine learning
methods. For example, Sarukkai used Markov chain for link prediction and path
analysis [16], and Zhu _et al._ extended Markov chain method to link
prediction for adaptive networks [17]. With the development of representation
learning, graph embedding (i.e., graph representation) has been widely used in
link prediction [18, 19]. Graph embedding is able to map graph into low-
dimension vector space, and at the same time, keep the structure feature and
inherent attribute of graph [20, 21, 22]. Commonly, its pivotal is to sample
enough structural information by random walks on graph (or network). DeepWalk
[23] is one of the most popular random-walk-based graph embedding models. The
link prediction method based on DeepWalk is shown to better predict the
possible incidence of MicroRNA genetic disease [24, 25], as well as individual
multiple interests or attributes [26, 27].

Although these random-walk-based graph embedding models succeed in link
prediction in many real networks, they involve critical experience-tuned
parameters, such like the sampling length of a random walk and the number of
random walks [23]. In a typical scenario, they may only specify a locally
maximum value within a finite interval of experience-tuned parameters, and the
error accumulation of multiple parameters would also lead to hindered link
prediction performance in sparse networks. Therefore, in the framework of
graph embedding, we propose a novel graph embedding model, Sparse Structure
Network Embedding (SSNE), to solve the above-mentioned problems of random-
walk-based graph embedding models. The SSNE includes two important operations,
the first is to transform the adjacency matrix of sparse network into a
general mutual information matrix based on the algorithm of the Sum of
Normalized H-order Adjacency Matrix (SNHAM) and the second is to map the SNHAM
matrix into a $d$-dimensional feature matrix for effective node representation
via a neural network model. Its details will be introduced in the Section IV.

To verify the effectiveness of SSNE for link prediction, we design experiments
over a variety of datasets. The experimental results based on sparse networks
show that the link prediction method based on SSNE outperforms other methods
based on structural similarity indexes, matrix optimization and other graph
embedding models, respectively. As for relatively dense and better structured
networks, SSNE still shows comparable performance as structural similarity
indexes, which is better than matrix optimization. In short, in this paper, we
make the following contributions:

  * •

We propose a novel graph embedding model that overcomes the drawbacks in the
popular random-walk-based graph embedding models. The SNHAM algorithm is used
to transform the adjacency matrix into a new matrix of theoretical co-
occurrence probability between each pair of nodes, which substitutes the
sampling method in random-walk-based graph embedding models. Meanwhile, we
testify the mapping operation of SNHAM algorithm to be an equivalent variation
of the singular value decomposition (SVD), which greatly improves the
computational efficiency of obtaining feature matrix.

  * •

We construct a link prediction method based SSNE. The results in the testing
experiments based on six real networks and two types of artificial network
models show the good performance of link prediction in sparse networks.

  * •

We verify the algorithmic stability of link prediction method based SSNE by
selecting different combinations of parameters. The results show that the
proposed method is generally insensitive to parameters.

The remainder of this paper is organized as follows. In the Section II, we
briefly survey related works. The Section III gives the problem definition.
The Section IV presents the whole framework of link prediction method based on
SSNE in detail. Section V contains experimental material and evaluation. The
Section VI presents the experimental result and discusses the effectiveness of
adjustable parameters in link prediction performance. We finally conclude the
paper in the Section VII.

##  II Related works

This section briefly illustrates the related works in two aspects. On the one
hand, we introduce some classical link prediction methods based on structural
similarity index and discuss corresponding research achievements in recent
years. On the other hand, we also discuss some popular graph embedding models
based on representation learning for the link prediction in complex networks.

###  II-A Link prediction based on structural similarity index

The structural similarity indexes are defined by the similarity between two
nodes based on their corresponding local or global structural measurements.
The common structural similarity indexes include the CN index, the AA index,
the RA index, the Katz index, the RWR index, the SimRank index _etc_ (refer to
[7] for more indexes). The CN index calculates the similarity of a pair of
nodes by counting the number of their common neighbors. The AA index [11] and
RA index [12], based on the CN index, punish the common neighbors with large
degrees respectively by the inverse logarithm and the reciprocal of common
neighbors’ degrees. The Katz index [13] can distinguish different influences
of different neighboring nodes and assign different weights to the neighboring
nodes, causing the weights to be biased towards short paths. The RWR index is
an alternative method of the PageRank algorithm and is defined by the
probability of random walk between a pair of nodes in the steady state [14].
The SimRank index [15] also involves with the random walk process and measures
how soon two random walkers respectively starting from the endpoints are
expected to meet a certain node. These structural similarity indexes have been
widely used to infer linking probability for link prediction in complex
networks. Herein, we illustrate some very recent works on the link prediction
methods based on structural similarity indexes.

Inspired by the above common structural similarity indexes, sever recent works
synthesized more structural measurements to form complicated structural
similarity index for link prediction. In [28], Zhou _et al._ replaced the
degree with H index to form H-index-based link prediction methods, which
significantly improve link prediction accuracy. In [29, 30], Zhu et al.
discussed the roles of degree, H-index, and coreness in link prediction in
complex networks, and form a hybrid similarity index synthesizing degree and
H-index. In [31], Pech et al. proposed a simple assumption that the likelihood
of the existence of a link between two nodes can be unfolded by a linear
summation of neighboring nodes’ contributions, and obtained the optimal
likelihood matrix that shows remarkably better predicting prediction
performance.

###  II-B Link prediction based on graph embedding

Graph embedding is used to map network structure into low-dimensional vector
space that is indicated by a feature matrix of nodes reduced from an adjacency
matrix. Based on the feature matrix of nodes, the similarity index is defined
by the similarity between feature vectors of nodes. Thus, the link prediction
method based on graph embedding strongly correlates with the graph embedding
models. Herein, we illustrate some previous works on embedding graph models.

In recent years, graph embedding models have attracted more attention. In
[23], Perozzi _et al._ proposed DeepWalk, where the random walk sampling
processes produce linear sequences of nodes and these sequences are used to
calculate co-occurrence probabilistic matrix of nodes that is mapped into a
feature matrix by a neural network model. In [32], Tang et al. explicitly
defined two objective functions, 1st-order proximity and 2nd-order proximity,
to obtain the topological information of network structure. They then used the
linear combination of such proximity to represent the global proximity. In
[33], Grover et al. proposed Node2Vec, which maintained the high order
proximity between nodes by maximizing the probability of subsequent nodes in
the random traversal graph. Compared with DeepWalk, Node2Vec has made some
improvements in random walk sampling process.

##  III Problem Definition

###  III-A Pre-Definition

An undirected unweighted network is represented by $G=\langle V,E\rangle$
where the node set V={ $v_{1}$,$v_{2}$,…,$v_{n}$ } and the edge set
$E=\left\\{e_{i,j}\right\\}(i,j\in V)$. The dynamic evolution of network
structure is represented by multiple snapshots of network, i.e.,
$\mathbb{G}=\left\\{G_{1},G_{2},\dots,G_{t-1},G_{t},\dots,G_{N-1},G_{N}\right\\}$.
At current time $t$, $G_{t}=<V_{t},E_{t}>$ denotes a snapshot of network.
Naturally, $G_{t-1}=<V_{t-1},E_{t-1}>$ shows a previous snapshot of network.
W͡e assume that the node set is stable, _i.e.,_
$V_{1}=V_{2}=\cdots=V_{t-1}=V_{t}$, but the edge set is dynamically changing,
which characterizes the dynamic evolution of network structure.

For simplicity, any two different nodes are indicated by symbols $u$ and $v$,
and the adjacency matrix of network is expressed by symbol $A$. Obviously, if
there exists an edge between nodes $u$ and $v$, $A(u,v)=1$, otherwise
$A(u,v)=0$. For a node $u$, its adjacency vector is $A_{u}^{1\times
n}=A(u,:)$. We assume that the feature matrix $R$ for node representations is
obtained from the dimensionality reduction of adjacency matrix. In the similar
way, for a node $u$, its $d$-dimensionality vector is $R_{u}^{1\times
d}=R(u,:)$.

We illustrate the important symbols involved in the model of SSNE. In the
SNHAM algorithm, the output is defined by matrix $SNHAM\in\mathbb{R}^{n\times
n}$ and the specific order is set as $h$. Because the elements of SNHAM matrix
can reflect the co-occurrence probability for each pair of nodes, a node $u$
has $n$-dimension vector of co-occurrence probability, $SNHAM_{u}^{1\times
n}=SNHAM(u,:)$. In the single-hidden layer feedforward neural network model,
the input is defined by the matrix $X\in\mathbb{R}^{n\times n}$, and the
kernel and activation function between the input and hidden layers is
respectively set as $W_{1}$ and $f_{1}(x)$; the output is defined by the
matrix $Y\in\mathbb{R}^{n\times n}$, and the kernel and activation function
between the hidden and output layers is respectively set as $W_{2}$ and
$f_{2}(x)$; the feature matrix $R\in\mathbb{R}^{n\times d}$ of node
representations is obtained in the hidden layer. The more details of symbol
description are summarized in Table 1.

TABLE I: Notation note of each symbol Symbol | Comments
---|---
$u$ or $v$ | denote a node
$n$ | the number of nodes, i.e., the full dimension
$d$ | the target of dimensionality reduction
$h$ | the orders of SNHAM algorithm
$t$ | the current time
$t-1$ | the previous time
$G_{t}$ | a snapshot of network at time $t$
$V_{t}$ | the node set in network $G_{t}$
$E_{t}$ | the edge set in the network $G_{t}$
$W_{i}$ | the kernel of $i$th-layer network
$f_{i}(x)$ | the activation function of $i$th-layer network
$\mathbb{G}$ | multiple snapshots of network
$A$ | the adjacency matrix of network
$A_{u}^{1\times n}$ | the adjacency vector for node $u$
$SNHAM$ | the output matrix of the SNHAM algorithm
$SNHAM_{u}^{1\times n}$ | the SNHAM vector for node $u$
$R$ | the feature matrix of node representations
$R_{u}^{1\times d}$ | the $d$-dimension vector for node $u$
$X$ | the input matrix of neural network model
$X_{u}^{1\times n}$ | the input vector for node $u$
$Y$ | the output matrix of neural network model
$Y_{u}^{1\times n}$ | the output vector for node $u$

###  III-B Problem Statement

Obviously, studying the whole dynamic evolution of network $\mathbb{G}$ is a
complicated and challenging task. In order to simplify the process of the
derivation, we herein only consider the relationship between the current $t$
and previous time $t-1$, that is $\mathbb{G}=\left\\{G_{t-1},G_{t}\right\\}$.
Therefore, inferring dynamic evolution of network structure from $t-1$ to $t$
realized by the link prediction based on $G_{t-1}$ and $G_{t}$. The train set
and test set can be set by $G_{t-1}$ and $G_{t}-G_{t-1}$, respectively. Note
that the real (benchmark) networks in the testing experiments aren’t temporal
(i.e., absent of the time scale). We thus assume the original network as
$G_{t}$ and hide a proportion of its edges to assume the residual network as
$G_{t-1}$. Based on $G_{t-1}$, our task is to get the feature matrix of node
representations which meets the lower dimension, but involves with large
number of topological information of network structure, and then apply the
feature matrix for predicting the hidden edges.

##  IV SSNE for link prediction

In this section, we describe the model of SSNE in details. As shown in Figure
1, the SSNE consists of two steps. First, we introduce the SNHAM algorithm to
obtain its corresponding matrix $SNHAM$ that can reflect the theoretical
values of co-occurrence probability between each pair of nodes. Then, we
design a neural network model to calculate the corresponding co-occurrence
probability (i.e., the output matrix $Y$). According to the difference between
the matrices $SNHAM$ and $Y$, the error of loss function is established. By
using the stochastic gradient descent approach to minimize the loss function,
we can get the optimal kernels and determine the feature matrix $R$ in the
hidden layer. However, the stochastic gradient descent approach has high time
complexity in its iterative operation. We then find an alternative method to
directly map the $\log(\mathrm{SNHAM})$ matrix into the feature matrix of node
representations, and demonstrate that the mapping operation is an equivalent
variation of SVD. Finally, we apply the results of the SSNE to construct the
similarity index for link prediction.

Figure 1: (Color online) The schematic description of the SSNE.

###  IV-A SNHAM matirx

Although using the random walk sampling process is effective to converting the
topological information of network structure into liner sequences, it has been
found that this method has obvious drawbacks. As mentioned above, graph
embedding models based on random walk need to determine some parameters of the
random walk, such as the sampling length of a random walk and the number of
random walks, so that they are sensitive to such parameters. More importantly,
we can only determine empirically optimal parameters of the random walk (i.e.,
local best at a finite parameter interval) . In further, the finite-length
linear sequences collected by the random walks have large errors in the
process of representing the boundary nodes. Thus, the accumulative errors of
the multiple parameters greatly affect the accuracy of link predication in
complex networks.

To solve the above problem existing in the graph embedding models based on
random walk, we propose the SNHAM algorithm which can directly capture the
topological information of network structure. We label the nodes in the
network and order them to obtain a adjacency matrix $A$. First, we normalize
the adjacency matrix by row to get the 1st-order transition probability matrix
$S_{1}$. The row normalization function is set as $\mathrm{Normal}(X)$, so the
above operation can be expressed by

$S_{1}=\mathrm{Normal}(A).$ (1)

In a similar way, we calculate the $h$-order transition probability matrix
$S_{h}$ by the $h$-order reachable matrix $A^{h}$ ($h=1,2,\cdots,h$),

$S_{h}=\mathrm{Normal}\left(A^{h}\right),\text{where}\
A^{h}=\overbrace{A\times A\cdots\times A}^{h}.$ (2)

Then, we define the $h$-order similar probabilistic co-occurrence matrix
$SPCO_{h}$, which is the sum of probability transition moments with restart.
The restart probability is set as $\alpha$. Thus, $SPCO_{h}$ is described by

$SPCO_{h}=\sum_{i=1}^{h}\left((1-\alpha)S_{i}+\alpha S_{1}\right).$ (3)

We consider that the restart process is excluded (i.e., $\alpha=0$), that is,
$SPCO_{h}$ can be reduced to the following form,

$SPCO_{h}=\sum_{i=1}^{h}S_{i}.$ (4)

Finally, we normalize the rows of $SPCO_{h}$ matrix, and the final result is
denoted as the SNHAM matrix, which can be expressed as follows:

$SNHAM=\mathrm{Normal}\left(SPCO_{h}\right).$ (5)

The SNHAM algorithm can efficiently obtain the locally topological information
of network structure, and effectively solve the drawbacks of random walk
sampling process. As the restart process is excluded, the adjusting parameter
in the SNHAM algorithm is only the order $h$. Obviously, the single parameter
can avoid the accumulative errors of multiple parameters in the random walk
sampling process. At the same time, the SNHAM algorithm to obtain the network
structure is no longer to transform the network structure into linear
sequences of nodes, so that there doesn’t exist the large errors in the
process of representing the boundary node of each linear sequence. We show the
pseudocode of the SNHAM algorithm in  Algorithm 1.

Algorithm 1 SNHAM

1:adjacency matrix $A$ of $G_{t-1}$; order index $h$

2:SNHAM matrix $SNHAM$

3:Initializing an $n\times n$ matrix $SPCO$;

4:for each $i\in[1,h]$ do

5: Calculating $A^{i}=\overbrace{A\times A\cdots\times A}^{i}$;

6: Normalizing matrix $A^{i}$ by row, $S_{i}=$ Normal $\left(A^{i}\right)$;

7: $SPCO_{h}$=$\sum_{i=1}^{h}$ $S_{i}$

8:end for

9:Normalizing matrix $SPCO_{h}$ by row, SNHAM= Normal $(SPCO_{h})$;

###  IV-B Neural network model

Neural network is widely used to study multi-level feature representation, and
the results obtained from representation learning are proved to be
successfully applied in different fields. Herein, we use the single-hidden
layer feedforward neural network model to construct high-quality and low-
dimensional feature representation based on the SNHAM matrix because it is
assumed to be a potential nonlinear mapping relationship between the vector
representation space of SNHAM matrix and the low-dimensional feature
representation space. In SSNE, the single-hidden layer feedforward neural
network model based on the SNHAM matrix is designed to realize the calculation
of co-occurrence probability matrix.

Specifically, for a given node $u$, we input its one-hot coding vector
$X_{u}^{1\times n}$, and map $X_{u}^{1\times n}$ into a low-dimensional vector
$R_{u}^{1\times d}$ through the kernel $W_{1}$ and activation function
$f_{1}(X)=\mathrm{X}$,

$R_{u}^{1\times d}=f_{1}\left(X_{u}^{1\times n}\cdot W_{1}^{n\times
d}\right)=X_{u}^{1\times n}\cdot W_{1}^{n\times d}.$ (6)

Referring to the neural network model, the low-dimensional vector
$R_{u}^{1\times d}$ is able to be mapped into the co-occurrence probabilistic
vector $Y_{u}^{1\times n}$ through the kernel $W_{2}$ and activation function
$f_{2}(X)=\operatorname{Softmax}(\mathrm{X})=\frac{e^{x}}{\Sigma_{i}e^{i}}$,

$Y_{u}^{1\times n}=\operatorname{Softmax}\left(R_{u}^{1\times d}\cdot
W_{2}^{d\times n}\right).$ (7)

We use the theoretical co-occurrence probabilistic vector $SNHAM_{u}^{1\times
n}$ of node $u$ obtained from the SNHAM matrix, and compare it with
$Y_{u}^{1\times n}$ via the Euclid measurement. The loss function
${L}\left(W_{1},W_{2}\right)$ is built by summing the errors across over all
nodes,

${L}\left(W_{1},W_{2}\right)=\sum_{u\in
V}\mathrm{d}\left(SNHAM_{u},\operatorname{Softmax}\left(\left(X_{u}\cdot
W_{1}\right)\cdot W_{2}\right)\right).$ (8)

By minimizing the loss function, the kernels $W_{1}$ and $W_{2}$ are obtained
through the stochastic gradient descent approach. We focus the low-dimensional
feature matrix in the hidden layer, which is described by

$\text{R}^{n\times d}=X^{n\times n}\cdot W_{1}^{n\times n}.$ (9)

As the stochastic gradient descent approach is high computational complexity,
we provide an alternative method in the following subsection to greatly
improve the computational efficiency of obtaining feature matrix.

###  IV-C SVD for log(SNHAM)

The above-mentioned optimization procedure of minimizing the loss function
${L}\left(W_{1},W_{2}\right)$ is actually equivalent to make
$\operatorname{Softmax}\left(\left(X\cdot W_{1}\right)\cdot W_{2}\right)$
approximate $SNHAM$ by adjusting the kernels $W_{1}$ and $W_{2}$. An ideal
situation is that ${L}\left(W_{1},W_{2}\right)=0$, which satisfies the
condition

$SNHAM_{i}=\operatorname{Softmax}\left(\left(X_{i}\cdot W_{1}\right)\cdot
W_{2}\right),\text{where}\ i=1,2,\cdots,n.$ (10)

We further simplify the variable $\left({X_{i}}\cdot W_{1}\right)\cdot W_{2}$.
Since the input matrix $X$ encoded for one-hot form is actually an identity
matrix, we can write $W_{1}\cdot W_{2}$ as the product matrix $W$. Then,
equation (10) can be rewritten as

$\text{S}NHAM_{i}=\operatorname{Softmax}\left(W_{i}\right),\text{where}\
i=1,2,\cdots,n.$ (11)

Supposing equation (11) has an inverse function, $W_{i}$ can be written as

$\text{W}_{i}=\operatorname{Softmax}^{-1}\left(SNHAM_{i}\right).$ (12)

Naturally, the main task turns to obtain such inverse function. We set a input
vector $x_{i}=\left(x_{i,1},x_{i,2},\ldots,x_{i,j},\ldots,x_{i,n}\right)$, and
the output vector via Softmax function is denoted as
$y_{i}=\left(y_{i,1},y_{i,2},\ldots,y_{i,j},\ldots,y_{i,n}\right)$. Without
loss of generality, each value $x_{i,j}$ producing a corresponding $y_{i,j}$
satisfies an equation,

${y}_{i,j}=\operatorname{Softmax}\left(x_{i,j}\right)=\frac{e^{x_{i,j}}}{\sum
e^{x_{i,j}}}.$ (13)

When the input vector is determined, $\sum e^{x_{i,j}}$ is a constant that is
set as $k_{i}$. Obviously, the conditions are satisfied with

$\text{ s.t. }\quad\left\\{\begin{array}[]{c}{\sum_{j=1}^{n}y_{i,j}=1}\\\ \\\
{\sum_{j=1}^{n}e^{x_{i,j}}=k_{i}}\end{array}\right.$ (14)

and used to obtain a variation of equation (13),

$x_{i,j}=\log\left(k_{i}\cdot y_{i,j}\right).$ (15)

Inspired by equation (15), we assume the inverse function with a formula,

$\operatorname{Softmax}^{-1}\left(y_{i,j}\right)=x_{i,j}=\log\left(c_{i}\cdot
y_{i,j}\right).$ (16)

For a certain $x_{i}$, equation (16) is determined only when $c_{i}$ is
constant. In further, we verify the above-mentioned assumption. Equation (16)
is generalized as,

$\left\\{\begin{aligned} x_{i,1}&=\log\left(c_{i}\cdot y_{i,1}\right)\\\
x_{i,2}&=\log\left(c_{i}\cdot y_{i,2}\right)\\\ ......\\\
x_{i,n}&=\log\left(c_{i}*y_{i,n}\right)\end{aligned}\right.$ (17)

which is equivalent to the following formula,

$\left\\{\begin{aligned} e^{x_{i,1}}&=c_{i}\cdot y_{i,1}\\\
e^{x_{i,2}}&=c_{i}\cdot y_{i,2}\\\ ......\\\ e^{x_{i,n}}&=c_{i}\cdot
y_{i,n}\end{aligned}\right.$ (18)

We sum the left and right terms in equation (18), and obtain

$\sum e^{x_{i,j}}=c_{i}\cdot\sum y_{i,j}.$ (19)

According to the conditions in equation (14), we obtain $c_{i}=k_{i}$ from
equation (19), which implies that for a certain $x_{i}$, $c_{i}$ is a
constant. Thus, the specific formula of inverse function is independent of
$c_{i}$. To make it easy to calculate, we set all $k_{i}$ to 1 by assuming the
independence of input vectors, so that the inverse function is simplified as

$\text{ Softmax }^{-1}\left(y_{i,j}\right)=x_{i,j}=\log\left(y_{i,j}\right).$
(20)

Turning to equation (12), it is specified as

$W_{i}=\operatorname{Softmax}^{-1}\left({SNHAM}_{i}\right)=\log\left({SNHAM}_{i}\right).$
(21)

Considering the zero value of co-occurrence probability in the SNHAM matrix,
we uniformly add a very small positive $\sigma$ ($\sigma=10^{-8}$ in the
testing experiments). We finally obtain the inverse function with the formula

$W=\log(SNHAM+\sigma)=\log(SNHAM^{\prime}).$ (22)

Through equation (22), the specific matrix $W$ is also acquired.

We have known $W=W_{1}\cdot W_{2}$, and divide $W$ by SVD to get $W_{1}$,
$W_{2}$ easily. Obviously, the SVD procedure of $\log(SNHAM^{\prime})$ is
approximately equivalent to the optimization procedure of neural network
model. Without loss of generality, we denote the decomposition process as

$\log(SNHAM^{\prime})=U\Sigma V^{T}$ (23)

We choose the first $d$ largest singular value, and approach
$log(SNHAM^{\prime})$ to $log(SNHAM^{\prime})_{d}$

$\log({SNHAM^{\prime}})\approx\log({SNHAM^{\prime}})_{d}=U_{d}\Sigma_{d}V_{d}^{T}$
(24)

According to equation (24), we easily obtain $W_{1}=U_{d}\Sigma_{d}$ and
$W_{2}=V_{d}^{T}$. Finally, according to equation (9), the $d$-dimensional
feature matrix $R$ can be expressed

$\text{R}^{n\times d}=U_{d}\Sigma_{d}.$ (25)

###  IV-D Similarity index based on feature matrix

After the original network topology is represented by the $d$-dimension
feature matrix by the SSNE, we use such feature matrix to construct similarity
index for link prediction. For any unknown edge $e_{u,v}$ between a pair of
nodes $u$ and $v$, its potential probability is quantified by the similarity
index of these two node. The similarity measurement is used by the Euclidean
distance between the feature vectors of $u$ and $v$, which is described as

$D(e_{u,v})=\sqrt{\left(x_{1v}-x_{1u}\right)^{2}+\left(x_{2v}-x_{2u}\right)^{2}+\cdots+\left(x_{dv}-x_{du}\right)^{2}}$
(26)

Considering the inverse correlation that the greater the distance is, the
lower the similarity is, we take its reciprocal and add 1 to $D(e_{u,v})$ to
prevent the case that $D(e_{u,v})$ is zero or too small. Finally, the
similarity index is constructed by

$S(e_{u,v})=\frac{1}{1+D(e_{u,v})}.$ (27)

In the link prediction in complex networks, it thinks that higher the
similarity index is, the higher potential possibility the unknown edge will be
linked with. We show the pseudocode of the link prediction method based on
SSNE in  Algorithm 2.

Algorithm 2 Link prediction based on SSNE

1:SNHAM Matrix $SNHAM$; dimension $d$

2:Evaluation Index $AUC$

3:Calculating $\log(\mathrm{SNAHM}+\sigma)$;

4:$\mathrm{U}\Sigma V^{T}=SVD(\log(SPMI+\sigma))$

5:Choosing $d$ largest singular value,
$\mathrm{U}_{d}\Sigma_{d}V_{d}^{T}\approx\mathrm{U}\Sigma V^{T}$

6:Obtaining feature matrix R, $R^{n\times d}=U_{d}\Sigma_{d}$

7:Calculating Euclidean distance,

8:$D(e_{u,v})=\sqrt{\left(x_{1v}-x_{1u}\right)^{2}+\cdots+\left(x_{dv}-x_{du}\right)^{2}}$

9:Calculating similarity index, $S(e_{u,v})=\frac{1}{1+D(e_{u,v})}$

10:Initializing sampling parameter of AUC, $N=672400$, $N^{\prime}=0$,
$N^{\prime\prime}=0$;

11:for each $i\in[1,N]$ do

12: if $S_{i}\left(e_{a,b}\right)>S_{i}\left(e_{c,d}\right)$ then
$N^{\prime}+=1$

13: else if $S_{i}\left(e_{a,b}\right)=S_{i}\left(e_{c,d}\right)$ then
$N^{\prime\prime}+=1$

14: else other cases

15: end if

16:end for

17:Calculating AUC, $AUC=\frac{N^{\prime}+0.5\cdot N^{\prime\prime}}{672400}$

##  V Experimental material and evaluation

We design testing experiments based on six real networks and two types of
artificial network models to validate the effectiveness of SSNE for the link
prediction in complex networks. In this section, the specific descriptions of
real networks, two types of artificial network models and the evaluation are
illustrated, respectively.

###  V-A Real networks

We show six real networks that are described as:

Brain[34]: It is the neuronal connection network of macaques or rhesus. The
nodes and edges represent neurons and fiber bundles among these neurons,
respectively. In this network, there are 242 neurons, 3,054 fiber bundles, and
the average degree of the network is 25.24.

Yeast[35]: It is the biological network in budding yeast. The nodes and edges
represent proteins and interactions among these proteins. In this network,
there are 2,375 proteins and 11,693 known interactions, and the average degree
of network is 9.85.

Air[34]: It is the traffic network of air control. The nodes and edges
represents airports or service centers and preferred air route among these
airports or service centers recommended by the NFDC (National Flight Data
Center), respectively. In this network, there are 1,226 airports or services
centers, 2,410 preferred air routes, and the average degree of the network is
3.93

Road[34]: It is the road network in Minnesota state. The nodes and edge
represents voluntary buildings and direct road between these buildings. In
this network, there are 2,642 buildings, and there are 3,303 direct roads, and
the average degree of the network is 2.50.

Twitter[36]:It is the forwarding network of twitter users about Obama’s re-
election as President of the United States in 2012. The nodes and edges
represent twitter users and retweeting relationships between these users,
respectively. In this network, there are 3,212 twitter users, 3,423 retweeting
relationships, and the average degree of the network is 2.13.

Power[34]:It is the west power grid in the U.S. The nodes and edges represent
substations or converters and high-voltage lines among these substations or
converters. In this network, there are 4,941 substations or converters, 6,594
high-voltage lines, and the average degree of the network is 2.70.

We also summarize the basic topological information of six real networks,
including the number of nodes and edges, the edge sparsity, the average
degree, the clustering coefficient and the degree heterogeneity, which are
shown in Table 2.

TABLE II: Basic topological information of six real networks. $|V|$ and $|E|$ indicate the number of nodes and edges, $\langle k\rangle$ is the average degree, $ES$ is the edge sparsity, $\langle d\rangle$ is the average distance, $C$ is the clustering coefficient, and $H=\frac{\langle k^{2}\rangle}{{\langle k\rangle}^{2}}$ is the degree heterogeneity. Nets | $|V|$ | $|E|$ | $\langle k\rangle$ | $ES$ | $\langle d\rangle$ | $C$ | $H$
---|---|---|---|---|---|---|---
$Brain$ | 242 | 3,054 | 25.24 | 0.1047 | 2.22 | 0.450 | 1.53
$Yeast$ | 2,375 | 11,693 | 9.85 | 0.0041 | 5.09 | 0.306 | 3.48
$Air$ | 1,226 | 2,410 | 3.93 | 0.0032 | 5.92 | 0.068 | 1.88
$Road$ | 2,642 | 3,303 | 2.50 | 0.0009 | 35.35 | 0.016 | 1.09
$Twitter$ | 3,212 | 3,423 | 2.13 | 0.0006 | 7.31 | 0.004 | 19.16
$Power$ | 4,941 | 6,594 | 2.67 | 0.0005 | 18.99 | 0.080 | 1.45

###  V-B Artificial network models

We have known that the BA and WS networks models are widely used to simulate
real complex networks because they characterizes the stylized facts of real
complex networks. Herein, we show two types of artificial network models that
are used in the following research, which are described as:

Barabasi-Albert network model[37]: The BA network model proposed by Barabasi
and Albert characterizes the scale-free property of real complex networks. By
using mean filed approximation, it can be proved that the resulted BA network
has a power-law degree distribution with an scale exponent 3. In the
simulating process, the number of nodes and edges are adjustable according to
the actual need.

Watts-Strogatz network model[38]: The WS network model proposed by Watts and
Strogatz characterizes the small-world property of real complex networks. The
resulted WS network has a larger cluster coefficient and shorter average
distance, however its degree distribution is Poisson. In the simulating
process, the number of nodes and edges, as well as the rewired probability are
adjustable according to the actual need.

###  V-C Evaluation

The common measuring index for evaluating link prediction method is AUC that
refers to the area under the curve of the receiver operating characteristic
(ROC) [39]. In the calculation of AUC, we don’t draw the specific ROC curve,
especially when the samples are very large. Alternatively, we generally use
the sampling method to obtain its approximate value. Once the partition of
training set and the testing set is determined, there are two kinds of unknown
edges in the training set, that is, one corresponds to the nonexistent edges
(i.e., they doesn’t exist in both training and testing sets) and the other
corresponds to the hidden edges (i.e., they only exist in the testing set).
For a given link prediction method, each unknown edge is given a similarity
index. AUC is equivalent to the probability that the similarity index of
randomly selected hidden edge in the testing set is higher than that of
randomly selected nonexistent edge [40].

That is to say, we randomly select a hidden edge and a nonexistent edge in the
testing set. If the similarity index of hidden edge is higher than that of
nonexistent edge, the AUC value is added by $1$, and if these two similarity
indexes are equal, the AUC value is added by $0.5$. The sampling process is
repeated with $N$ times. We assume that there are $N^{\prime}$ and
$N^{\prime\prime}$ times of the sampling processes that meet the above-
mentioned two cases, respectively. The AUC value is calculated as,

$AUC=\frac{N^{\prime}+0.5\cdot N^{\prime\prime}}{N}$ (28)

According to [41], when $N\geq 672400$, we can guarantee with 90% confidence
that the absolute error of AUC will not exceed one thousandth.

##  VI Experimental result and discussion

In this section, we mainly present the performance of link prediction method
based on SSNE, and compare the proposed method with other baselines. The
20$\%$ edges of current network $G_{t}$ is hidden to obtain the previous
network $G_{t-1}$. There are no isolated nodes in both $G_{t-1}$ and $G_{t}$.
Furthermore, we explore the effectiveness of adjustable parameters in the
proposed method according to the experimental results based on real networks
and artificial networks. Finally, we summarize the optimal AUC values obtained
from the proposed method and the mainstream methods based on six real networks
and two types of artificial network models.

###  VI-A Link prediction in real networks

Herein, we first examine the performance of link prediction method based on
SSNE, and compare the proposed method with several mainstream methods based on
structural similarity index and graph embedding, such as CN, AA, RA, RWR and
DeepWalk. More other methods are shown in the following summary of
experimental result. AUC is used to evaluate the link prediction performance
of these methods. The order $h$ and dimension $d$ are considered as the
adjustable parameters, which regulate the link prediction method based on
SSNE. Because the full dimension $n$ is different from each network, $d$ is
dependent on $n$, _i.e._ , $d=p\times n$. Note that $p$ is an alternative
parameter of $d$ that indicates the proportion of dimension reduction to
network size. Figure 2 presents the performance comparison of different link
prediction methods for six real networks. It suggests that except the Yeast,
the link prediction method based on SSNE (short of $SSNE(h,p)$) behaves better
than these mainstream methods.

Figure 2: (Color online) Performance comparison of different link prediction
methods for six real networks. Except the Yeast, the link prediction method
based on SSNE (short of $SSNE(h,p)$) behaves better than these mainstream
methods.

More concretely, as shown in Figure 2, it is found that in these network with
relatively large average degree (e.g., Brain, Yeast), the link prediction
performance of the proposed method is similar to that of method based on
DeepWalk, and both of them do not significantly outperform other methods based
on structural similarity index. However, when the average degree is relatively
small (e.g., Road, Twitter, Power), the proposed method performs the best.
Thus, we think that the proposed method is more suitable to solve the link
prediction problem of sparse network. Note that such observation will be
further verified by the artificial networks in the following subsection.

Figure 3: (Color online) The Influence of both the order $h$ and the
proportion $p$ on the link prediction method based on SSNE. For each network,
we show the AUC values in respect to different combinations of $h$ and $p$.
The results verify the stability of the proposed method because they share a
similar trend in respect to $h$ and $p$.

At the same time, it is also found that the proposed method is affected by the
adjustable parameters. We use different combinations of order $h$ and
proportion $p$ to comprehensively analyze the link prediction for six real
networks. Figure 3 presents the influence of both $h$ and $p$ on the link
prediction performance based on six real networks. The best AUC values of six
real networks are 0.938 of Yeast, 0.856 of Brain, 0.834 of Air, 0.952 of Road,
0.616 of Twitter, and 0.928 of Power. We find that the proposed method is not
particularly sensitive to the changes of $h$ and $p$. More concretely, for a
given $h$, the link prediction performance is nearly unchanged when $p$ varies
from 0.1 to 0.9, which suggests that the size of dimension reduction has
little influence on the proposed method. While for a given $p$, the link
prediction performance changes largely when $h$ gradually increases in a small
range and then becomes approximately stable with the convergence of SNHAM
matrix , which implies that the SNHAM matrix with a small order (at least 10)
contains most of the topological information of network structure.
Furthermore, the results of each network reveal a similar trend, which verify
the stability of the proposed method. After the above-mentioned analysis, we
observe that when $h=10$ and $p=0.1$, the proposed method almost converges to
the optimal link prediction performance. It roughly suggests the default set
of adjustable parameters can correspond to $h=10$ and $p=0.1$ for obtaining
better link prediction performance.

###  VI-B link prediction in artificial networks

We also testify the proposed method based on artificial networks. The
artificial networks are generated by the BA and WS network models,
respectively. Each type is composed of multiple artificial networks with
different average degrees and sizes of nodes. Specifically, the sizes of nodes
in the BA (or WS) networks vary from 1000 to 5000, and for the BA (or WS)
networks with the fixed size, their average degrees vary from 2 to 10 with a
step-length 2 by adding edges which indicates the changes of edge sparsity. We
try to study the relationship between the network sparsity and link prediction
performance (i.e., AUC) obtained from the proposed method.

Figure 4: (Color online) The link prediction performance of the proposed
method based on the BA and WS networks with different average degrees and
sizes of nodes. In the left panel, the AUC values as a function of the average
degree shows that proposed method is much more suitable for the BA network
with relatively smaller average degree and lower edge sparsity. In the right
panel, the AUC values as a function of the average degree shows that the
result is contrary to that found in the BA networks.

Figure 4 presents the AUC values as a function of the average degree, which
are obtained from the link prediction in the BA and WS networks. As shown in
the left panel of Figure 4, we can see that the link prediction performance is
better realized by the proposed method when the BA networks have relatively
smaller average degree and lower edge sparsity (e.g. $\langle k\rangle=2$ and
$N=5000$). In particular, no matter of the sizes of nodes, the AUC values are
optimal when the average degrees of BA networks are $\langle k\rangle=2$,
which suggests that the proposed method is sensitive to the average degree.
Meanwhile, as shown in the right panel of Figure 4, we can see that when the
average degrees in the WS networks increase, the link prediction performance
become much better, which is contrary to the results found in the BA networks.
Nevertheless, when the average degrees increase, the differences of link
prediction performance between the proposed method and those methods based on
structural similarity indexes become smaller (see in Table 3). That is, in the
sparse WS networks ($\langle k\rangle=2$), the proposed method show much
better link prediction performance. Thus, these results to some extent show
that the proposed method is more suitable for the link prediction in sparse
networks.

###  VI-C Discussion

As we comprehensively analyze the link prediction performance of the proposed
method based on real networks and artificial networks, we further discuss the
performance comparison between the proposed method and more mainstream methods
by conducting extensive experiments. Note that the average degrees of
artificial networks are set as 2 and 6, and their sizes are both 5000 nodes.

Table 3 presents the performance comparison of all link prediction methods by
the AUC values. The mainstream methods are divided into three types, each of
them is based on structural similarity indexes including CN, Salton, Jaccard,
AA, RA, RWR LHN-I, LHN-II, Katz, SimRank, matrix optimization including CLMC
[42] and other graph embedding models including DeepWalk (DW), Diff2Vec (D2V)
[43], Struc2Vec (S2V) [44]. The specific parameter sets are illustrated: RWR
with parameter $c=0.8$; SimRank with parameter $\lambda=0.8$; CLMC with
parameters $\alpha_{1}=0.001$, $\alpha_{2}=0.01$ and $\alpha_{3}=100$;
DeepWalk with parameters $Windows=10$, $length=40$, $times=30$, $d=128$;
Diff2Vec with parameters $\alpha=0.025$, $Windows=10$,
$vertexsetcardinality=40$, $numdiffusions=10$, $d=128$; Struc2Vec with
parameters $times=20$, $Windows=5$, $length=40$, $d=64$; SSNE with default
parameters $h=10$, $rate=0.1$.

More concretely, in Table 3, the first and second best AUC values are
underlined by bold characters. We can find that for these networks with
relatively large average degree (e.g., Brain and Yeast the link prediction
performance obtained by structural similarity indexes is better than other two
types of link prediction methods. However, except the second best AUC values
in the Air network, we can see that the proposed method (i.e., SSNE) achieves
the best AUC values for the link prediction in these real networks with
relatively small average degree. And, in artificial networks, we can see that
the proposed method performs the best AUC values for the link prediction in
these artificial sparse networks (i.e., $\langle k\rangle=2$).

TABLE III: The performance comparison of link prediction methods

Nets | CN | Salton | Jaccard | AA | RA | RWR | LHN-I | LHN-II | Katz | SimRank | CLMC | DW | D2V | S2V | SSNE
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$Brain$ | 0.884 | 0.885 | 0.905 | 0.891 | 0.889 | 0.886 | 0.771 | 0.651 | 0.886 | 0.754 | 0.784 | 0.839 | 0.793 | 0.652 | 0.846
$Yeast$ | 0.911 | 0.910 | 0.913 | 0.912 | 0.912 | 0.967 | 0.905 | 0.962 | 0.967 | 0.674 | 0.928 | 0.954 | 0.947 | 0.271 | 0.934
$Air$ | 0.620 | 0.619 | 0.626 | 0.620 | 0.620 | 0.813 | 0.619 | 0.781 | 0.813 | 0.674 | 0.621 | 0.757 | 0.863 | 0.350 | 0.823
$Road$ | 0.514 | 0.514 | 0.515 | 0.514 | 0.514 | 0.791 | 0.514 | 0.779 | 0.779 | 0.931 | 0.202 | 0.764 | 0.940 | 0.529 | 0.942
$Twitter$ | 0.496 | 0.496 | 0.507 | 0.496 | 0.496 | 0.483 | 0.496 | 0.465 | 0.483 | 0.432 | 0.479 | 0.346 | 0.469 | 0.135 | 0.558
$Power$ | 0.614 | 0.614 | 0.614 | 0.614 | 0.614 | 0.809 | 0.614 | 0.809 | 0.809 | 0.919 | 0.505 | 0.860 | 0.907 | 0.450 | 0.927
$BA(2)$ | 0.499 | 0.499 | 0.500 | 0.499 | 0.499 | 0.477 | 0.499 | 0.477 | 0.477 | 0.346 | 0.354 | 0.374 | 0.370 | 0.181 | 0.581
$BA(6)$ | 0.519 | 0.518 | 0.526 | 0.519 | 0.519 | 0.599 | 0.518 | 0.497 | 0.599 | 0.646 | 0.507 | 0.582 | 0.640 | 0.242 | 0.499
$WS(2)$ | 0.501 | 0.501 | 0.501 | 0.501 | 0.501 | 0.498 | 0.501 | 0.498 | 0.498 | 0.551 | 0.516 | 0.645 | 0.693 | 0.561 | 0.713
$WS(6)$ | 0.769 | 0.770 | 0.771 | 0.770 | 0.770 | 0.826 | 0.770 | 0.828 | 0.826 | 0.306 | 0.753 | 0.832 | 0.835 | 0.482 | 0.843

##  VII Conclusion

As graph embedding is recently used for link prediction in complex networks,
this paper proposes a novel link prediction method based on SSNE that is
constructed in the framework of graph embedding. We comprehensively describe
the procedure of SSNE from two aspects, the SNHAM matrix and the neural
network model. The SNHAM matrix contains the $h$-order structural information
of adjacency matrix, and the neural network model is used to learn the
$d$-dimensional representation of the SNHAM matrix. Through the SSNE, we can
effectively obtain the graph representation of network structure. Note that
the graph embedding procedure of SSNE is irrelevant to specific network
structure. Most importantly, in the SSNE, the adjustable parameters have been
extremely reduced into two variables. Thus, the SSNE overcomes the critical
drawbacks of the random-walk-based graph embedding models by avoiding directly
random walk sampling process.

Meanwhile, to reduce the computational complexity of neural network model, we
assume that the optimization procedure of mining the loss function is actually
equivalent to make the output matrix approximate the SNHAM matrix by adjusting
the kernels of neural network model. The output matrix is denoted by the
product matrix of the kernels. Then, we formalize the association between the
SNHAM matrix and the product matrix through the Softmax function. By verifying
the assumption of the inverse Softmax function, we obtain the product matrix
indicated by the logarithmic SNHAM matrix. Finally, we use the SVD to solve
the product matrix and obtain the $d$-dimensional feature matrix.

The link prediction method based on feature matrix is constructed by
calculating the similarity indexes among feature vectors. We use six real
networks and two types of artificial network models to test the link
prediction performance of the proposed method. The testing experiments are
designed in three aspects. We first verify the effectiveness of the proposed
method on the link prediction in diverse real networks and the sensitivity of
the adjustable parameters to the proposed method. It has found that the
proposed method is more suitable for the link prediction in relatively sparse
network and only partially sensitive to the order of the SNHAM matrix.
Then,the effectiveness of the proposed method on the link prediction in sparse
network is further verified based on artificial networks. Finally, we discuss
the results of the comparison of proposed method with a lot of mainstream
methods based on structural similarity indexes, matrix optimization and other
graph embedding models. It suggests that the proposed method show better link
prediction performance in relatively sparse network.

## References

  * [1] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, “Complex networks: Structure and dynamics,” _Physics Reports_ , vol. 424, no. 4-5, pp. 175–308, 2006.
  * [2] E. Bullmore and O. Sporns, “Complex brain networks: graph theoretical analysis of structural and functional systems,” _Nature Reviews Neuroscience_ , vol. 10, no. 3, pp. 186–198, 2009.
  * [3] S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardenes, M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin, “The structure and dynamics of multilayer networks,” _Physics Reports_ , vol. 544, no. 1, pp. 1–122, 2014.
  * [4] L. Getoor and C. P. Diehl, “Link mining: a survey,” _ACM SIGKDD Explorations Newsletter_ , vol. 7, no. 2, pp. 3–12, 2005.
  * [5] D. Liben-Nowell and J. Kleinberg, “The link-prediction problem for social networks,” _Journal of the American Society for Information Science and Technology_ , vol. 58, no. 7, pp. 1019–1031, 2007.
  * [6] R. N. Lichtenwalter, J. T. Lussier, and N. V. Chawla, “New perspectives and methods in link prediction,” in _Proceedings of the 16th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2010, pp. 243–252.
  * [7] L. Lü and T. Zhou, “Link prediction in complex networks: A survey,” _Physica A_ , vol. 390, no. 6, pp. 1150–1170, 2011.
  * [8] V. Martínez, F. Berzal, and J.-C. Cubero, “A survey of link prediction in complex networks,” _ACM Computing Surveys_ , vol. 49, no. 4, p. 69, 2017\.
  * [9] S. Haghani and M. R. Keyvanpour, “A systemic analysis of link prediction in social network,” _Artificial Intelligence Review_ , vol. 52, no. 3, pp. 1961–1995, 2019.
  * [10] F. Lorrain and H. C. White, “Structural equivalence of individuals in social networks,” _The Journal of Mathematical Sociology_ , vol. 1, no. 1, pp. 49–80, 1971.
  * [11] L. A. Adamic and E. Adar, “Friends and neighbors on the web,” _Social Networks_ , vol. 25, no. 3, pp. 211–230, 2003.
  * [12] T. Zhou, L. Lü, and Y.-C. Zhang, “Predicting missing links via local information,” _The European Physical Journal B_ , vol. 71, no. 4, pp. 623–630, 2009.
  * [13] L. Katz, “A new status index derived from sociometric analysis,” _Psychometrika_ , vol. 18, no. 1, pp. 39–43, 1953.
  * [14] S. Brin and L. Page, “The anatomy of a large-scale hypertextual web search engine,” _Computer networks and ISDN systems_ , vol. 30, no. 1-7, pp. 107–117, 1998.
  * [15] G. Jeh and J. Widom, “Simrank: a measure of structural-context similarity,” in _Proceedings of the 8th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2002, pp. 538–543.
  * [16] R. R. Sarukkai, “Link prediction and path analysis using markov chains,” _Computer Networks_ , vol. 33, no. 1-6, pp. 377–386, 2000.
  * [17] J. Zhu, J. Hong, and J. G. Hughes, “Using markov chains for link prediction in adaptive web sites,” in _International Conference on Soft Issues in the Design, Development, and Operation of Computing Systems_. Springer, 2002, pp. 60–73.
  * [18] H. Xiao, M. Huang, and X. Zhu, “From one point to a manifold: Knowledge graph embedding for precise link prediction,” _arXiv preprint arXiv:1512.04792_ , 2015.
  * [19] R.-M. Cao, S.-Y. Liu, and X.-K. Xu, “Network embedding for link prediction: The pitfall and improvement,” _Chaos_ , vol. 29, no. 10, p. 103102, 2019\.
  * [20] Y. Gu, Y. Sun, Y. Li, and Y. Yang, “Rare: Social rank regulated large-scale network embedding,” in _Proceedings of the 27th International Conference on World Wide Web Conference_. ACM, 2018, pp. 359–368.
  * [21] J. Ni, S. Chang, X. Liu, W. Cheng, H. Chen, D. Xu, and X. Zhang, “Co-regularized deep multi-network embedding,” in _Proceedings of the 27th International Conference on World Wide Web Conference_. ACM, 2018, pp. 469–478.
  * [22] G. H. Nguyen, J. B. Lee, R. A. Rossi, N. K. Ahmed, E. Koh, and S. Kim, “Continuous-time dynamic network embeddings,” in _Companion Proceedings of the The Web Conference 2018_. ACM, 2018, pp. 969–976.
  * [23] B. Perozzi, R. Al-Rfou, and S. Skiena, “Deepwalk: Online learning of social representations,” in _Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2014, pp. 701–710.
  * [24] G. Li, J. Luo, Q. Xiao, C. Liang, P. Ding, and B. Cao, “Predicting microrna-disease associations using network topological similarity based on deepwalk,” _IEEE Access_ , vol. 5, pp. 24 032–24 039, 2017.
  * [25] Z. Chen, X. Wang, P. Gao, H. Liu, and B. Song, “Predicting disease related microrna based on similarity and topology,” _Cells_ , vol. 8, no. 11, p. 1405, 2019.
  * [26] Z. Jin, R. Liu, Q. Li, D. D. Zeng, Y. Zhan, and L. Wang, “Predicting user’s multi-interests with network embedding in health-related topics,” in _2016 International Joint Conference on Neural Networks_. IEEE, 2016, pp. 2568–2575.
  * [27] Z. Yu, Z. Zhang, H. Chen, and J. Shao, “Structured subspace embedding on attributed networks,” _Information Sciences_ , vol. 512, pp. 726–740, 2020\.
  * [28] W. Zhou, J. Gu, and Y. Jia, “h-index-based link prediction methods in citation network,” _Scientometrics_ , vol. 117, no. 1, pp. 381–390, 2018.
  * [29] X. Zhu, Y. Yang, L. Li, and S. Cai, “Roles of degree, h-index and coreness in link prediction of complex networks,” _International Journal of Modern Physics B_ , vol. 32, no. 16, p. 1850197, 2018.
  * [30] X. Zhu, W. Li, H. Tian, and S. Cai, “Hybrid influence of degree and h-index in the link prediction of complex networks,” _EPL_ , vol. 122, no. 6, p. 68003, 2018.
  * [31] R. Pech, D. Hao, Y.-L. Lee, Y. Yuan, and T. Zhou, “Link prediction via linear optimization,” _Physica A_ , vol. 528, p. 121319, 2019.
  * [32] J. Tang, M. Qu, M. Wang, M. Zhang, J. Yan, and Q. Mei, “Line: Large-scale information network embedding,” in _Proceedings of the 24th International Conference on World Wide Web_. ACM, 2015, pp. 1067–1077.
  * [33] A. Grover and J. Leskovec, “node2vec: Scalable feature learning for networks,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2016, pp. 855–864.
  * [34] J. Kunegis, “Konect – the koblenz network collection,” 2013.
  * [35] V. Batagelj and A. Mrvar, “Pajek datasets,” 2009.
  * [36] R. A. Rossi, D. F. Gleich, A. H. Gebremedhin, and M. A. Patwary, “What if clique were fast? maximum cliques in information networks and strong components in temporal networks,” _arXiv preprint arXiv:1210.5802_ , pp. 1–11, 2012.
  * [37] A.-L. Barabási and R. Albert, “Emergence of scaling in random networks,” _Science_ , vol. 286, no. 5439, pp. 509–512, 1999.
  * [38] D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’ networks,” _Nature_ , vol. 393, no. 6684, p. 440, 1998.
  * [39] J. A. Hanley and B. J. McNeil, “The meaning and use of the area under a receiver operating characteristic (roc) curve.” _Radiology_ , vol. 143, no. 1, pp. 29–36, 1982.
  * [40] T. Fawcett, “An introduction to roc analysis,” _Pattern Recognition Letters_ , vol. 27, no. 8, pp. 861–874, 2006.
  * [41] Q.-M. Zhang, X.-K. Xu, Y.-X. Zhu, and T. Zhou, “Measuring multiple evolution mechanisms of complex networks,” _Scientific Reports_ , vol. 5, p. 10350, 2015.
  * [42] Z. Zhang, D. Kang, C. Gao, and J. Shao, “Semisync: Semi-supervised clustering by synchronization,” in _International Conference on Database Systems for Advanced Applications_. Springer, 2019, pp. 358–362.
  * [43] B. Rozemberczki and R. Sarkar, “Fast sequence-based embedding with diffusion graphs,” in _International Workshop on Complex Networks_. Springer, 2018, pp. 99–107.
  * [44] L. F. Ribeiro, P. H. Saverese, and D. R. Figueiredo, “struc2vec: Learning node representations from structural identity,” in _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2017, pp. 385–394.


# Disentangling the spectral properties of the Hodge Laplacian:
Not all small eigenvalues are equal

###### Abstract

The rich spectral information of the graph Laplacian has been instrumental in
graph theory, machine learning, and graph signal processing for applications
such as graph classification, clustering, or eigenmode analysis. Recently, the
Hodge Laplacian has come into focus as a generalisation of the ordinary
Laplacian for higher-order graph models such as simplicial and cellular
complexes. Akin to the traditional analysis of graph Laplacians, many authors
analyse the smallest eigenvalues of the Hodge Laplacian, which are connected
to important topological properties such as homology. However, small
eigenvalues of the Hodge Laplacian can carry different information depending
on whether they are related to curl or gradient eigenmodes, and thus may not
be comparable. We therefore introduce the notion of persistent eigenvector
similarity and provide a method to track individual harmonic, curl, and
gradient eigenvectors/-values through the so-called persistence filtration,
leveraging the full information contained in the Hodge-Laplacian spectrum
across all possible scales of a point cloud. Finally, we use our insights (a)
to introduce a novel form of topological spectral clustering and (b) to
classify edges and higher-order simplices based on their relationship to the
smallest harmonic, curl, and gradient eigenvectors.

Index Terms—  Spectral Clustering, Topological Signal Processing, Persistent
Homology, Topological Data Analysis

##  1 Introduction

In many areas of Mathematics, Physics, and Signal Processing, the eigenvalues
of a linear operator or matrix encode important properties of the system
governed by the operator. Accordingly, harmonic and spectral analysis
approaches have been proven to be useful in many application areas. In the
context of the analysis of networks and graphs, the idea of exploiting the
spectral properties of an algebraic representation of the system is deeply
engrained in spectral graph theory [1, 2]. While the initial focus of spectral
graph theory has historically been on the adjacency matrix, the combinatorial
and normalized Laplacian matrices (and variants thereof) have replaced the
adjacency matrix as the operator typically considered, which may be explained
by their favourable structure (e.g., positive semidefiniteness) and their deep
connections to continuous mathematics. Intuitively, the Graph Laplacian is an
operator that encodes a notion of “local difference” between a node and its
neighbours and may be seen as a discrete equivalent of the continous Laplacian
operator. Its $0$-eigenvalues encode connected components, whereas remaining
smaller eigenvalues carrying information on graph connectivity [3] are used
for spectral clustering [4] and graph classification and matching [5]. Further
examples of spectral characteristics that are informative for applications
include the spectral radius of the adjacency matrix, which governs the
epidemic spreading threshold [6], the second eigenvalue of the Laplacian,
which characterizes the convergence of consensus dynamics [7], and the
spectral gap of the transition matrix, which is paramount for the mixing
behaviour of Markov chains [8].

In these archetypal examples, the relevant feature of an eigenvalue is its
numerical value; eigenvalues are “equal” in other regards. In this paper we
will argue that this mantra does not hold for the spectrum of the Hodge
Laplacian, a higher order generalisation of the Graph Laplacian for simplicial
complexes (SCs). SCs are higher-order generalisations of graphs which allow to
incorporate complex relationships and geometric and topological information
beyond pairwise interactions. SCs have enjoyed a renaissance in the context of
topological signal processing in the last years [9, 10]. Importantly, the set
of eigenvector/eigenvalue pairs of the Hodge Laplacian consists of three
distinct types of eigenvectors: _harmonic_ , _gradient_ , and _curl_
eigenvectors, which correspond to the Hodge decomposition of the simplicial
signal spaces. Accordingly, rather than treating the spectrum of the Hodge
Laplacian simply as a set of numbers whose numerical values are of interest,
in this paper we more explicitly consider the gradient, curl and harmonic
eigenvalues as separate to leverage the rich additional information provided
by the Hodge Laplacian spectrum.

Contributions and Outline. In Section 2, we give a brief overview over the
concepts of topological signal processing, simplicial complexes,
$\alpha$-filtrations, and persistent Laplacians. In Figure 2, we track and
analyse the evolution of the eigenvectors of the persistent Laplacian across
the stages of an $\alpha$-filtrations theoretically and experimentally. We
introduce the concepts of persistent eigenvector similarity and matchings to
give a new form of persistent eigenvalue diagrams. We finally apply these
insights on data sets in Figure 3, developing the notions of _topological
spectral clustering_ (Figure 3) and the HGC-values for inferring simplex roles
(Figure 4).

Related Work. The persistent Laplacians were first defined in [11]. However,
the authors did not distinguish between the gradient and curl eigenvalues in
their analysis. Further work on the persistent Laplacians includes [12, 13,
14]. The Hodge Laplacian has already been studied in the topological signal
processing literature [10]. Ideas of spectral clustering [4] have been applied
in some form to simplicial complexes. However, Refs. [15, 16] only consider
harmonic eigenvectors, Ref. [17] is only interested in the support of the
eigenvectors and disregards the gradient/curl differences, and Refs. [18, 19]
are only interested in point clustering. The gradient and curl spaces of the
Hodge Laplacian have been explored in [20] without the analysis of the
connection to filtrations of simplicial complexes.

Fig. 1: Smallest eigenvectors at three stages of the $\alpha$-filtration.
_Blue:_ Harmonic, _Green:_ Gradient, _Red:_ Curl eigenvectors. The intensity
of the colour denotes the absolute value of the corresponding entry in the
eigenvector, where we scale the maximal entry to 1 for visualisation purposes.
The three stages of the filtration correspond to the rows and represent three
different regimes: 1st row: There are many small “holes” in the data set.
Accordingly, the first five eigenvectors are harmonic and strongly localised.
In the subsequent eigenvalues, gradient components dominate, which are being
localised on the bridges between the clusters. The curl eigenvalues are large
in comparison. 2nd row: There is a single harmonic eigenvector corresponding
to the hole in the middle of the clusters, which is stronger on the bridges
and the edges closer to the centre. The $\updownarrow$ and $\leftrightarrow$
gradient vectors are next and are localised at the bridges. We have a curl
flow for each of the four clusters, and finally a diagonal gradient flow. 3rd
row: There is a $1:1$ correspondence to the types of the previous vectors.
However, the order of the eigenvectors has changed with the curl eigenvalues
getting smaller. Both gradient and curl vectors are less localised.

##  2 Simplicial complexes and persistence

Simplicial Complexes. Simplicial complexes [21, 22] are generalisations of
graphs consisting of nodes ($0$-simplices), edges ($1$-simplices), filled-in
triangles ($2$-simplices), tetrahedra ($3$-simplices), etc.

######  Definition 2.1 (Simplicial complex).

A simplicial complex $\mathcal{S}$ consists of a set of vertices $X$ and a set
of finite non-empty subsets of $X$, called simplices $S$, such that (i) $S$ is
closed under taking non-empty subsets and (ii) for every $x\in X$, the
singleton $\\{x\\}$ is contained in $S$. For simplicity, we often identify
$\mathcal{S}$ with its set of simplices and use $\mathcal{S}_{n}$ to denote
the subset of simplices with $n+1$ elements.

The orientation of an simplex $\sigma\in\mathcal{S}$ is represented by an
ordering of its vertices. We consider two orientations equivalent if they
differ by an even permutation. Hence, for $k>0$ every $k$-simplex admits two
possible orientations. We denote by
$C_{k}\coloneqq\mathbb{R}^{|\mathcal{S}_{k}|}$ the $k$-th simplicial signal
space, i.e. the space of real signals supported on the $k$-simplices. An
orientation on the simplices allows for bookkeeping in form of _boundary
matrices_ : The boundary matrix $B_{k}(\mathcal{S})$ of an SC records the
incidence relations between the $k$-simplices and the $(k-1)$-simplices with
respect to the chosen orientations. For $k=1$, this is just the ordinary edge-
vertex incidence matrix of a graph. We define the boundary matrix $B_{0}$ to
be the empty matrix.

Hodge Laplacians. The Hodge Laplacians of an simplicial complex $\mathcal{S}$
are important linear maps with close connection to many properties of the
underlying SC. The $k$-th Hodge Laplacian $L_{k}$ is a square matrix indexed
over the $k$-simplices of $\mathcal{S}$.

######  Definition 2.2 (Hodge Laplacian).

For an SC $\mathcal{S}$ with associated boundary matrices
$\left(B_{i}\right)$, we denote by $L_{k}^{\text{up}}\coloneqq
B_{k+1}B_{k+1}^{\top}$ the $k$-th Up-Laplacian, by
$L_{k}^{\text{down}}\coloneqq B_{k}^{\top}B_{k}$ the $k$-th Down-Laplacian,
and finally by their sum $L_{k}\coloneqq
L_{k}^{\text{down}}+L_{k}^{\text{up}}$ the $k$-th _Hodge Laplacian_. We note
that $L_{0}$ recovers the graph Laplacian.

######  Theorem 2.3 (Hodge Decomposition [23, 9, 10, 24]).

For an SC $\mathcal{S}$ with boundary matrices $\left(B_{i}\right)$, Hodge
Laplacians $\left(L_{i}\right)$, and simplicial signal spaces
$\left(C_{i}\right)$, we have that

$\mathbb{R}^{|\mathcal{S}_{k}|}\cong C_{k}=\hbox
to0.0pt{$\displaystyle\overbracket{\phantom{\underbrace{\operatorname{Im}B_{k}^{\top}}_{\text{gradient
space}}\oplus\overbrace{\underbrace{\ker L_{k}}_{\text{harmonic
space}}}}}^{\ker
L_{k}^{\text{up}}}$\hss}\underbrace{\operatorname{Im}B_{k}^{\top}}_{\text{gradient
space}}\oplus\overbracket{\underbrace{\ker L_{k}}_{\text{harmonic
space}}\oplus\underbrace{\operatorname{Im}B_{k+1}}_{\text{curl
space}}}^{\phantom{\ker L_{k}^{\text{down}}}\ker L_{k}^{\text{down}}}.$

For $k=1$, the _harmonic_ space characterises edge flows around the holes of
the SC; the _gradient_ space characterises flows that can be generated by
assigning “potentials” to every vertex and computing their difference along
every edge; the curl space characterises flow “swirling” within the SC, and
any curl flow can be created by assigning a (magnetic) potential to every
2-simplex in the SC and then computing the induced circulations around the
boundaries of the 2-simplices. See Figure 1 for a visualisation of different
flow types. This structure has been exploited in many applications, with the
harmonic space receiving the largest attention ([25] for trajectory outlier
prediction, [15] for simplicial clustering, and [16] for topological point
cloud clustering).

$\alpha$-Complexes and the Delaunay Triangulation. Given a point cloud
$X\subset\mathbb{R}^{2}$, a _Delaunay_ triangulation is a triangulation such
that no point in $X$ is contained inside the circumcircle of any of the
triangles. We can extend this to an SC $\text{DT}(X)$ by taking $X$ as the
vertices, the edges as the $1$-simplices, and the triangles as the
$2$-simplices. This construction naturally generalises to higher dimensions by
dividing the convex hull of the points in $\mathbb{R}^{n}$ into $n$-simplices.
This construction provides a way to construct a sparse high-dimensional
complex on top of data points in $\mathbb{R}^{n}$, on which we can then use
ideas of topological signal processing. However, we sometimes want to factor
in a notion of “scale”, depending on which we might decide to connect only
vertices which are “close”. This leads to the notion of $\alpha$-complexes
[26], where $\alpha$ is essentially the threshold distance for a pair of
points to be considered “close”:

######  Definition 2.4 ($\alpha$-complexes).

Given a point cloud $X\in\mathbb{R}^{n}$ with associated generalised Delaunay
triangulation $\text{DT}(X)$ and a parameter $\alpha$. Then, the
$\alpha$-complex is the SC $\mathcal{S}_{\alpha}(X)$ with vertices $X$ and
simplices $S_{\alpha}(X)=\left\\{\sigma\in\text{DT}(X):|x-y|\leq\alpha\forall
x,y\in\sigma\right\\}.$

For $\alpha\leq\alpha^{\prime}$, we have the canonical inclusion
$\mathcal{S}_{\alpha}(X)\subset\mathcal{S}_{\alpha^{\prime}}(X)$. Because of
their origin in the Delaunay triangulation, the $\alpha$-complexes contain
only very few edges in comparison to the popular Vietoris–Rips complexes,
while retaining the same homotopy type. This results in good computational
performance in downstream tasks.

Persistent Homology and the Persistent Laplacians. In many applications, we
have no a-priori knowledge of how to pick the scale-parameter $\alpha$
correctly. Varying $\alpha$ in the $\alpha$-complexes of $X$ enables us to
extract features of $X$ across all possible scales. Different methods of
defining, analysing, and subsuming these features have been proposed. One of
the most popular tools in this context is Persistent Homology, a tool from
algebraic topology that tracks the evolution of the “shape” of a point cloud
across the scales [27], in terms of the number of connected components, holes,
cavities, etc. However, some interesting properties are not as black-and-white
as the existence of holes: For graphs, the second-smallest eigenvalue
$\lambda_{2}$ of the graph Laplacian is a measure of how easily the graph can
be divided into separate clusters [4, 2, 3]. The extreme value of $0$ means
that the graph is disconnected, while the maximum of $n$ is only achieved in
complete graphs; all other graphs range somewhere in between. Hence these
“continuous” spectral quantities can contain more information than their
binary counterparts from algebraic topology. This has motivated the study of
the persistent Laplacians [11], which are in a special case the Hodge
Laplacians $L(\mathcal{S}_{\alpha}(X))$ of the $\alpha$-complexes.

##  3 Eigenvectors of the Hodge Laplacian

Fig. 2: Evolution of the smallest eigenvalues of Hodge Laplacian across the
steps of $\alpha$-filtration of the point cloud of Figure 1. _Top:_ We witness
two distinct phases in eigenvalue behaviour. In the first regime of low
filtration values, the associated $\alpha$-complexes are sparsely connected
and exhibit many holes (large number of number of harmonic eigenvalues
$\lambda=0$). There is a V-shaped pattern of gradient eigenvalues. This
correspond to the fact that the clusters first grow, which creates very “thin”
bridges (small gradient eigenvalues) between growing clusters. Once the
clusters start to grow together the connecting bridges become “wider” (more
connecting edges), and the gradient eigenvalues increase again. In the second
regime of medium to high filtration values, the associated $\alpha$-complexes
grow denser and the holes are filled. Only one harmonic eigenvalue remains,
whereas three structural gradient eigenvalues remain relatively small.
However, the most distinct feature of the second regime is the cohort of
decreasing curl eigenvectors/-values overtaking the gradient vectors of the
initial V-shape. The smallest four of these eigenvalues are separated from the
rest by a significant margin. They represent the curl flow generated around
each of the four clusters. _Bottom:_ Number of eigenvalues in smallest $40$
being curl, gradient, and harmonic.

Harmonic, gradient, and curl eigenvectors. In the previous section, we have
introduced the concept of three different types of eigenvectors of Hodge
Laplacians: harmonic, curl, and gradient eigenvectors. In Figure 1 we have
depicted the eigenvectors corresponding to the smallest $8$ eigenvalues at
three steps of an $\alpha$-filtration on a point cloud sampled from four disks
arranged in a square. It is immediately clear that a small eigenvalue can
encode very different properties for the curl or the gradient case: A small
gradient eigenvalue arises from two clusters of the data connected by a small
“bridge”. A small curl eigenvalue arises from a densely connected cluster with
large radius. Finally, a harmonic eigenvalue just means that there exists a
hole in the data set, without any information on its size.

Tracking the eigenvectors across multiple filtration steps. Figure 1 reveals
that there is a strong connection between eigenvectors in different steps of
the filtration. However, these eigenvectors equivalence classes frequently
swap order. Thus, simply treating the _$n$ -th smallest eigenvalue/vector
pair_ as a fixed instance does not respect the underlying structure of the
filtration. Accordingly, we propose to connect the eigenvectors based on
similarity. However, defining such a similarity is not straight-forward
because the Hodge Laplacians of different filtration steps operate on
different simplicial signal spaces. Furthermore, being an eigenvector is
invariant under multiplication by a scalar and hence two similar eigenvectors
can differ by scalar multiplication by $-1$.

######  Definition 3.1 (Persistent eigenvector similarity).

Let $X\in\mathbb{R}^{n}$ be a point cloud, $k\in\mathbb{Z}_{\geq 0}$ a degree,
$0<\alpha<\alpha^{\prime}$ two filtration values, and $L_{\alpha}$ and
$L_{\alpha^{\prime}}$ the $k$-th Hodge Laplacians of the $\alpha$-complexes
$\mathcal{S}_{\alpha}(X)$ and $\mathcal{S}_{\alpha^{\prime}}(X)$. Furthermore,
let $v$ and $v^{\prime}$ be two eigenvectors of $L_{\alpha}$ and
$L_{\alpha^{\prime}}$ respectively. Finally, let
$\iota\colon\mathbb{R}^{|\mathcal{S}_{\alpha,k}(X)|}\hookrightarrow\mathbb{R}^{|\mathcal{S}_{\alpha^{\prime},k}(X)|}$
be the inclusion of the simplicial signal spaces induced by the inclusion of
$\alpha$-complexes
$\mathcal{S}_{\alpha}\hookrightarrow\mathcal{S}_{\alpha^{\prime}}$. We then
define the persistent eigenvector similarity (PES) of $v$ and $v^{\prime}$ by
$\text{PES}(v,v^{\prime})\coloneqq|\iota(v)^{\top}v^{\prime}|/\left(\lVert
v\rVert_{2}\cdot\lVert v^{\prime}\rVert_{2}\right).$

######  Definition 3.2 (Persistent eigenvector matching (PEM)).

Let $X\in\mathbb{R}^{n}$ be a point cloud, $k\in\mathbb{Z}_{\geq 0}$ a degree,
$0<\alpha<\alpha^{\prime}$ two filtration values, and $L_{\alpha}$ and
$L_{\alpha^{\prime}}$ the $k$-th Hodge Laplacians of the $\alpha$-complexes
$\mathcal{S}_{\alpha}(X)$ and $\mathcal{S}_{\alpha^{\prime}}(X)$. Furthermore,
let $V$ and $V^{\prime}$ be two sets of eigenvectors of $L_{\alpha}$ and
$L_{\alpha^{\prime}}$. We then define the PEM:

$\displaystyle\text{PEM}_{X,\alpha,\alpha^{\prime}}$
$\displaystyle(V,V^{\prime})\coloneqq\smash{\left\\{(v,v^{\prime})\in V\times
V^{\prime}:\phantom{\min_{\hat{v}\in
V}\text{PES}(\hat{v},v^{\prime})=\text{PES}(v,v^{\prime})=\min_{\hat{v}^{\prime}\in
V^{\prime}}\text{PES}(v,\hat{v}^{\prime})}\right.}$
$\displaystyle\left.\min_{\hat{v}\in
V}\text{PES}(\hat{v},v^{\prime})=\text{PES}(v,v^{\prime})=\min_{\hat{v}^{\prime}\in
V^{\prime}}\text{PES}(v,\hat{v}^{\prime})\right\\}$

For non-degenerate cases, i.e. when no two eigenvector pairs have the same
PES, persistent eigenvector matching provides a one-to-one matching (possibly
not matching some of the eigenvectors). We can then use this matching to track
the evolution of different eigenvalues across the stages of the filtration,
see Figure 2.

##  4 Applications

Algorithm 1 Topological Spectral Clustering (TSC)

Input: Simplicial (_/cellular/cubical_) complex $\mathcal{S}$, dimension $n$,
#eigenvectors $h$, #clusters $c$, selection of $C\in\\{\text{gradient,
harmonic, curl, total}\\}$ eigenmodes.

1. Compute smallest $h$ $C$-eigenvectors $v_{1},\dots,v_{h}$ of $L_{n}(\mathcal{S})$.

2. Set $v(\sigma_{n})\coloneqq\left(v_{1}\left[\sigma_{n}\right],\dots,v_{h}\left[\sigma_{n}\right]\right)$.

3. Set $\hat{v}(\sigma_{n})\coloneqq\text{sgn}\left(1^{\top}v(\sigma_{n})\right)v(\sigma_{n})$.

4. Cluster $\mathcal{S}_{n}$ using $k$-means on $\\{\hat{v}(\sigma_{n}):\sigma_{n}\in\mathcal{S}_{n}\\}\subset\mathbb{R}^{|\mathcal{S}_{n}|}$.

 _Optional:_ Cluster points in $\mathcal{S}_{0}$ wrt. adjacent $n$-simplices.

Output: Cluster-ids of $k$-simplices $\mathcal{S}_{n}$ of step 4.

(a) Gradient Clustering.

(b) Curl Clustering.

(c) Harmonic Clustering.

(d) Total Clustering.

Fig. 3: Topological Spectral Clustering. In contrast to ordinary
$0$-dimensional spectral clustering, higher-order spectral clustering allows
for the subcategories of _gradient_ , _curl_ and _harmonic_ clustering,
highlighting different features of the simplicial complex.

In the previous sections, we have argued why providing additional information
on the often considered eigenvalues of Hodge Laplacians is crucial. We now
give a range of possible applications building on top of these considerations.

Topological spectral clustering. We generalise the traditional notion of
spectral clustering to higher-order simplicial complexes, see Algorithm 1.
Given a simplicial complex $\mathcal{S}$ with $n$-th Hodge Laplacian $L$ and
eigenvectors $v_{1},\dots,v_{k}$, we will cluster an $n$-simplex
$\sigma_{n}\in\mathcal{S}_{n}$ based on its associated eigenvector entries
$v(\sigma_{n})\coloneqq\left(v_{1}\left[\sigma_{n}\right],\dots,v_{k}\left[\sigma_{n}\right]\right)$.
However, in the construction of the Hodge Laplacian we have to fix an
arbitrary orientation of the simplices, and changing this arbitrary
orientation leads to sign changes in the eigenvector entries of $\sigma_{n}$.
Hence we need to cluster the set
$\\{v(\sigma_{n}):\sigma_{n}\in\mathcal{S}_{n}\\}/(v\sim-v)$ inside
$\mathbb{R}^{|\mathcal{S}_{n}|}$ where we identify $v$ and $-v$ for all
$v\in\mathbb{R}^{|\mathcal{S}_{n}|}$. However, clustering algorithms work best
on plain real vector spaces. Thus we use $k$-means on
$\\{\text{sgn}\left(1^{\top}v(\sigma_{n})\right)v(\sigma_{n}):\sigma_{n}\in\mathcal{S}_{n}\\}$.
Note that we have made the embedding of $\sigma_{n}$ invariant under
orientation change. Building on our previous considerations, we can choose a
subset of the curl, gradient, and harmonic eigenvectors to highlight different
properties of our point cloud and our simplicial complex, see Figure 3. By
aggregating the cluster-ids of the incident edges for each node, we can obtain
an induced clustering of the points similar to [16].

Fig. 4: Edge roles using HGC-values. Relative strength of _blue:_ harmonic,
_red:_ curl, _green:_ gradient contributions. _Top:_ $\alpha$-complexes on two
3d point clouds sampled on real world objects of ModelNet40 data set [28].
_Bottom:_ $\alpha$-complexes on two synthetic point clouds.

Edge roles using HGC-values. Edges ($n$-simplices) can have many different
roles inside an SC: Some edges connect nodes in the centre of a cluster, some
are located in the outer rim of a cluster, or connect different clusters, and
again others are located on the edge of “empty spaces” in the SC, or have no
relationship to any of the clusters. This role based clustering is orthogonal
to TSC. Edge role classification plays an important role in applications on
brain networks in neuroscience [29], link analysis in social networks [30],
etc. We provide a metric to extract these relevant information:

######  Definition 4.1 (HGC-values).

Given a simplicial complex $\mathcal{S}$, a dimension $n$ and a number $k$,
let $V\coloneqq v_{1},\dots,v_{k}$ be the eigenvectors associated to the $k$
smallest eigenvalues of the $n$-th Hodge Laplacian $L_{n}$ of $\mathcal{S}$.
Furthermore, we denote by $V_{\text{harm}}\coloneqq\\{v_{i}\in
V:L_{n}v_{i}=0\\}$, $V_{\text{curl}}\coloneqq\\{v_{i}\in
V:L^{\text{up}}_{n}v_{i}\not=0\\}$, $V_{\text{grad}}\coloneqq\\{v_{i}\in
V:L^{\text{down}}_{n}v_{i}\not=0\\}$ the harmonic, curl, and gradient
eigenvectors. Let $e_{\text{max}}$ denote the entry with maximum absolute
value across all $v\in V$. Setting the max over $\emptyset$ to be $0$, we can
define the HGC-value of the simplex $\sigma_{n}$ $\text{HGC}(\sigma_{n})$ to
be

$\left(\max_{v\in V_{\text{harm}}}(|v(\sigma_{n})|),\max_{v\in
V_{\text{grad}}}(|v(\sigma_{n})|),\max_{v\in
V_{\text{curl}}}(|v(\sigma_{n})|)\right)/|e_{\text{max}}|.$

Intuitively, the HGC-values measure the relevance of an edge for harmonic,
gradient, and curl eigenvectors. Combinations of these then correspond to the
edge roles discussed above, see Figure 4 for an illustration. Although
normalising by the largest entry of any eigenvector worked best in practice,
other normalisations are possible as well. For a different approach to edge
roles based on harmonic, curl, and gradient vectors see [20].

Discussion. In this paper, we have argued why disentangling the roles and
properties of different small eigenvalues of the Hodge Laplacian is beneficial
in theory and applications. Building on this insight, we have introduced a
method to track individual eigenvalues through an $\alpha$-filtration,
introduced a novel tunable spectral clustering algorithm on higher-order
simplicial complexes, and introduced the notion of HGC-values to extract
simplicial roles. Finally, we showed the expressiveness of topological
spectral clustering and the HGC-values in experiments on synthetic and real-
world data.

## References

  * [1] Chris Godsil and Gordon F Royle,  Algebraic graph theory, vol. 207,  Springer Science & Business Media, 2001.
  * [2] Fan RK Chung,  Spectral graph theory, vol. 92,  American Mathematical Soc., 1997.
  * [3] Daniel Spielman,  “Spectral graph theory,”  Combinatorial scientific computing, vol. 18, pp. 18, 2012.
  * [4] Ulrike von Luxburg,  “A tutorial on spectral clustering,”  Statistics and Computing, vol. 17, no. 4, pp. 395–416, Dec 2007\.
  * [5] Marcelo Fiori and Guillermo Sapiro,  “On spectral properties for graph matching and graph isomorphism problems,”  Information and Inference: A Journal of the IMA, vol. 4, no. 1, pp. 63–76, 2015.
  * [6] Sudip Saha, Abhijin Adiga, B Aditya Prakash, and Anil Kumar S Vullikanti,  “Approximation algorithms for reducing the spectral radius to control epidemic spread,”  in Proceedings of the 2015 SIAM International Conference on Data Mining. SIAM, 2015, pp. 568–576.
  * [7] Tuncer Can Aysal, Boris N Oreshkin, and Mark J Coates,  “Accelerated distributed average consensus via localized node state prediction,”  IEEE Transactions on signal processing, vol. 57, no. 4, pp. 1563–1576, 2008.
  * [8] Stephen Boyd, Persi Diaconis, and Lin Xiao,  “Fastest mixing markov chain on a graph,”  SIAM review, vol. 46, no. 4, pp. 667–689, 2004.
  * [9] Sergio Barbarossa and Stefania Sardellitti,  “Topological signal processing over simplicial complexes,”  IEEE Transactions on Signal Processing, vol. 68, pp. 2992–3007, 2020.
  * [10] Michael T. Schaub, Yu Zhu, Jean-Baptiste Seby, T. Mitchell Roddenberry, and Santiago Segarra,  “Signal processing on higher-order networks: Livin’ on the edge… and beyond,”  Signal Processing, vol. 187, pp. 108149, 2021.
  * [11] Rui Wang, Duc Duy Nguyen, and Guo-Wei Wei,  “Persistent spectral graph,”  International journal for numerical methods in biomedical engineering, vol. 36, no. 9, pp. e3376, 2020.
  * [12] Facundo Mémoli, Zhengchao Wan, and Yusu Wang,  “Persistent laplacians: Properties, algorithms and implications,”  SIAM Journal on Mathematics of Data Science, vol. 4, no. 2, pp. 858–884, 2022.
  * [13] Thomas Davies, Zhengchao Wan, and Ruben J Sanchez-Garcia,  “The persistent Laplacian for data science: Evaluating higher-order persistent spectral representations of data,”  in Proceedings of the 40th International Conference on Machine Learning, Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, Eds. 23–29 Jul 2023, vol. 202 of Proceedings of Machine Learning Research, pp. 7249–7263, PMLR.
  * [14] Jian Liu, Jingyan Li, and Jie Wu,  “The algebraic stability for persistent laplacians,”  arXiv preprint arXiv:2302.03902, 2023.
  * [15] Stefania Ebli and Gard Spreemann,  “A notion of harmonic clustering in simplicial complexes,”  in 2019 18th IEEE International Conference On Machine Learning And Applications (ICMLA), 2019, pp. 1083–1090.
  * [16] Vincent P Grande and Michael T Schaub,  “Topological point cloud clustering,”  in Proceedings of the 40th International Conference on Machine Learning, Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, Eds. 23–29 Jul 2023, vol. 202 of Proceedings of Machine Learning Research, pp. 11683–11697, PMLR.
  * [17] Sanjukta Krishnagopal and Ginestra Bianconi,  “Spectral detection of simplicial communities via hodge laplacians,”  Physical Review E, vol. 104, no. 6, pp. 064303, Dec 2021.
  * [18] Thummaluru Siddartha Reddy, Sundeep Prabhakar Chepuri, and Pierre Borgnat,  “Clustering with simplicial complexes,”  arXiv preprint arXiv:2303.07646, 2023.
  * [19] Kiya W Govek, Venkata S Yamajala, and Pablo G Camara,  “Clustering-independent analysis of genomic data using spectral simplicial theory,”  PLoS computational biology, vol. 15, no. 11, pp. e1007509, 2019\.
  * [20] Michael T Schaub, Jörg Lehmann, Sophia N Yaliraki, and Mauricio Barahona,  “Structure of complex networks: Quantifying edge-to-edge relations by failure-induced flow redistribution,”  Network Science, vol. 2, no. 1, pp. 66–89, 2014.
  * [21] Allen Hatcher,  Algebraic Topology,  Cambridge University Press, Cambridge, 2002.
  * [22] G.E. Bredon, J.H. Ewing, F.W. Gehring, and P.R. Halmos,  Topology and Geometry,  Graduate Texts in Mathematics. Springer, New York, 1993.
  * [23] Lek-Heng Lim,  “Hodge Laplacians on graphs,”  vol. 62, no. 3, pp. 685–715, 2020.
  * [24] T Mitchell Roddenberry, Nicholas Glaze, and Santiago Segarra,  “Principled simplicial neural networks for trajectory prediction,”  in International Conference on Machine Learning. PMLR, 2021, pp. 9020–9029.
  * [25] Florian Frantzen, Jean-Baptiste Seby, and Michael T Schaub,  “Outlier detection for trajectories via flow-embeddings,”  in 2021 55th Asilomar Conference on Signals, Systems, and Computers. IEEE, 2021, pp. 1568–1572.
  * [26] Herbert Edelsbrunner,  “Alpha shapes-a survey,”  in Tessellations in the Sciences: Virtues, Techniques and Applications of Geometric Tilings. 2011.
  * [27] Herbert Edelsbrunner and John Harer,  “Persistent homology—a survey,”  Contemporary mathematics, vol. 453, no. 26, pp. 257–282, 2008.
  * [28] Zhirong Wu, Shuran Song, Aditya Khosla, Fisher Yu, Linguang Zhang, Xiaoou Tang, and Jianxiong Xiao,  “3d shapenets: A deep representation for volumetric shapes,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 1912–1920.
  * [29] Joshua Faskowitz, Richard F. Betzel, and Olaf Sporns,  “Edges in brain networks: Contributions to models of structure and function,”  Network Neuroscience, vol. 6, no. 1, pp. 1–28, 02 2022.
  * [30] Pasquale De Meo, Emilio Ferrara, Giacomo Fiumara, and Angela Ricciardello,  “A novel measure of edge centrality in social networks,”  Knowledge-Based Systems, vol. 30, pp. 136–150, 2012.