Capturing the interplay of dynamics and networks through parameterizations of Laplacian operators

Random walks

One of the most-widely studied dynamical processes on networks is the random walk. The simplest is the discrete time unbiased random walk (URW), where a walker at vertex i follows one of the edges with a probability proportional to the weight of the edge (Ross, 2014; Aldous & Fill, 2002). In this case, the state vector $\mathbf{\theta }$⁴ forms a distribution whose expected value follows the update equation: ${\displaystyle {\theta }_i\left(t+1\right)=\sum _j{P}_{ij}{\theta }_j\left(t\right).}$ Here P is a stochastic matrix whose entry P_ij is the transition probability for a walker to go from the vertex j to i, P_ij = a_ij∕d_j.

The update equation of an unbiased random walk leads to the difference equation ${\displaystyle \Delta {\theta }_i={\theta }_i\left(t+1\right)-{\theta }_i\left(t\right)=\sum _j{P}_{ij}{\theta }_j\left(t\right)-{\theta }_i\left(t\right)=-\sum _j{L}_{ij}^{RW}{\theta }_j\left(t\right),}$ where L^RW is the normalized random walk Laplacian matrix with $L^{RW}=\mathit{I}-\mathit{A}{\mathit{D}}_{\mathit{A}}^{-1}$.

To go from a discrete time synchronous random walk to a continuous time dynamics, we introduce a waiting time function for the asynchronous jumps performed by the walk (Ross, 2014). Assuming a simple Poisson process where the waiting times between jumps are exponentially distributed as the PDF $f\left(t,\tau \right)=\frac{1}{{\tau }_i}{e}^{-\frac{t}{{\tau }_i}}$, we can rewrite the above difference equations as differential equations, ${\displaystyle \frac{d{\theta }_i}{dt}=-\sum _j\frac{L_{ij}^{RW}}{{\tau }_j}{\theta }_j.}$ The solution to the above differential equations gives the state vector of the random walk at any time t: ${\displaystyle \mathbf{\theta }\left(t\right)=e^{-L^{RW}{\mathit{T}}^{-1}t}\cdot \mathbf{\theta }\left(0\right),}$ where $\mathit{T}$ is the n × n diagonal matrix with the mean waiting time τ_i as entries. If the dynamical process converges, then regardless of its initial value $\mathbf{\theta }\left(0\right)$, the stationary distribution π_i has the following density: (1)${\displaystyle {\pi }_i=\underset{t\to \infty }{lim}{\theta }_i\left(t\right)=\frac{d_i{\tau }_i}{\sum _j{d}_j{\tau }_j}.}$ Intuitively, the stationary distribution is proportional to the product of vertex degree and the mean waiting time.

A natural extension of the process is to bias the random walk towards specific vertices, making it a biased random walk (BRW). According to Lambiotte et al. (2011), any biased random walk defined with the transition probability P_ij ∝ b_ia_ij (where b_i is the bias towards vertex i) can be reduced to a URW on a re-weighted “interaction network” with the adjacency matrix ${\displaystyle \mathit{W}=\mathit{B}\mathit{A}\mathit{B},}$ where $\mathit{B}$ is a diagonal matrix with ${\mathit{B}}_{ii}=b_i$. The above symmetric re-weighting ensures that ${\displaystyle P_{ij}=\frac{b_i{a}_{ij}{b}_j}{\sum _i{b}_i{a}_{ij}{b}_j}\propto b_i{a}_{ij},P_{ji}=\frac{b_j{a}_{ji}{b}_i}{\sum _j{b}_j{a}_{ji}{b}_i}\propto b_j{a}_{ji}.}$

In one class of BRWs previously studied in network communications (Ling et al., 2013; Fronczak & Fronczak, 2009; Gómez-Gardeñes & Latora, 2008), bias b_i has a power-law dependence on degree: $P_{ij}\propto d_i^{\beta }{a}_{ij}$. The exponent β controls the strength of bias. The URW is recovered with β = 0; When β > 0, biases toward high degree vertices are introduced, and when β < 0, the random walk is more likely to jump to a lower degree neighbor.

