The Derivative of Kemeny's Constant as a Centrality Measure in Undirected Graphs

TL;DR

This paper introduces new centrality measures based on the directional derivatives of Kemeny's constant in undirected graphs, providing tools for analyzing connectivity, sensitivity, and link prediction in network structures.

Contribution

It defines novel edge and non-edge centrality measures using derivatives of Kemeny's constant, with explicit formulas and algorithms, extending previous work and applicable to various network analyses.

Findings

Derived explicit formulas involving the inverse Laplacian for centrality measures.

Connected the new measures to existing ones in prior literature.

Applied the measures to road networks and link prediction tasks.

Abstract

Kemeny's constant quantifies a graph's connectivity by measuring the average time for a random walker to reach any other vertex. We introduce two concepts of the directional derivative of Kemeny's constant with respect to an edge and use them to define centrality measures for edges and non-edges in the graph. Additionally, we present a sensitivity measure of Kemeny's constant. An explicit expression for these quantities involving the inverse of the modified graph Laplacian is provided, which is valid even for cut-edges. These measures are connected to the one introduced in [Altafini et al., SIMAX 2023], and algorithms for their computation are included. The benefits of these measures are discussed, along with applications to road networks and link prediction analysis. For one-path graphs, an explicit expression for these measures is given in terms of the edge weights.

Tables

Table 1. Table 1: The correlation coefficients between pairs of the five measures.

	J	A	R	C	K
Jaccard	1	0.83	0.76	0.56	-0.15
Adamic–Adar	0.83	1	0.96	0.64	-0.17
Resource allocation	0.76	0.96	1	0.56	-0.15
Common neighbour	0.56	0.64	0.56	1	-0.59
Kemeny derivative	-0.15	-0.17	-0.15	-0.59	1

Table 2. Table 2: The node labels for the top ten new links predicted by the new measure.

1	RAVASZ, E	YOOK, S
2	NEDA, Z	YOOK, S
3	GOMEZGARDENES, J	PACHECO, A
4	SCHUBERT, A	YOOK, S
5	FARKAS, I	YOOK, S
6	FERRERICANCHO, R	MONTOYA, J
7	VALVERDE, S	MONTOYA, J
8	RAVASZ, E	TOMBOR, B
9	VESPIGNANI, A	WEIGT, M
10	BARRAT, A	BOGUNA, M

Equations

$$-1\le {\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_{n-1}<{\lambda }_n=1,$$

$$\kappa (G)=trace(I-P+1{\mathbf{v}}^{\top }{)}^{-1}-1,\mathbf{v}\in {\mathbb{R}}^n,{\mathbf{v}}^{\top }1=1,$$

$$I-P+\frac{1}{\Vert \mathbf{d}{\Vert }_1}1{\mathbf{d}}^{\top }=D^{-\frac{1}{2}}\left(I-D^{-\frac{1}{2}}A{D}^{-\frac{1}{2}}+\frac{1}{\Vert \mathbf{d}{\Vert }_1}{\mathbf{d}}^{\frac{1}{2}}{\mathbf{d}}^{\frac{1}{2}\top }\right)D^{\frac{1}{2}},$$

$$\begin{array}{rl}& H=I-D^{-\frac{1}{2}}A{D}^{-\frac{1}{2}}+\frac{1}{\Vert \mathbf{d}{\Vert }_1}{\mathbf{d}}^{\frac{1}{2}}{\mathbf{d}}^{\frac{1}{2}\top }=D^{-\frac{1}{2}}S{D}^{-\frac{1}{2}},\\ & S=L+\frac{1}{\Vert \mathbf{d}{\Vert }_1}\mathbf{d}{\mathbf{d}}^{\top },L=D-A.\end{array}$$

$$\kappa (G)=trace(H^{-1})-1=trace(D^{\frac{1}{2}}{S}^{-1}{D}^{\frac{1}{2}})-1.$$

$$\kappa (G)=\sum _{i=1}^{n-1}\frac{1}{1-{\lambda }_i}.$$

$$\begin{array}{rl}c(p,q)& =\kappa (P_{pq})-\kappa (P)=trace(D^{\frac{1}{2}}(S_{pq}^{-1}-S^{-1})D^{\frac{1}{2}}),\\ S_{pq}& =S-a_{pq}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top },\end{array}$$

$$c(p,q)=a_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1}D{S}_{pq}^{-1}{\mathbf{w}}_{pq}.$$

$$c_r(p,q)=a_{pq}{\mathbf{w}}_{pq}^{\top }(rD+S{)}^{-1}D(rD+S_{pq}{)}^{-1}{\mathbf{w}}_{pq}.$$

$$c^F(p,q)=\underset{r\to 0}{\lim }\left(\frac{1}{r}-c_r(p,q)\right),$$

$$A_{pq}(t)=A+t{a}_{pq}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top },{\mathbf{w}}_{pq}={\mathbf{e}}_p-{\mathbf{e}}_q.$$

$$\begin{array}{rl}& H_{pq}(t)=H-t{a}_{pq}{D}^{-\frac{1}{2}}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top }{D}^{-\frac{1}{2}},\\ & S_{pq}(t)=S-t{a}_{pq}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top },\\ & L_{pq}(t)=L-t{a}_{pq}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top },\\ & P_{pq}(t)=P+t{a}_{pq}{D}^{-1}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top }.\end{array}$$

$$\frac{\kappa (P_{pq}(t))-\kappa (P)}{t}$$

$$\mu (p,q)=\underset{t\to 0^{+}}{\lim }\frac{\kappa (P_{pq}(t))-\kappa (P)}{t}.$$

$$\bar{\mu }(p,q)=\underset{t\to 0^{+}}{\lim }\frac{\kappa (P+t{D}^{-1}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top })-\kappa (P)}{t}.$$

$$\mu (p,q)=a_{pq}\bar{\mu }(p,q).$$

$$\alpha (p,q)={\mathbf{w}}_{pq}^{\top }{S}^{-1}{\mathbf{w}}_{pq},\beta (p,q)={\mathbf{w}}_{pq}^{\top }(S^{-1}D{S}^{-1}){\mathbf{w}}_{pq};$$

$$\begin{array}{rl}& \kappa (P_{pq}(t))=\kappa (P)+\frac{t{a}_{pq}}{1-t{a}_{pq}\alpha (p,q)}\beta (p,q),\\ & {\kappa }^{(i)}(P_{pq}(t))=\frac{i!(a_{pq}{)}^i}{(1-t{a}_{pq}\alpha (p,q){)}^{i+1}}\alpha (p,q{)}^{i-1}\beta (p,q),\end{array}$$