Another type of BRW is the maximum-entropy random walk (Burda et al., 2009; Lambiotte et al., 2011), defined as ${\displaystyle {\theta }_i\left(t+1\right)=\sum _j\frac{\stackrel{⃗}{v_{{\mathit{A}}_i}}{a}_{ij}}{{\lambda }_{max}\stackrel{⃗}{v_{{\mathit{A}}_j}}}{\theta }_j\left(t\right),}$ where $\overrightarrow{v_{\mathit{A}}}$ is the eigenvector of $\mathit{A}$ associated with its largest eigenvalue λ_max: $\mathit{A}\overrightarrow{v_{\mathit{A}}}={\lambda }_{max}\overrightarrow{v_{\mathit{A}}}$. Again, an unbiased random walk on the interaction network $\mathit{W}={\mathit{V}}_{\mathit{A}}\mathit{A}{\mathit{V}}_{\mathit{A}}$ is equivalent to biased random walk on the original network $\mathit{A}$ (the entries of diagonal matrix ${\mathit{V}}_{\mathit{A}}$ is the components of the eigenvector $\overrightarrow{v_{\mathit{A}}}$). In particular, the stationary distributions of both can be written as ${\pi }_i=\frac{\stackrel{⃗}{v_{{\mathit{A}}_i}^2}}{{\sum }_i{\stackrel{⃗}{v_{{\mathit{A}}_i}}}^2}.$

Consensus and opinion dynamics

Another closely related class of discrete time dynamical processes is the so-called the “consensus process” (DeGroot, 1974; Lambiotte et al., 2011; Olfati-Saber, Fax & Murray, 2007; Krause, 2008). Consensus process models coordination across a network where each vertex updates its “belief” based on the average “beliefs” of its neighbors. Unlike random walks, which conserves total state value throughout the network (since the state vector is always a distribution), the consensus process follows the following update equation ${\displaystyle {\theta }_i\left(t+1\right)=\frac{1}{d_i}\sum _j{a}_{ij}{\theta }_j\left(t\right).}$ This leads to the difference equation ${\displaystyle \Delta {\theta }_i={\theta }_i^{t+1}-{\theta }_i^t=-\sum _j{L}_{ij}^{CON}{\theta }_j\left(t\right)}$ where L^CON is the consensus Laplacian matrix with $L^{CON}=\mathit{I}-{\mathit{D}}_{\mathit{A}}^{-1}\mathit{A}$. For an undirected graph with a symmetric $\mathit{A}$, L^CON = [L^RW]^T.

Consensus can also be turned into asynchronous continuous time dynamics. Again, assuming a Poisson process where the update interval at each vertex is exponentially distributed as ${\mathrm{\tau }}_i\left(t\right)=\frac{1}{{\tau }_i}{e}^{-\frac{t}{{\tau }_i}}$, we can rewrite the above difference equations as differential equations, ${\displaystyle \frac{d{\theta }_i}{dt}=-\sum _j\frac{L_{ij}^{CON}}{{\tau }_i}{\theta }_j.}$

The consensus process always converge to a uniform “belief” state with the value, (2)${\displaystyle {\pi }_i=\frac{1}{\sum _j{d}_j{\tau }_j}\sum _i{\theta }_i\left(0\right)d_i{\tau }_i.}$

Just like the URW, unbiased consensus can also be generalized by introducing a weight when averaging over neighbors’ values. This opens the door to consensus dynamics such as opinion dynamics (Krause, 2008), and linearized approach to synchronization models (Lerman & Ghosh, 2012; Motter, Zhou & Kurths, 2005; Arenas, Díaz-Guilera & Pérez-Vicente, 2006).

Communities and conductance

In network clustering and community detection, previous work has focused on identifying subsets of vertices S⊆V that interact more frequently with vertices in the same community than vertices in other subsets (Fortunato, 2010; Porter, Onnela & Mucha, 2009). A standard approach to clustering defines an objective function that measures the quality of a cluster. For a subset S⊆V, let $\bar{S}=V\setminus S$ denote the complement of S, which consists of vertices that are not in S. Let $\mathrm{cut}\left(S,\bar{S}\right)={\sum }_{i\in S,j\in \bar{S}}{a}_{i,j}$ denote the total interaction strength of all edges used by S to connect with the outside world. Let $\mathrm{vol}\left(S\right)={\sum }_{i\in S}{d}_i={\sum }_{i\in S,j\in V}{a}_{i,j}$ denote the volume of weighted “importance” for all vertices in S.

One popular measure of the quality of a subset S as a potential good cluster (or a community) (Kannan, Vempala & Vetta, 2004; Spielman & Teng, 2004; Chung, 1997) is to use the ratio of these two quantities: (3)${\displaystyle \varphi \left(S\right)=\frac{\mathrm{cut}\left(S,\bar{S}\right)}{min\left(\mathrm{vol}\left(S\right),\mathrm{vol}\left(\bar{S}\right)\right)}.}$ For example, a subset that (approximately) minimizes this quantity—the conductance of S—is a desirable cluster, as it maximizes the fraction of affinities within the subset. If interactions among vertices are proportional to their affinity weights, then a set with small conductance also means that its members interact significantly more with each other than with outside members. The smallest achievable ratio over all possible subsets is also known as the isoperimetric number. As an important measure for mixing time in classic Markov chains, conductance has proven mathematical bounds in terms of the second eigenvalue of its Laplacian (Cheeger, 1970; Jerrum & Sinclair, 1988; Lawler & Sokal, 1988). Other well-known quality functions are normalized cut (Shi & Malik, 2000) and ratio-cut, given respectively by ${\displaystyle \frac{\mathrm{cut}\left(S,\bar{S}\right)}{\mathrm{vol}\left(S\right)}+\frac{\mathrm{cut}\left(S,\bar{S}\right)}{\mathrm{vol}\left(\bar{S}\right)}\text{and}\frac{\mathrm{cut}\left(S,\bar{S}\right)}{min\left(\left|S\right|,\left|\bar{S}\right|\right)}.}$

Algorithmically, once a quality function is selected, one can then perform a graph partitioning algorithm or any community detection algorithm to find clusters that optimize the objective. The optimization, however, is usually a combinatorial problem. To address this problem on large networks, efficient approximate solutions have been developed, such as Spielman & Teng (2004), Andersen, Chung & Lang (2007), and Andersen & Peres (2009). Others took a machine learning approach, proposing efficient approximations by enforcing various smoothness and regularization conditions (Avrachenkov et al., 2011; Bertozzi & Flenner, 2012).

While most community detection algorithms do not explicitly model the dynamical process that defines the interactions between vertices, the connection between conductance and unbiased random walks is quite well studied (Kannan, Vempala & Vetta, 2004; Spielman & Teng, 2004; Chung, 1997). In particular, Chung’s work on heat kernel page rank and Cheeger inequality, where a dynamical system is built using the normalized Laplacian, provides a theoretical framework for provably good approximations to the isoperimetric number (Chung, 2007). Intuitively, the relationship between clustering and dynamics can be captured as: a community is a cluster of vertices that “trap” a random walk for a long period of time before it jumps to other communities (Lovász, 1996; Shi & Malik, 2000; Rosvall & Bergstrom, 2008; Spielman & Teng, 2004). Therefore, the presence of a good cluster based on conductance implies that it will take a random walk a long time to reach its stationary distribution. Similar interplays with community structures can also be generalized to richer dynamical processes, with different time scale, biases and locality settings (Lambiotte, Delvenne & Barahona, 2008; Lambiotte et al., 2011; Jeub et al., 2015).

Similarity transformations

Changing ρ in Eq. (5) leads to different representations of the same linear operator, unifying seemingly unrelated dynamics, such as “consensus” and “random walk.” To see this, we refer to the idea of matrix similarity.

In linear algebra, similarity is an equivalence relation for square matrices. Two n × n matrices X and Y are similar if (6)${\displaystyle X=QY{Q}^{-1},}$ where the invertible n × n matrix Q is called the change of basis matrix. Similar matrices share many key properties, including their rank, determinant and eigenvalues. Eigenvectors are also transforms of each other under a change of basis.