$$S_{pq}(t{)}^{-1}=S^{-1}+\frac{t{a}_{pq}}{1-t{a}_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1}{\mathbf{w}}_{pq}}{S}^{-1}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1},$$

$$\begin{array}{rl}\kappa (P_{pq}(t))& =trace(D^{\frac{1}{2}}{S}^{-1}{D}^{\frac{1}{2}})\\ & +\frac{t{a}_{pq}}{1-t{a}_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1}{\mathbf{w}}_{pq}}trace(D^{\frac{1}{2}}{S}^{-1}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1}{D}^{\frac{1}{2}})-1\\ & =\kappa (P)+\frac{t{a}_{pq}}{1-t{a}_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1}{\mathbf{w}}_{pq}}{\mathbf{w}}_{pq}^{\top }(S^{-1}D{S}^{-1}){\mathbf{w}}_{pq},\end{array}$$

$$\begin{array}[]{ll}\mu(p,q)=a_{pq}\beta(p,q),&\mu^{(i)}(p,q)=i!a_{pq}^{i}\alpha(p,q)^{i-1}\beta(p,q),\\[6.45831pt] \bar{\mu}(p,q)=\beta(p,q),&\bar{\mu}^{(i)}(p,q)=i!a_{pq}^{i-1}\alpha(p,q)^{i-1}\beta(p,q),\end{array}$$

$$\begin{array}{rl}& 2\underset{i}{\min }{\zeta }_i\le \bar{\mu }(p,q)\le 2\underset{i}{\max }{\zeta }_i,\\ & 2\underset{i}{\min }{d}_i\underset{i}{\min }\left({\nu }_i^{-2}\right)\le \bar{\mu }(p,q)\le 2\underset{i}{\max }{d}_i\underset{i}{\max }\left({\nu }_i^{-2}\right),\\ & 2\underset{i}{\min }\left(d_i^{-1}\right)\min \left\{\frac{1}{(1-{\lambda }_1{)}^2},1\right\}\le \bar{\mu }(p,q)\le 2\underset{i}{\max }\left(d_i^{-1}\right)\max \left\{\frac{1}{(1-{\lambda }_{n-1}{)}^2},1\right\}\end{array}$$

$$\bar{\mu }(p,q)={\mathbf{w}}_{pq}^{\top }(S^{-1}D{S}^{-1}){\mathbf{w}}_{pq}=2\frac{{\mathbf{w}}_{pq}^{\top }(S^{-1}D{S}^{-1}){\mathbf{w}}_{pq}}{{\mathbf{w}}_{pq}^{\top }{\mathbf{w}}_{pq}}.$$

$$\bar{\mu }(p,q)=2\frac{{\mathbf{y}}^{\top }D\mathbf{y}}{{\mathbf{y}}^{\top }\mathbf{y}}\cdot \frac{{\mathbf{w}}_{pq}^{\top }{S}^{-2}{\mathbf{w}}_{pq}}{{\mathbf{w}}_{pq}^{\top }{\mathbf{w}}_{pq}}.$$

$$\bar{\mu }(p,q)=2\frac{{\mathbf{z}}^{\top }{H}^{-2}\mathbf{z}}{{\mathbf{z}}^{\top }\mathbf{z}}\cdot \frac{{\mathbf{w}}_{pq}^{\top }{D}^{-1}{\mathbf{w}}_{pq}}{{\mathbf{w}}_{pq}^{\top }{\mathbf{w}}_{pq}},\mathbf{z}=D^{-\frac{1}{2}}{\mathbf{w}}_{pq}.$$

$$c(p,q)=\frac{\mu (p,q)}{1-a_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1}{\mathbf{w}}_{pq}}.$$

$$S_{pq}^{-1}=S^{-1}\left(I+\frac{a_{pq}}{1-a_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1}{\mathbf{w}}_{pq}}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top }{S}^{-1}\right).$$

$${\mu }_r(p,q)=a_{pq}{\mathbf{w}}_{pq}^{\top }{S}_r^{-1}D{S}_r^{-1}{\mathbf{w}}_{pq},S_r=S+rD,$$

$$(rD+S_{pq}{)}^{-1}=S_r^{-1}+\frac{a_{pq}}{1-a_{pq}{\mathbf{w}}_{pq}^{\top }{S}_r^{-1}{\mathbf{w}}_{pq}}{S}_r^{-1}{\mathbf{w}}_{pq}{\mathbf{w}}_{pq}^{\top }{S}_r^{-1},$$

$$\begin{array}{rl}{c}_r(p,q)=& {\mu }_r(p,q)\left(1+\frac{a_{pq}}{1-a_{pq}{\mathbf{w}}_{pq}^{\top }{S}_r^{-1}{\mathbf{w}}_{pq}}\cdot {\mathbf{w}}_{pq}^{\top }{S}_r^{-1}{\mathbf{w}}_{pq}\right)\\ =& \frac{1}{1-a_{pq}{\mathbf{w}}_{pq}^{\top }{S}_r^{-1}{\mathbf{w}}_{pq}}{\mu }_r(p,q).\end{array}$$

Full text

Keywords: Kemeny Constant; Graph; Edge Centrality; Network Centrality; Link Prediction.

1 Introduction

The significance of edges or nodes in a connected graph $G$ is generally assessed by computing a function of the adjacency matrix [4, 19] or by modeling a random walk on the graph using a Markov chain and analyzing its performance metrics. One such metric is Kemeny’s constant $\kappa(G)$, which represents the expected number of steps for the chain to travel in the graph from a starting state $i$ to a randomly selected destination state, chosen according to the chain’s stationary distribution. This expected time remains the same regardless of the initial state $i$, hence the name ”constant” [21]. Algebraically, Kemeny’s constant is given by the trace of a pseudoinverse of the matrix $I-P$, where $P$ is the transition probability matrix of the Markov chain [20].

In the context of graphs, Kemeny’s constant is used, for example, to evaluate the efficiency of robotic surveillance strategies [26], to partition nodes into clusters [7], to control road network dynamics [13], and more. The effects on Kemeny’s constant of the removal or addition of edges in tree graphs are studied in [23], while a structured condition number for Kemeny’s constant is introduced in [11].

In [3], the variation of Kemeny’s constant after removing an edge is used to measure the importance of that edge. More specifically, given an undirected connected weighted graph $G=(V,E)$, where $V=\{1,\ldots,n\}$ is the set of vertices and $E\subset V\times V$ is the set of edges, the measure $c(e)=\kappa(\widehat{G})-\kappa(G)$ is introduced, where $e=\{p,q\}\in E$, $\widehat{G}$ is the graph obtained from $G$ by removing the edge $\{p,q\}$ and adding two loops $\{p,p\}$ and $\{q,q\}$. The addition of the two loops is introduced to avoid the Braess paradox [3]. The definition of this measure does not apply if $e$ is a cut-edge, that is, an edge whose removal disconnects the graph, since, in this case, $\widehat{G}$ is not connected and $\kappa(\widehat{G})$ is infinite. To overcome this issue, the authors of [3] introduced a regularized and filtered version of the measure $c(e)$; this version is easy to compute, but, for cut-edges, it is prone to numerical cancellation issues and requires cut-edge detection. The former difficulty is addressed in [9], where an equivalent formula for $c(e)$ is provided for cut-edges only. However, this workaround requires the previous detection of cut-edges in the graph and a different treatment for normal edges and cut-edges.

In this paper, instead of considering the variation of Kemeny’s constant after an edge removal and the addition of suitable loops, we analyze its infinitesimal variation by introducing a suitable directional derivative of Kemeny’s constant. A definition of the derivative of Kemeny’s constant has already been introduced in [7] for cluster separation in graphs, and in [2] for estimating the effect of the removal of edges in particular graphs. Here, we introduce two definitions of directional derivative of Kemeny’s constant with respect to an edge of the graph: a weighted and an unweighted derivative.

The weighted derivative with respect to the edge $e=\{p,q\}\in E$, denoted by $\mu(p,q)$, is shown to be well defined and always positive. It quantifies how quickly Kemeny’s constant increases as the weight of edge $e$ decreases. Roughly speaking, a high value of $\mu(p,q)$ indicates that a small reduction in the edge’s weight leads to a significant rise in Kemeny’s constant, making it a useful measure of the edge’s importance. Unlike for the measure introduced in [3], $\mu(p,q)$ has several key advantages: It is universally defined, i.e., it is well-defined for all edges, including cut-edges, so there is no need to identify them or apply any special rules as filtering and regularization; it is numerically stable since its computation avoids cancellation making it reliable; it is efficient, i.e., it can be computed with the same arithmetic cost as the measure from [3].

The unweighted derivative with respect to the edge $e=\{p,q\}$, denoted by $\bar{\mu}(p,q)$, is shown to be well defined and always positive also for non-edges, i.e., if $e\not\in E$. Moreover, we show that $\mu(p,q)=a_{pq}\bar{\mu}(p,q)$, where $a_{pq}$ is the weight of edge $\{p,q\}$, and the importance of non-edges is also defined through $\bar{\mu}(p,q)$. This latter property finds applications in link prediction in network analysis [28], where the goal is to predict the creation of new edges in a given graph.

Finally, from the definition of the unweighted derivative, we derive a global measure of the sensitivity of Kemeny’s constant to small perturbations of the edge weights. This measure can also be related to the analysis of graph robustness, investigated in [5], and to the graph perturbation needed to enforce a specific centrality measure [12].

We introduce effective algorithms for the computation of $\mu(p,q)$ and $\bar{\mu}(p,q)$, based on the Cholesky factorization of the graph Laplacian matrix. Moreover, we provide explicit expressions of $\mu(p,q)$ for the class of one-path graphs in terms of the weights of the edges.

Numerical experiments conducted on the Pisa roadmap demonstrate that the measure $\mu(p,q)$ effectively identifies roads with high centrality values in practice. The measure $\bar{\mu}(p,q)$ is applied to the link prediction problem in star, one-path, and cyclic graphs, and in a binary tree, to identify new links that would increase the graph’s circulability; moreover, it is applied in a collaboration network to predict new collaborations. The global measure of the sensitivity of Kemeny’s constant is applied to star, one-path, and circulant graphs. The numerical experiments show a behaviour which is expected from the topology of the tested graphs.

The paper is organized as follows. After recalling in Section 2 some preliminary results, we introduce in Section 3 the new measures $\mu(p,q)$ and $\bar{\mu}(p,q)$, prove their positivity, relate them to a suitable rank-1 perturbation of the Laplacian of the graph, and give explicit expressions in terms of the centralities introduced in [3]. Moreover, in Section 3.1 we comment on the condition number of Kemeny’s constant, and in Section 3.2 we introduce a global measure of sensitivity. In Section 4, we analyze some computational issues of these measures and provide an algorithm for their computation. In Section 5, an explicit expression for $\mu(p,q)$ is given for one-path graphs. In Section 6 we present some numerical experiments that show the better qualities of the new measure with respect to those of [3], together with some applications to the link prediction problem and to the global measure of the sensitivity of Kemeny’s constant. In Section 7, we draw some conclusions.

2 Notation and preliminaries

Here and hereafter, we use the following notation: $\bm{1}=(1,1,\ldots,1)^{\top}$ and $\bm{e}_{i}$ is the $i$th column of the identity matrix $I$. A real matrix $S$ is said to be stochastic if it has nonnegative entries and $S\bm{1}=\bm{1}$. Any stochastic matrix $S$ has spectral radius 1, 1 is an eigenvalue, and if $S$ is irreducible then the eigenvalue 1 is simple, by the Perron-Frobenius theorem.

Given an undirected connected graph $G=(V,E)$, where $V=\{1,\ldots,n\}$ is the set of vertices and $E\subset V\times V$ is the set of edges, denote $A=(a_{ij})$ the associated (symmetric) adjacency matrix so that $a_{ij}=a_{ji}>0$ is the weight of the edge that connects $i$ and $j$ for $\{i,j\}\in E$, and $a_{ij}=0$ if $\{i,j\}\not\in E$. We denote the edges with unordered pairs $\{i,j\}$. The size of $A$ is $n=|V|$, where $|V|$ denotes the cardinality of $V$. In the sequel, we denote $\bm{d}=A\bm{1}=(d_{i})\in\mathbb{R}^{n}$ and $D=\hbox{diag}(\bm{d})\in\mathbb{R}^{n\times n}$.

Since $G$ is connected, then $A$ is irreducible, so that $d_{i}>0$ for any $i$, and we may define the irreducible stochastic matrix $P=D^{-1}A$, whose entries are the transition probabilities in a random walk over the graph $G$. The transition probability matrix $P$ defines the Markov chain that models this random walk. Denote by $\lambda_{i}$, $i=1,\ldots,n$, the eigenvalues of $P$. Since $P$ is stochastic, then $|\lambda_{i}|\leq 1$. Moreover, since $D^{\frac{1}{2}}PD^{-\frac{1}{2}}=D^{-\frac{1}{2}}AD^{-\frac{1}{2}}$ is symmetric, the eigenvalues of $P$ are real and can be ordered as

[TABLE]

where the latter strict inequality follows from the irreducibility of $P$ and the Perron-Frobenius theorem.