Recall that under our framework, the symmetric version of the parameterized Laplacian matrix is ${\displaystyle {\mathcal{L}}^{SYM}={\mathit{T}}^{-1∕2}{\mathit{D}}_{\mathit{W}}^{-1∕2}\left({\mathit{D}}_{\mathit{W}}-\mathit{W}\right){\mathit{D}}_{\mathit{W}}^{-1∕2}{\mathit{T}}^{-1∕2}.}$ We can rewrite the operator describing random walk dynamics as: (7)${\displaystyle {\mathcal{L}}^{RW}=\left({\mathit{D}}_{\mathit{W}}-\mathit{W}\right){\left({\mathit{D}}_{\mathit{W}}\mathit{T}\right)}^{-1}={\left({\mathit{D}}_{\mathit{W}}\mathit{T}\right)}^{1∕2}{\mathcal{L}}^{SYM}{\left({\mathit{D}}_{\mathit{W}}\mathit{T}\right)}^{-1∕2}.}$ Thus, continuous time random walk with delay factors $\mathit{T}$ is similar to the symmetric normalized Laplacian. Similarly, we can rewrite the continuous time consensus dynamics under our framework as (8)${\displaystyle {\mathcal{L}}^{CON}={\left({\mathit{D}}_{\mathit{W}}\mathit{T}\right)}^{-1}\left({\mathit{D}}_{\mathit{W}}-\mathit{W}\right)={\left({\mathit{D}}_{\mathit{W}}\mathit{T}\right)}^{-1∕2}{\mathcal{L}}^{SYM}{\left({\mathit{D}}_{\mathit{W}}\mathit{T}\right)}^{1∕2}={{\mathcal{L}}^{RW}}^T.}$ The fact that “consensus,” “symmetric” and “random walk” operators are similar means that they model the same dynamics on a network, provided that we observe them in a consistent basis.

The random walk Laplacian matrix provides a physical intuition for our framework. An unbiased random walk on the interaction graph $\mathit{W}$ is equivalent to a biased random walk on the original adjacency matrix $\mathit{A}$ (Lambiotte et al., 2011). On the other hand, τ_i specifies the mean delay time of the random walk on vertex i before a transition. This interpretation reveals the orthogonal nature of the parameters: namely $\mathit{W}$ controls the distribution of walk trajectories while $\mathit{T}$ controls the delay time of vertex transitions along each trajectory.

While we use symmetric operators for mathematical convenience in definitions and proofs and abuse the notation $\mathcal{L}={\mathcal{L}}^{SYM}$, it is often more intuitive to think from the random walk or consensus perspective. In the following subsections, we will use the random walk formulation (ρ = − 1∕2) as examples, but all results apply to arbitrary ρ values under a simple change of basis. More discussion about the similarity transformation follows after we introduce a few properties of the parameterized Laplacian.

Scaling transformations

Reweighing transformations

The last parameterization we explore is one that transforms the adjacency matrix of a graph, $\mathit{A}$, to the interaction matrix $\mathit{W}$. Given an adjacency matrix $\mathit{A}$, the choice of $\mathit{W}$ is a rather flexible design option. In fact, we can arbitrarily manipulate the adjacency matrix as long as the result is still a positive and symmetric matrix, for any perceived dynamics.

In this paper, we limit our attention to bias transformations of the original adjacency matrix $\mathit{A}$. We call them the reweighing transformations. Whereas the scaling transformation changes the delay time at each vertex, the reweighing transformation changes the trajectory of the dynamic process. Note that this transformation also changes the degree matrix ${\mathit{D}}_{\mathit{W}}$.

As described in ‘Background and Related Work’, a biased random walk with transition probability P_ij ∝ b_ia_ij is equivalent to an unbiased random walk on an “interaction graph,” represented by the reweighed adjacency matrix: (11)${\displaystyle w_{ij}=b_i{a}_{ij}{b}_j,}$ where we constrain b_i > 0. This transformation allows the parameterized Laplacian to model many different types of dynamic processes by transforming them into a unbiased walk on the reweighted interaction graph.

Special cases

The simple parameterization of the Laplacian in terms of $\mathit{T}$ and $\mathit{W}$ allows us to model a variety of dynamic processes,⁵ including those described by the Laplacian and normalized Laplacian, as well as a continuous family of new operators that are not as well studied. It also contains operators for modeling some types of epidemic processes. The consideration of this family of operators is also partially motivated by recent experimental work in understanding network centrality (Ghosh & Lerman, 2011; Lerman & Ghosh, 2012).

Unbiased Laplacian.