The Kemeny constant $\kappa(G)$ of $G$ represents the average travel time in the network described by the graph [13]. In other words, $\kappa(G)$ measures the difficulty of circulating in the graph $G$. An algebraic expression is given by

[TABLE]

and this expression is valid independently of the choice of the vector $\bm{v}$ [22, 27].

Observe that, choosing $\bm{v}=\frac{1}{\|\bm{d}\|_{1}}\bm{d}$, we have

[TABLE]

where $\bm{d}^{\frac{1}{2}}=D^{\frac{1}{2}}\bm{1}$. That is, the matrix $I-P+\frac{1}{\|\bm{d}\|_{1}}\bm{1}\bm{d}^{\top}$ is similar to the symmetric matrix $H$, where

[TABLE]

Note that $S$ is a rank-1 perturbation of the graph Laplacian $L$. This way, we may write

[TABLE]

Since $(I-P)\bm{1}=0$, from the Brauer theorem [10], we deduce that the eigenvalues of $I-P+\bm{1}\bm{v}^{\top}$ coincide with the eigenvalues of $I-P$, i.e., $1-\lambda_{i}$, except for the eigenvalue $1-1=0$, which is replaced by 1. Therefore the eigenvalues of $H^{-1}$ are $\gamma_{i}=\frac{1}{1-\lambda_{i}}$, for $i=1,2,\ldots,n-1$ and $\gamma_{n}=1$, so that $\kappa(G)$ has the equivalent expression

[TABLE]

It is convenient to write $\kappa(P)$ in place of $\kappa(G)$ if we want to point out the dependence of Kemeny’s constant on the matrix $P$.

We recall the definition of the Kemeny-based centrality $c(p,q)$ of the edge $\{p,q\}$ of the graph $G$ given in [3]. Let $A_{pq}=A+a_{pq}\bm{w}_{pq}\bm{w}_{pq}^{\top}$, where $\bm{w}_{pq}=\bm{e}_{p}-\bm{e}_{q}$. The matrix $A_{pq}$ is the adjacency matrix of the graph $G_{pq}$ obtained from $G$ by removing the edge $\{p,q\}$ and adding the two loops $\{p,p\}$, $\{q,q\}$ both with weight $a_{pq}$. In [3] it is proved that if $\{p,q\}$ is not a cut edge, $c(p,q)$ defined as

[TABLE]

is always positive. If $\{p,q\}$ is a cut-edge, then the matrix $S_{pq}$ is not invertible, and definition (2.5) does not apply. Theorem 4.3 in [3] states that if $\{p,q\}$ is not a cut-edge then

[TABLE]

In the same paper [3] the regularized centrality measure $c_{r}(p,q)$ is defined for each $r>0$ as

[TABLE]

This formula is also valid for cut-edges since $r>0$, and for a non-cut-edge $\lim_{r\to 0}c_{r}(p,q)=c(p,q)$. For a cut-edge $\{p,q\}$, the same authors introduce the filtered measure defined as

[TABLE]

and show that the limit exists and is positive. We refer to [9] for a detailed analysis of this quantity.

3 Directional derivative of Kemeny’s constant

For $p<q$ such that $a_{pq}\neq 0$, i.e., the edge $e=\{p,q\}$ exists in the graph, consider the new graph ${G}_{pq}(t)$ obtained by changing the weight $a_{pq}$ of the edge $e$ to $(1-t)a_{pq}$ and by adding a loop in vertices $p$ and $q$ with weight $ta_{pq}$, where $0<t\leq 1$. The adjacency matrix $A_{pq}(t)$ of $G_{pq}(t)$ can be obtained by means of a rank-one correction of $A$, that is,

[TABLE]

Since $\bm{w}_{pq}^{\top}\bm{1}=0$, then $A_{pq}(t)\bm{1}=A\bm{1}=\bm{d}$, so that the stochastic matrix $P_{pq}(t)$ associated with $G_{pq}(t)$ is $P_{pq}(t)=D^{-1}A_{pq}(t)$, and the symmetric matrices $H_{pq}(t)$, $S_{pq}(t)$, and $L_{pq}(t)$ corresponding to $H$, $S$, and $L$ in (2.3), are

[TABLE]

Observe that for $t=1$, the graph $G_{pq}(t)$ is obtained by removing the edge $\{p,q\}$ and by adding a loop in the vertex $p$ and one in the vertex $q$ of weight $a_{pq}$. This way, for $t=1$, the difference $\kappa(P_{pq}(t))-\kappa(P)$ coincides with the centrality measure $c(p,q)$ of the edge $\{p,q\}$ of formula (2.5) introduced in [3].

On the other hand, observe that the ratio

[TABLE]

is nothing else but the incremental ratio of the function $\kappa(P_{pq}(t))$, along the direction $a_{pq}D^{-1}\bm{w}_{pq}\bm{w}_{pq}^{\top}$, at $P_{pq}(0)=P$ with incremental step $t$, and its limit for $t\to 0^{+}$ is the directional derivative of $\kappa(P)$.

Observe also that the stochastic matrix $P_{pq}(t)=D^{-1}A_{pq}(t)$ can be expressed as a convex combination of the stochastic matrices $P$ and $P_{pq}(1)$, that is, $P_{pq}(t)=(1-t)P+tP_{pq}(1)$. A similar formulation of $P_{pq}(t)$ is given in [7] to define a directional derivative of Kemeny’s constant along a different suitable direction that is used there to split the graph into connected components.

We introduce a new measure of the edge $\{p,q\}$ given by the directional derivative of $\kappa(G_{pq}(t))$ as a function of $t$, namely

[TABLE]

We will show that $\mu(p,q)$ is well defined and positive for any edge $\{p,q\}$, even when $\{p,q\}$ is a cut-edge. We also introduce an unweighted measure, obtained by computing the derivative along the unweighted direction $D^{-1}\bm{w}_{pq}\bm{w}_{pq}^{\top}$, namely,

[TABLE]

We will show that this latter measure is positive also in the case of a non-edge, $a_{pq}=0$, and that

[TABLE]

We refer to $\mu(p,q)$ and $\bar{\mu}(p,q)$ as the weighted and unweighted centrality measure, respectively.

The following theorem provides an explicit expression of $\kappa(P_{pq}(t))$ and of all its derivatives $\kappa^{(i)}(P_{pq}(t))$ with respect to $t$, given in terms of the quantities

[TABLE]

it will be used to find explicit expressions for the new measures $\mu(p,q)$ and $\bar{\mu}(p,q)$.

Theorem 1**.**

Assume that the graph $G$ is connected. Then we have

[TABLE]

where $\alpha(p,q)$ and $\beta(p,q)$ are defined in (3.5) and $\kappa^{(i)}(P_{t})$ denotes the $i$-th derivative of $\kappa(P_{pq}(t))$ with respect to $t$. In particular, we have $\kappa^{\prime}(P_{pq}(t))|_{t=0}=a_{pq}\beta(p,q)$.

Proof.

If $a_{pq}=0$ the result is trivial. Otherwise, applying the Sherman-Woodbury-Morrison formula to $S_{pq}(t)$ in (3.2) yields

[TABLE]

whence we get

[TABLE]

which provides the first part of the theorem. The second part, concerning derivatives of $\kappa(P_{pq}(t))$, follows by taking the derivatives of the rational function $ta_{pq}/(1-ta_{pq}\bm{w}_{pq}^{\top}S^{-1}\bm{w}_{pq})$. ∎

In the case of the unweighted increment, where in the expression (3.1) of $A_{pq}(t)$, the increment $ta_{pq}\bm{w}_{pq}\bm{w}_{pq}^{\top}$ is replaced by $t\bm{w}_{pq}\bm{w}_{pq}^{\top}$, we get analogous expressions to those of $\kappa(P_{pq}(t))$ and its derivatives where $a_{pq}$ is replaced by 1.

From the above result, we may deduce explicit expressions for $\mu(p,q)$ and $\bar{\mu}(p,q)$ together with their derivatives.

Corollary 1**.**

Under the assumption of Theorem 1, we have

[TABLE]

where $\alpha(p,q)$ and $\beta(p,q)$ are defined in (3.5).

In other words, the entries of the matrix $\beta(p,q)$, defined in (3.5), where $p,q=1,\ldots,n$, provide the values of the centrality measure of the edges and the non-edges (in the case of the unweighted measure $\bar{\mu}(p,q)$) of a graph.

Relying on the symmetry of $H$ and $S$, we may deduce lower and upper bounds to the value of $\bar{\mu}(p,q)$. In this regard, we have the following corollary.

Corollary 2**.**

Under the hypotheses of Theorem 1, we have

[TABLE]

where $\zeta_{i}$, $\lambda_{i}$ and $\nu_{i}$ are the eigenvalues of $S^{-1}DS^{-1}$, $P$ and $S=D-A+\frac{1}{\|\bm{d}\|_{1}}\bm{d}\bm{d}^{\top}$, respectively.

Proof.

From (3.5) and Corollary 1, since $\bm{w}_{pq}^{\top}\bm{w}_{pq}=2$, we get

[TABLE]

Since $\bm{w}_{pq}^{\top}\bm{w}_{pq}=2$, this implies the first pair of inequalities. Moreover, setting $\bm{y}=S^{-1}\bm{w}_{pq}$ we get

[TABLE]

The first Rayleigh quotient is in the interval $[\min_{i}d_{i},\max_{i}d_{i}]$, the second one is in the interval $[\min_{i}\left(\nu_{i}^{-2}\right),\max_{i}\left(\nu_{i}^{-2}\right)]$, where $\nu_{i}$ are the eigenvalues of $S$. From this, we obtain the first two inclusions. To prove the last inclusion, we rely on the identity $S=D^{\frac{1}{2}}HD^{\frac{1}{2}}$ of (2.3), so that we may write $\bar{\mu}(p,q)=\bm{w}_{pq}^{\top}(D^{-\frac{1}{2}}H^{-2}D^{-\frac{1}{2}})\bm{w}_{pq}$. Applying the same argument as before yields

[TABLE]

The first Rayleigh quotient is in the range $[\min\gamma_{i}^{2},\max_{i}\gamma_{i}^{2}]$, where $\gamma_{i}$ are the eigenvalues of $H^{-1}$. Recall that $H$ is similar to $D^{-1}S=I-P+\frac{1}{\|\bm{d}\|_{1}}\bm{1}\bm{d}^{\top}$ so that $\gamma_{i}=\frac{1}{1-\lambda_{i}}$, $i=1,\ldots,n-1$, $\gamma_{n}=1$. This concludes the proof. ∎

An interesting observation stems from comparing the expression of $\mu(p,q)$ given in Corollary 1 with the centrality measure given in [3]. Indeed, Theorem 4.3 in [3], rephrased in our notation, states that if $\{p,q\}$ is not a cut-edge, then the centrality measure given in [3] can be rewritten as in (2.6). This expression of $c(p,q)$ allows us to relate $c(p,q)$ to $\mu(p,q)$ as shown in the following theorem.

Theorem 2**.**

Assume that the graph $G$ is connected, $a_{pq}\neq 0$, and $\{p,q\}$ is not a cut-edge. Then, we have

[TABLE]

Proof.

Rewriting $S_{pq}=(I-a_{pq}\bm{w}_{pq}\bm{w}_{pq}^{\top}S^{-1})S$ and applying the Sherman-Woodbury-Morrison formula yields

[TABLE]

Plugging the latter expression into (2.6), because of Theorem 1, yields the result. ∎

Remark 1*.*

As shown in [3], if $\{p,q\}$ is a cut-edge, then $c(p,q)=\infty$. Therefore, from Theorem 2 we find that, for a cut-edge $\{p,q\}$, necessarily one has $a_{pq}\bm{w}_{pq}^{\top}S^{-1}\bm{w}_{pq}=1$.

We may also relate the measure $\mu(p,q)$ with the filtered centrality $c^{F}(p,q)$ defined in [3], see (2.8).

Consistently with the definition of $c_{r}(p,q)$, we define the regularized derivative-based centrality $\mu_{r}(p,q)$ as

[TABLE]

so that $\mu_{0}(p,q)=\mu(p,q)$. A simple formal manipulation relying on the second equation in (2.5) and the Sherman-Woodbury-Morrison formula shows that

[TABLE]

so that from (2.7) we obtain

[TABLE]

The above expression allows us to relate the filtered centrality measure $c^{F}(p,q)$, defined in (2.8), to our measure $\mu(p,q)$. We have the following

Theorem 3**.**

If $G$ is connected, $a_{pq}\neq 0$, and $\{p,q\}$ is a cut-edge, then

[TABLE]

Proof.

Since $S_{r}^{-1}=S^{-1}(I+rDS^{-1})^{-1}=S^{-1}-rS^{-1}DS^{-1}+O(r^{2})\doteq S^{-1}-rS^{-1}DS^{-1}$, where $\doteq$ means equality up to $O(r^{2})$ terms, then we may write

[TABLE]

where the last equality follows from Remark 1. Whence, in view of (3.6) we get

[TABLE]