Reweighing each edge by the inverse of the square root of the endpoint degrees gives the what is known as the normalized adjacency matrix $\mathit{W}={\mathit{D}}_{\mathit{A}}^{-1∕2}\mathit{A}{\mathit{D}}_{\mathit{A}}^{-1∕2}$ (Chung, 1997). Then, the degree of vertex i of the reweighted graph is ${d_{\mathit{W}}}_i={\sum }_{j\in V}\mathit{W}\left[i,j\right]$. With $\mathit{T}={d_{\mathit{W}}}_{max}{D}_{\mathit{W}}^{-1}$ we define the unbiased Laplacian matrix: ${\displaystyle \mathcal{L}=\frac{1}{{d_{\mathit{W}}}_{max}}\left({\mathit{D}}_{\mathit{W}}-\mathit{W}\right).}$

Unbiased Laplacian is an example of the degree based biased random walk with $P_{ij}\propto d_i^{-1∕2}{a}_{ij}$ (‘Background and Related Work’). An URW on the reweighed adjacency matrix $\mathit{W}$ is equivalent to a BRW on the original adjacency matrix of the following dynamics (13)${\displaystyle {\theta }_i\left(t+1\right)=\sum _j\frac{d_i^{-1∕2}{a}_{ij}}{\sum _k{d}_k^{-1∕2}{a}_{kj}}{\theta }_j\left(t\right).}$

The stationary distribution for this class of BRWs in general is ${\pi }_i=\frac{{\sum }_i{d}_i^{\beta }{a}_{ij}{d}_j^{\beta }}{{\sum }_{ij}{d}_i^{\beta }{a}_{ij}{d}_j^{\beta }}.$

Equivalent to the (scaled) graph Laplacian of the normalized adjacency matrix, the diagonal matrix $\mathit{T}{\mathit{D}}_{\mathit{W}}$ of the unbiased Laplacian is also effectively a scalar. As a result, the “random walk” and “consensus” formulations are exactly the same as the symmetric formulation.

These four special cases are related to each other through various transformations introduced earlier in this section, which are captured by Fig. 1.

Stationary distribution of a random walk

A vertex has high centrality with respect to a random walk if it is visited frequently by it. This is specified by the distribution of the dynamic process at time t: (14)${\displaystyle \mathbf{\theta }\left(t\right)=e^{-{\mathcal{L}}^{RW}t}\cdot \mathbf{\theta }\left(0\right)=\sum _{k=0}^{\infty }\frac{{\left(-t\right)}^k}{k!}{{\mathcal{L}}^{RW}}^k\mathbf{\theta }\left(0\right),}$ where $\mathbf{\theta }\left(0\right)$ is the state vector describing the initial distribution of the random walk. The stationary distribution of the random walk: (15)${\displaystyle \underset{t\to \infty }{lim}\mathbf{\theta }\left(t\right)=\mathbf{\pi }\text{with}{\pi }_i=\frac{{d_{\mathit{W}}}_i{\tau }_i}{\sum _j{d_{\mathit{W}}}_j{\tau }_j},}$ because ${\displaystyle \left({\mathit{D}}_{\mathit{W}}-\mathit{W}\right){\left(\mathit{T}{\mathit{D}}_{\mathit{W}}\right)}^{-1}\Pi =\left({\mathit{D}}_{\mathit{W}}-\mathit{W}\right)\vec{1}=\vec{0},}$ with $\mathbf{\pi }$ being the vector with π entries and Π being the diagonal matrix with the same elements. By convention, $\mathbf{\pi }$ is the standard centrality measure in conservative processes, including random walks (Ghosh & Lerman, 2012).

If we define centrality as the stationary distribution of a random walk, the importance of a vertex can be thought of as the total time a random walk spends at the vertex in the steady state. This is proportional to both vertex degree and delay factor, which we will later relate to the volume measure. If ${\mathcal{L}}^{RW}$ is a normalized Laplacian, this centrality measure is exactly the heat kernel page rank (Chung, 2009), which is identical to degree centralities since $\mathit{W}=\mathit{A}$ and $\mathit{T}=\mathit{I}$.

Stationary distribution of consensus dynamics