Now, it remains to expand $\mu_{r}(p,q)$ as a power series in $r$. Since $(rD+S)^{-1}\doteq S^{-1}-rS^{-1}DS^{-1}$, we have

[TABLE]

This completes the proof. ∎

Observe that, setting $\bm{y}_{pq}=D^{\frac{1}{2}}S^{-1}\bm{w}_{pq}$ we get

[TABLE]

moreover, the Rayleigh quotient $\bm{y}_{pq}^{\top}H\bm{y}_{pq}/\bm{y}_{pq}^{\top}\bm{y}_{pq}$ is in between the maximum and the minimum eigenvalue of $H^{-1}=D^{\frac{1}{2}}S^{-1}D^{\frac{1}{2}}$. Recalling (2.1), these eigenvalues are $\gamma_{i}=\frac{1}{1-\lambda_{i}}$, $i=1,\ldots,n-1$, $\gamma_{n}=1$, and we get the following bounds.

Corollary 3**.**

If $A$ is irreducible and $\{p,q\}$ is a cut-edge, then

[TABLE]

3.1 Condition number

Observe that in the graph with adjacency matrix $A-t\bm{w}_{pq}\bm{w}_{pq}^{\top}$, the edge $\{p,q\}$ has a weight $a_{pq}+t$ and the edges $p$ and $q$ have weights $a_{pp}-t$ and $a_{qq}-t$, respectively. The derivative with respect to $t$ of the Kemeny constant of this graph expresses the absolute variation consequent to the absolute perturbation $t$ of $a_{pq}$ and $a_{qp}$ and to the perturbation $-t$ in $a_{pp}$ and in $a_{qq}$, i.e., in the loops in the nodes $p$ and $q$. This way, we may look at this derivative at $t=0$ as the condition number of the Kemeny constant with respect to an absolute perturbation on $a_{pq}$, when this perturbation is applied also to the loops in $p$ and $q$ with negative weight. Relying on Theorem 1, it is easy to verify that this derivative at $t=0$ is given by $-\beta(p,q)=-\bm{w}_{pq}^{\top}(S^{-1}DS^{-1})\bm{w}_{pq}$. This shows that increasing the weight of $a_{pq}$ and $a_{qp}$, and decreasing by the same amount the weights $a_{pp}$ and $a_{qq}$ leads to a decrease of Kemeny’s constant. Moreover, the larger $\beta(p,q)$, the stronger is the variation, and consequently the more ill-conditioned is Kemeny’s constant.

3.2 A global measure of sensitivity

By following the way the total communicability is defined [6, 4], we can take the arithmetic mean of the measure $\bar{\mu}(p,q)$ over each pair $(p,q)$ (including non-edges and self-loops, where we assume $\bar{\mu}(p,p)=0$), namely

[TABLE]

as a global measure of the sensitivity of Kemeny’s constant of a graph $G$ to perturbations of its edges. A direct computation based on Theorem 1 shows that

[TABLE]

Observe that the above expression involves the trace and the sum of the entries of the matrix $S^{-1}DS^{-1}$; the trace resembles the Estrada index [4], while the sum of the entries resembles total communicability [6].

4 Computational issues

This section discusses the issues encountered in the actual computation of the quantity $\mu(p,q)$.

From Corollary 1 and equation (3.5), we deduce that

[TABLE]

where $\bm{w}_{pq}=\bm{e}_{p}-\bm{e}_{q}$, the matrix $S$ is defined in (2.4), and $\sigma_{ij}$ are the entries of $S^{-1}$. Moreover,

[TABLE]

The major effort in computing $\mu(p,q)$ through (4.1) consists in computing the vector $\bm{x}=S^{-1}\bm{w}_{pq}$, that is, solving the system $S\bm{x}=\bm{w}_{pq}$. If the goal is to compute the centrality measure of a single edge $\{p,q\}$, then one can compute the Cholesky factorization $S=R^{\top}R$ of $S$, and then solve the two triangular systems with matrix $R^{\top}$ and $R$, respectively. The cost of this approach is dominated by the computation of the Cholesky factorization of $S$, that is, $\frac{1}{3}n^{3}+O(n^{2})$ arithmetic operations (ops) for a general matrix $A$ [19, Appendix C]. Alternatively, one can apply any (preconditioned) iterative method to solve the system $S\bm{x}=\bm{w}_{pq}$, taking advantage of the possibly sparse structure of $A$. However, the cost of the former approach based on Cholesky factorization can be reduced by exploiting the possible sparsity of the matrix $A$ as follows.

First, a symmetric sparse matrix can be reduced to a banded form by a permutation of rows and columns [14]. Therefore, without loss of generality, we may assume that $A$ is banded with bandwidth $2b-1$, i.e., $a_{ij}=0$ if $|i-j|\geq b$. Secondly, we may exploit that $L=D-A$ is a weakly diagonally dominant and irreducible M-matrix that is also singular. This fact implies that its leading principal submatrix $L_{n-1}$ of $L$, with size $n-1$, is positive definite; therefore, it admits the Cholesky factorization $L_{n-1}=R_{n-1}^{\top}R_{n-1}$ where $R_{n-1}$ is an $(n-1)\times(n-1)$ upper triangular matrix.

Thus, we may write the Cholesky factorization of $L$ as

[TABLE]

This expression allows us to compute the entries $v_{1},v_{2},\ldots,v_{n-1}$ of the last column the Cholesky factor $R$ of $L$ by solving the linear system

[TABLE]

where $(\ell_{ij})_{ij}$ denote the entries of $L$.

Now, the system $S\bm{x}=\bm{w}_{pq}$ can be equivalently rewritten as

[TABLE]

that is,

[TABLE]

On the other hand, since $\bm{1}^{\top}L=0$ and $\bm{1}^{\top}\bm{w}_{pq}=0$, by multiplying on the left (4.4) by $\bm{1}^{\top}$, we get $\xi=\bm{d}^{\top}\bm{x}=0$, so that the system (4.4) turns into

[TABLE]

Let $\bm{y}\in\mathbb{R}^{n-1}$ be the vector solving the nonsingular system

[TABLE]

From the condition $R^{\top}R\bm{x}=\bm{w}_{pq}$ we get \left[\begin{array}[]{c|c}R_{n-1}&\bm{v}\\ \end{array}\right]\bm{x}=\bm{y}, so that, adding the condition $\bm{d}^{\top}\bm{x}=0$, we may rewrite the system (4.5) as

[TABLE]

where $\bm{d}_{n-1}=(d_{1},\ldots,d_{n-1})^{\top}$. Setting $\bm{z}$ the solution of the system

[TABLE]

we arrive at

[TABLE]

We may conclude with the following expression for the solution $\bm{x}^{\top}=[\bm{x}_{n-1}^{\top}\ x_{n}]$:

[TABLE]

The overall computation is summarized in Algorithm 1.

Remark 2*.*

Since $L=(\ell_{ij})_{ij}$ is banded, then $\ell_{in}=0$ for $i=1,\ldots,n-b$, and since $R_{n-1}^{\top}$ is lower triangular, then $v_{i}=0$ for $i=1,\ldots,n-b$ and the vector $(v_{n-b+1},\ldots,v_{n-1})^{\top}$ solves a $(b-1)\times(b-1)$ triangular system.

Remark 3*.*

Since $\bm{w}_{pq}=\bm{e}_{p}-\bm{e}_{q}$, $p<q$, then the first $p-1$ components of $\bm{w}_{pq}$ are zero, and since $R_{n-1}^{\top}$ is lower triangular, then $w_{i}=0$ for $i=1,\ldots,p-1$ and the linear system (4.6) reduces to an $(n-p)\times(n-p)$ triangular system.

We can perform a complexity analysis of Algorithm 1 by evaluating the number of arithmetic operations required for the computation. In this analysis, we assume that the matrix $A$ has bandwidth $2b-1$.

Concerning the preprocessing stage, computing $\bm{d}$ costs $(b-1)(2n-(b-1))$ arithmetic operations (ops), while at most $n$ subtractions are needed to compute $L$. The Cholesky factorization of a positive semidefinite matrix, with bandwidth $2b-1$, costs $\frac{1}{3}nb^{2}$. We refer to the book [18, Chapter 10] for the nice computational and numerical features of the Cholesky factorization. We also point out that in our case, the matrix $L_{n-1}$ is a nonsingular M-matrix. This feature adds nice numerical properties to the algorithm for Cholesky factorization if complemented with the GTH trick [15] to avoid possible numerical cancellation.

The vector $\bm{v}$ in the preprocessing stage can be computed by solving a triangular system of size $b$ (see Remark 2), for the cost of $b^{2}$ ops, while the computation of $\bm{z}$ requires the solution of the triangular system (4.7) of size $n-1$ with bandwidth $b$ for the cost of $n(2b-1)-b^{2}+b$ ops. Finally, computing $\rho=\bm{z}^{\top}\bm{v}-d_{n}$ costs just $2b-1$ ops since $\bm{v}$ has only the last $b-1$ components nonzero. Thus, the overall cost of the preprocessing stage is $n(\frac{1}{3}b^{2}+4b-3)-b^{2}+5b-2$ ops and is dominated by the cost of Cholesky factorization.

Concerning the processing stage, because of Remark 3, step 1 requires the solution of an $(n-p)\times(n-p)$ triangular banded system that costs roughly $2b(n-p)$ ops; step 2 requires computing the scalar product $\bm{z}^{\top}\bm{y}$ which costs $2(n-p)$ ops in view of Remark 3. Step 3 requires solving an upper triangular system of size $n$ with bandwidth $b$; this costs roughly $2bn$ ops. That is, the overall cost of the processing stage, excluding preprocessing, is $2(n-p)(b+1)$ operations.

To compute the centrality of the edge $\{p,q\}$, it is enough to compute $a_{pq}\bm{x}^{\top}D\bm{x}$ that requires $3n$ additional ops. Therefore, the overall cost of computing $\mu(p,q)$, once the preprocessing stage has been performed, is $n(4b+5)-2p(b+1)$ ops.

To compute the centralities of all the edges, one can repeat the processing part for all the pairs $\{p,q\}$ corresponding to an edge. An upper bound to this computational cost is $n(4b+3)|E|+\frac{1}{3}nb^{2}$, where $|E|<nb$ is the number of edges in the graph.

In Section 5, we will provide explicit expressions of $\mu(p,q)$ for one-path graphs as performed in [9] for $c(e)$, and show the good behaviour of the new measure.

5 One-path graphs

Consider the case of a one-path graph formed by $n$ vertices connected as in Figure 1.

The adjacency matrix $A=(a_{ij})$ is tridiagonal with $a_{i,i+1}=a_{i+1,i}\neq 0$, $a_{ij}=0$, otherwise. Let $P=D^{-1}A$, where $D=\hbox{diag}(\bm{d})$, $\bm{d}=A\bm{1}$, then

[TABLE]

where $\theta_{0}=1-\lambda_{0}\geq 0$, $\theta_{i}=1-\lambda_{i}-\nu_{i}\geq 0$ for $i=1,\ldots,n-2$, $\theta_{n-1}=1-\lambda_{n-1}$, and $0<\lambda_{i}<1$, $0<\nu_{i}<1$. For this graph, we have $\theta_{i}\neq 0$ if and only if there is a loop $\{i,i\}$ at node $i$. From [8, Theorem 5.1], we find that Kemeny’s constant of $P$ is given by

[TABLE]

where $\bm{\pi}=(\pi_{i})_{i=0,\ldots,n-1}$ is the steady state vector of $P$, i.e., the stochastic solution of $\pi^{\top}P=\pi^{\top}$. Recall also that

[TABLE]

Now consider an integer $q$, $1\leq q\leq n-1$, and define

[TABLE]

Observe that $A_{q}(t)\bm{1}=\bm{d}$ so that denoting $P_{q}(t)=\hbox{diag}(\bm{d})^{-1}A_{q}(t)$ we have

[TABLE]

where $\lambda_{i}(t)=\lambda_{i}$, $\nu_{i+1}(t)=\nu_{i+1}$ are independent of $t$ for $i\neq q-1$, while $\lambda_{q-1}(t)=\lambda_{q-1}(1-t),\quad\nu_{q}(t)=\nu_{q}(1-t)$. These expressions, together with (5.3) imply that $\bm{\pi}^{\top}P_{q}(t)=\bm{\pi}^{\top}$, consequently

[TABLE]

The $k$th term in the above summation is independent of $t$ if $k\neq q-1$ while the $(q-1)$st term is

[TABLE]

so that

[TABLE]

We conclude with the following theorem.

Theorem 4**.**

For the centrality $\mu(q,q+1)$ of the edge $\{q,q+1\}$ in the one-path graph associated with the matrix $P$ of (5.1), we have

[TABLE]

where $\sigma_{k}=\sum_{j=0}^{k}\pi_{j}$, $\pi^{\top}P=\pi^{\top}$, and the components $\pi_{i}$ of $\pi$ satisfy (5.3). Moreover, in the case of unitary weights, we have $\lambda_{0}=\nu_{n-1}=1$, $\lambda_{q}=\nu_{n-q-1}=1/2$ for $q>0$, so that

[TABLE]

These values give

[TABLE]