In consensus processes, the state vector always converges to a uniform state, where each vertex has the same value of the dynamic variable. As a result, the stationary distribution is not an appropriate measure of vertex centrality, since it deem all vertices to be equally important. However, the final consensus value associated with each vertex is (16)${\displaystyle {\pi }_i=\frac{1}{\sum _j{d_{\mathit{W}}}_j{\tau }_j}\sum _{i\in \mathit{V}}{\theta }_i\left(0\right){d_{\mathit{W}}}_i{\tau }_i,}$ where weight of vertex i in this average is $\frac{{d_{\mathit{W}}}_i{\tau }_i}{{\sum }_j{d_{\mathit{W}}}_j{\tau }_j}.$

Intuitively, as a measure of importance, it make sense to define the centrality of a vertex in the consensus process as its contribution to the final value. This consistency between “consensus” and “random walk” leads us to define the parameterized centrality.

Parameterized centrality

As shown in ‘Similarity transformations,’ the matrices connected through a similarity transformation represent the same linear operator up to a change of basis. For example, the relationship between “consensus” and “random walk” dynamics are captured by Fig. 2.

The above equivalence applies to all state vectors at any time t, including the stationary state. To verify, we first rewrite the initial state vector in terms of the eigenvectors of $\mathcal{L}$ $\left\{\overrightarrow{v_1},\overrightarrow{v_2},\dots ,\overrightarrow{v_n}\right\}$, indexed by their corresponding eigenvalues in ascending order λ₁ < λ₂ < ⋯, with the smallest λ₁ as the dominant eigenvalue. (17)${\displaystyle \mathbf{\theta }\left(t\right)=e^{-\mathcal{L}t}\cdot \mathbf{\theta }\left(0\right)=\sum _{k=0}^{\infty }\frac{{\left(-t\right)}^k}{k!}{\mathcal{L}}^k\mathbf{\theta }\left(0\right)=\sum _i\sum _{k=0}^{\infty }\frac{{\left(-t\right)}^k}{k!}{{\lambda }_i}^k{z}_i\overrightarrow{v_i}=\sum _i{z}_i{e}^{-{\lambda }_{i}t}\overrightarrow{v_i}=\sum _i{\mathit{u}}_i^T\mathbf{\theta }\left(0\right)e^{-{\lambda }_{i}t}\overrightarrow{v_i}=\mathbb{V}{e}^{-\Lambda t}{\mathbb{U}}^T\mathbf{\theta }\left(0\right),}$ where in the last step we used matrices to simplify the notation, with Λ being the diagonal matrix of eigenvalues, $\mathbb{V}$ composed of $\left\{\overrightarrow{v_1},\overrightarrow{v_2},\dots ,\overrightarrow{v_n}\right\}$ as columns and ${\mathbb{U}}^T={\mathbb{V}}^{-1}$. One interesting observation is that by left multiplying both sides with ${\mathbb{U}}^T$, we have ${\displaystyle {\mathbb{U}}^T\mathbf{\theta }\left(t\right)={\mathbb{U}}^T\mathbb{V}{e}^{-\Lambda t}{\mathbb{U}}^T\mathbf{\theta }\left(0\right)=e^{-\Lambda t}{\mathbb{U}}^T\mathbf{\theta }\left(0\right).}$ Recall that ${\mathbb{U}}^T\mathbf{\theta }$ is a vector in the eigenbasis $\mathbb{V}$. Applying the operator $\mathcal{L}$ to any input vector simply re-scales it according to eigenvalues. Since the smallest eigenvalue of the parameterized Laplacian is always 0, we have ${\displaystyle {\mathit{u}}_1^T\mathbf{\theta }\left(t\right)=e^{-{\lambda }_{1}t}{\mathit{u}}_1^T\mathbf{\theta }\left(0\right)={\mathit{u}}_1^T\mathbf{\theta }\left(0\right),}$ which states that the state vector is conserved along the direction of the dominant eigenvector $\overrightarrow{v_1}$.

The state vector reaches a stationary distribution π (18)${\displaystyle \pi =\underset{t\to \infty }{lim}\mathbf{\theta }\left(t\right)=\underset{t\to \infty }{lim}{e}^{{\lambda }_{1}t}\mathbf{\theta }\left(t\right)=z_1{\left(\frac{e^{{\lambda }_1}}{e^{{\lambda }_1}}\right)}^t\overrightarrow{v_1}+z_2{\left(\frac{e^{{\lambda }_1}}{e^{{\lambda }_2}}\right)}^t\overrightarrow{v_2}+\cdots +z_n{\left(\frac{e^{{\lambda }_1}}{e^{{\lambda }_n}}\right)}^t\overrightarrow{v_n}\approx z_1\overrightarrow{v_1}.}$ Since all terms vanish as t → ∞, the stationary state vector π only depends on $\overrightarrow{v_1}$. $z_1\overrightarrow{v_1}$ qualifies as a time invariant, initialization-independent vertex centrality measure.

Table 2 summarizes the properties of the stationary distributions and centralities associated with different similarity transformation of the parameterized Laplacian. ${\left[\mathbf{\theta }\right]}_{\rho }$ represents the vector $\mathbf{\theta }$ under the basis specified by the ρ parameter, with the random walk vector under the standard basis being $\mathbf{\theta }\left(0\right)$.

The spectral theorem states that any symmetric real matrix, regardless if its rank, has an orthonormal basis $\mathbb{V}$ which consists of its eigenvectors. Under the parameterized Laplacian framework, the symmetric formulation with ρ = 0 falls into this category. In the above table, we have chosen the normalization of the orthonormal basis $\sqrt{{\sum }_j{d}_j{\tau }_j}$ as the common normalization for all formulations.

As the table shows, similarity transformations of the same operator give the same the state vector $\mathbf{\theta }$, as long as the input and output vectors are properly transformed into the correct basis. They represent the same dynamics in different coordinate systems. Since centrality is determined by the dynamic process on a given network, it should be unified across these similarity transformations. In theory, any coordinate system can be set as the standard. Here, following the intuitions described earlier, we define the unnormalized stationary state vector of the random walk as the parameterized centrality: (19)${\displaystyle c_i={d_{\mathit{W}}}_i{\tau }_i.}$

Another motivation behind this definition is to establish a direct connection between centrality and community measures, as we will later demonstrate with the notion of parameterized volume (23).

Transformations and special cases

Parameterized centrality includes many well known centrality measures as special cases. Below, we summarize the induced special cases discussed in the previous subsection.

Parameterized conductance

Recall that conductance is a community quality measure associated with unbiased random walks. (20)${\displaystyle \varphi \left(S\right)=\frac{{\mathrm{cut}}_{\mathit{A}}\left(S,\bar{S}\right)}{min\left(\mathrm{vol}\left(S\right),\mathrm{vol}\left(\bar{S}\right)\right)},}$ where $\mathrm{vol}\left(S\right)={\sum }_{i\in S}{d}_i$ and ${\mathrm{cut}}_{\mathit{A}}={\sum }_{i\in S,j\in \bar{S}}{a}_{ij}$.

We generalize this notion with a claim that every dynamic process has an associated function that measures the quality of the cluster with respect to that process. Optimizing the quality function leads to cohesive communities, i.e., groups of vertices that “trap” the specific dynamic process for a long period of time.

Consider a dynamic process defined by the spreading operator $\mathcal{L}={\mathit{T}}^{-1∕2}{\mathit{D}}_{\mathit{W}}^{-1∕2}\left({\mathit{D}}_{\mathit{W}}-\mathit{W}\right){\mathit{D}}_{\mathit{W}}^{-1∕2}{\mathit{T}}^{-1∕2}$. We define the parameterized conductance of a set S with respect to $\mathcal{L}$ as: (21)${\displaystyle h_{\mathcal{L}}\left(S\right)=\frac{{\mathrm{cut}}_{\mathit{W}}\left(S,\bar{S}\right)}{min\left({\mathrm{vol}}_{\mathcal{L}}\left(S\right),{\mathrm{vol}}_{\mathcal{L}}\left(\bar{S}\right)\right)}=\frac{\sum _{i\in S,j\in \bar{S}}{w}_{ij}}{min\left(\sum _{i\in \mathrm{S}}{d_{\mathit{W}}}_i{\tau }_i,\sum _{i\in \bar{\mathrm{S}}}{d_{\mathit{W}}}_i{\tau }_i\right)}.}$ The minimum over all possible S is the parameterized conductance of the graph, (22)${\displaystyle {\varphi }_{\mathcal{L}}\left(G\right)=\underset{S\in V}{min}{h}_{\mathcal{L}}\left(S\right).}$ Notice that we have also defined the parameterized volume of a set S⊆V as (23)${\displaystyle {\mathrm{vol}}_{\mathcal{L}}\left(S\right)=\sum _{i\in S}{c}_i=\sum _{i\in S}{d_{\mathit{W}}}_i{\tau }_i,}$ which is the sum of parameterized centralities of member vertices. Using the random walk perspective, the numerator measures the random jumps across communities, while the denominator ensures a balanced bisection. As previously pointed out, the presence of a good cut implies that it will take a random walk a long time to cross this boundary and reach its stationary distribution. This corresponds to a small numerator. The parameterized volume can be interpreted as the total time a random walk stays within a community after convergence, as it is proportional to both vertex degrees and vertex delay factors. This interpretation of the denominator coincides with our definition of parameterized centrality (19).