It is interesting to observe that the values given in (5.5) are located on the graph of a parabola as the values of the measure given in [9, Theorem 7], which essentially differ for a normalizing factor and lower order terms.

6 Numerical tests

We have performed some experiments to test our centrality measures. A first test concerns the good physical behavior of the measure for assessing the importance of roads in a given road map. In Figure 2, we compare our measure to the measure of [3] on the road map of Pisa. The two measures have been linearly normalized so that their values are in the range $[0,1]$, i.e., the transformation $x\to(x-a)/(b-a)$ has been applied, where $a$ and $b$ are the minimum and the maximum values taken by the measure. Regarding the association of the color with the numerical value of the measure, we have considered the square root of the measure to better highlight the intermediate values.

From this comparison, we can see that the two measures behave similarly; in particular, both highlight the most important connections, such as the bridges on the river Arno and the railway bridges. Moreover, the relative differences of the two measures after linear normalization are between $7.2\times 10^{-8}$ and $2.8\times 10^{-3}$. Additionally, the histograms of the two measures in Figure 3, which report, on a log scale in the $y$ axis, the number of values falling in each of 50 bins, appear very similar. However, the new measure assigns fewer edges to the rightmost bins, which corresponds to larger values, providing a more definite ranking of the most important roads. This clarifies that the bridges over the Arno river are the most crucial, while the less connected peripheral railroad bridges are assigned a lower importance. We believe this new measure offers a more realistic and useful importance rating.

Using MATLAB on a 12-core, 3.4 GHz parallel server, our new implementation significantly reduces computation time. For a smaller dataset with n=1404 nodes and m=1794 edges, the computation takes about 0.1 seconds. For the entire Tuscany region ($n=1,150,138$, $m=1,223,248$), the computation is completed in 62 minutes, a significant improvement over the 18 hours reported in the implementation given in [3]. This substantial speed-up is primarily due to our enhanced code implementation.

A second test concerns the computation of the centrality measure $\bar{\mu}(p,q)$ of the edges and the “non-edges”. In Figure 4, the cases of some specific graphs are considered. Namely, star graphs, one-path graphs, and cyclic graphs are analyzed with size $n=10$, together with a binary tree graph with 7 nodes. On the left, the graph picture is displayed, on the right, the values of the centralities of edges and non-edges are shown, accompanied by a background color that highlights the centrality level of the corresponding edge.

In the case of the star graph, all the non-edges have necessarily the same centrality value, which is slightly higher than twice the value of the existing edges. For the one-path graph, the most important non-edge is the one that connects the first node to the last one, and the importance of non-edges is roughly proportional to the number of nodes that they allow to ’skip’ over. Similarly, for the cyclic graph, the most important non-edges are the ones that connect a node to its direct opposite. Also, notice that the matrix with the centralities of edges and non-edges is circulant, due to the cyclicity of the problem. For the binary tree, the most important non-edges are those that connect the leaves having a large relative distance.

Evaluating the centrality of the non-edges is also related to the problem of link prediction in social networks [17], where the goal is to predict the likelihood of a future association between two nodes, knowing that there is no association between the nodes in the current state of the graph. For this problem, the directional derivative of Kemeny’s constant with respect to the non-edge $\{p,q\}$ can be interpreted in the following way: the higher the variation, the more unlikely is the possibility that the social relation between node $p$ and node $q$ will be activated.

We compared the measure $\bar{\mu}$ to several existing ones on a link prediction experiment with real-world data: the network of collaboration between network science authors studied in [25]111https://networks.skewed.de/net/netscience. This network is a weighted, undirected graph. We restrict our analysis to its largest connected component, with 379 nodes and 914 edges.

Most existing measures used for link prediction, such as the Jaccard and Adamic-Adar indices [24] and the resource allocation index [29], are based only on common neighbours: the probability of a new link between $p$ and $q$ is higher when they have many common neighbours. As a result, the score of any two nodes with no common neighbours is zero, and these measures cannot predict links between them. The so-called common neighbour centrality [1] is also based on common neighbours, but with a further additive term that is inversely proportional to the distance between the two nodes; hence, this centrality measure gives a nonzero score also to nodes that have no common neighbours. We compare these four indices to our new measure $\bar{\mu}$, using their implementations in the Python library NetworkX [16].

We display in Figure 5 the new predicted edges on a plot of the original graph. The correlation coefficients between the five measures are shown in Table 1.

The correlations between the Kemeny derivative and the other four measures are negative, as expected, because they are sorted in opposite ways: a low Kemeny derivative means a high probability of a new link. The three measures based on common neighbours are highly correlated with each other, and only poorly correlated with the new measure. The new measure is moderately correlated with the common neighbour centrality. The node labels for the top ten new links predicted by the new measure $\bar{\mu}$ are shown in Table 2. We can verify on literature databases that pairs #1, #6, #7 have a common publication that was not included in the original dataset; we consider this as a confirmation that our measure is effective at link prediction.

For the circulant, star, and one-path graphs represented in Figure 4, we have computed the global sensitivity $\zeta(G)$ defined in (3.7), as a function of the number of vertices $n$. From Figure 6, we observe that $\zeta(G)$ is almost constant for the star graph, while it increases very quickly for the circulant and one-path graph. This behavior is expected from the topology of the three kinds of graphs.

7 Conclusions

Edge centrality measures are important for assessing the relevance of roads in a road map. While the measure introduced in [3] is highly effective, its definition and computation must be adjusted for cut-edges [3, 9]. If not calculated correctly, this measure is susceptible to numerical instability.

We developed a novel, always-positive measure for graph centrality based on the directional derivative of Kemeny’s constant. This measure is applicable to all edges, including cut-edges, and avoids the computational issues of cancellation. We have derived an explicit formula for this measure, which relates it to existing measures from [3] and [9] through the inverse of a modified Laplacian matrix.

Our measure can also be computed for non-edges, making it useful for link prediction. We have provided explicit expressions for one-path graphs, showing our measure’s relationship to edge weights. Numerical experiments confirm its effectiveness as a tool for assessing road map centrality and as a strong alternative for link prediction analysis.

Funding

This work has been partially supported by: MUR Excellence Department Project awarded to the Department of Mathematics, University of Pisa, CUP I57G220007

00001; European Union - NextGenerationEU under the National Recovery and Resilience Plan (PNRR) - Mission 4 Education and research - Component 2 From research to business - Investment 1.1 Notice Prin 2022 - DD N. 104 2/2/2022, titled “Low-rank Structures and Numerical Methods in Matrix and Tensor Computations and their Application”, proposal code 20227PCCKZ – CUP I53D23002280006. The second and the third author are member of the INdAM GNCS group.