Transformations and special cases

We can use any transformation to produce new dynamics, and the corresponding parameterized conductance will be redefined according to Eq. (21), (24)${\displaystyle h_{\mathcal{L}}\left(S\right)=\frac{\sum _{i\in S,j\in \bar{S}}{w}_{ij}}{min\left(\sum _{i\in \mathrm{S}}{d_{\mathit{W}}}_i{\tau }_i,\sum _{i\in \bar{\mathrm{S}}}{d_{\mathit{W}}}_i{\tau }_i\right)}.}$

However, the effect of transformations on the resulting communities is not as obvious when compared with the parameterized centrality. Below, we elaborate the effect of transformations on the parameterized conductance measure in cases and examples.

First of all, the similarity transformation keeps both numerator and denominator the same, which makes the quality function of the same communities identical. This ultimately leads to identical parameterized conductances, which is the minimum over all possible bisections. Uniform scaling does change the denominator. However, because all possible bisections are scaled uniformly, the relative quality measure remain the same, leading to identical parameterized conductances communities.

From the algorithmic perspective, both similarity and uniform scaling transformations preserve spectral properties of the operator. Since the spectrum is the only input information our spectral dynamics clustering Algorithm 1 uses, we always expect to get the same solution after the transformations. This is not the case with non-uniform scaling and reweighing transformations.

With non-uniform scaling, the numerator remains unchanged. It is each vertex’s delay time change that scales the volume measures in the denominator, which in turn results in different optimal bisections because of the balance constraint.

The reweighing transformation is the most complex of all, changing both the numerator and denominator in Eq. (21). This trade-off between cut and balance can oftentimes be very complicated to analyze (as will be seen with real world networks).

Finally, we summarize the induced special cases.

Parameterized Cheeger inequality

Given the parameterized conductance measure, finding the best community bisection is still a combinatorial problem, which quickly becomes computationally intractable as the network grows in size. In this subsection we will extend the theorems for the classic Laplacian to our parameterized setting, ultimately leading to efficient approximate algorithms with theoretical guarantees. For mathematical convenience we will use the symmetric formulation and assume that ρ = 0 for $\mathcal{L}$. Cheeger inequality (Cheeger, 1970) states that ${\displaystyle {\varphi }^2\left(G\right)∕2\le {\lambda }_2\le 2\varphi \left(G\right)}$ where λ₂ is the second smallest eigenvalue of the symmetric normalized Laplacian, $\mathcal{L}=\mathit{I}-{\mathit{D}}^{-1∕2}{\mathit{WD}}^{-1∕2},$ and ϕ(G) is conductance. The relationship between conductance and spectral properties of the Laplacian enables the use of its eigenvectors for partitioning graphs, particularly the nearest-neighbor graphs and finite-element meshes (Spielman & Teng, 1996).

In this section, we generalize Cheeger inequality to any spreading operator under our framework and its associated parameterized conductance of the graph (given by Eq. (22)). Compared with classic results in Markov chain mixing times (Jerrum & Sinclair, 1988; Lawler & Sokal, 1988), we generalize Cheeger inequality to accommodate the asyncronized delay factors in $\mathit{T}$. It also comes with algorithmic consequences, leading to spectral partitioning algorithms that are efficient in finding low conductance cuts for a given operator.