Transfer operators on graphs: spectral clustering and beyond

In what follows, we will mainly consider directed graphs, but the derived methods can, of course, also be applied to undirected graphs. Relationships with well-known clustering techniques for undirected graphs will be discussed below.

Definition 2.1 (directed graph). A directed graph $\mathscr{G} = (\mathscr{V}, \mathscr{E})$ is given by a set of vertices $\mathscr{V} = \{{\mathscr{v}_1}, \dots, {\mathscr{v}_n} \}$ and edges $\mathscr{E} \subseteq \mathscr{V} \times \mathscr{V}$.

The weighted adjacency matrix $A \in \mathbb{R}^{n \times n}$ associated with a graph $\mathscr{G}$ is defined by

$\begin{align*} a_{ij} = \begin{cases} w\left({\mathscr{v}_i}, {\mathscr{v}_j}\right), & \textrm{if } \left({\mathscr{v}_i}, {\mathscr{v}_j}\right) \in \mathscr{E}, \\ 0, & \textrm{otherwise}, \end{cases} \end{align*}$

where $w({\mathscr{v}_i}, {\mathscr{v}_j}) \gt 0$ is the weight of the corresponding edge. If the graph is undirected, then the adjacency matrix is symmetric. Additionally, we introduce the row-stochastic transition probability matrix $S = D_{\mathscr{o}}^{-1} A \in \mathbb{R}^{n \times n}$, with

$\begin{align*} D_{\mathscr{o}} = \operatorname{diag}\left(\mathscr{o}\left({\mathscr{v}_1}\right), \dots, \mathscr{o}\left({\mathscr{v}_n}\right)\right) \quad \textrm{and} \quad \mathscr{o}\left({\mathscr{v}_i}\right) = \sum_{j = 1}^n a_{ij}. \end{align*}$

That is, s_ij is the probability that a random walker starting in ${\mathscr{v}_i}$ will go to ${\mathscr{v}_j}$ in one step. We refer to $D_{\mathscr{o}}$ as the matrix of out-degrees, and we will later also make use of the in-degree matrix $D_{\mathscr{i}}$, which is defined in a similar way, with

$\begin{align*} D_{\mathscr{i}} = \operatorname{diag}\left(\mathscr{i}\left({\mathscr{v}_1}\right), \dots, \mathscr{i}\left({\mathscr{v}_n}\right)\right) \quad \textrm{and} \quad \mathscr{i}\left({\mathscr{v}_i}\right) = \sum_{j = 1}^n a_{ji}. \end{align*}$

In order to generate benchmark problems with clearly defined clusters, we will use the so-called directed stochastic block model.

Definition 2.2 (directed stochastic block model). The graph $\mathscr{G}$ is said to be sampled from the directed stochastic block model (DSBM) $\mathbf{G}(r_b, n_b, E)$, with $E \in [0, 1]^{r_b \times r_b}$, if the adjacency matrix

$\begin{align*} A = \begin{bmatrix} A_{11} & \dots & A_{1 r_b} \\ \vdots & \ddots & \vdots \\ A_{r_b 1} & \dots & A_{r_b r_b} \end{bmatrix} \in \mathbb{R}^{\left(r_b \hspace{0.1em} n_b\right) \times \left(r_b \hspace{0.1em} n_b\right)} \end{align*}$

is a block matrix whose blocks $A_{ij} \in R^{n_b \times n_b}$ contain positive entries with probability F_ij . Here, r_b is the number and n_b the size of the blocks.

Our definition differs slightly from the DSBM described in [31], in that we allow different probabilities for all the blocks. Undirected graphs can be constructed in a similar way. While the DSBM is a frequently used model to generate benchmark problems, real-world graphs typically exhibit a more complicated structure. An algorithm to generate weighted directed graphs with heterogeneous node-degree distributions and cluster sizes for testing community detection methods is described in [32]. Briefly, the graphs in [32] are constructed by randomly assigning in-degrees to each node from a power law distribution and assigning out-degrees from a delta distribution. A mixing parameter is introduced for each of these degrees that is related to the quality of the graph partition—that is, the lower the mixing parameter the better the partitioning of the graph. Cluster sizes are also drawn from a power law distribution and a subgraph is constructed from each of these clusters. The subgraphs are finally joined randomly, with a rewiring process to ensure that the in-degree and out-degree distributions are preserved. We will use this algorithm to generate a set of graphs with different characteristics.

The spectral clustering approach for directed and time-evolving graphs proposed in [23] is based on a generalized Laplacian that contains information about coherent sets. The detection of such coherent sets requires the notion of transfer operators, which describe the evolution of probability densities or observables. Before we introduce these operators, let us first illustrate the definition of coherence—defined in [8, 33] for continuous dynamical systems—for directed graphs.

Definition 2.3 (coherent pair). Let S^τ be the flow associated with a dynamical system and τ a fixed lag time. Two sets $\mathbb{A}$ and $\mathbb{B}$ form a so-called coherent pair if $S^\tau(\mathbb{A}) \approx \mathbb{B}$ and $S^{-\tau}(\mathbb{B}) \approx \mathbb{A}$.

That is, the set $\mathbb{A}$ is almost invariant under the forward–backward dynamics and we call it a finite-time coherent set. To avoid that $(S^{-\tau} \circ S^\tau)(\mathbb{A}) = \mathbb{A}$ for all sets $\mathbb{A}$, a small random perturbation of the dynamics is required for deterministic systems, see [8]. In our setting, the dynamics are given by a random walk process on the graph and thus already non-deterministic.

Example 2.4. Figure 1 illustrates the difference between a coherent set (green) and a set that is dispersed by the random-walk dynamics (red). The example shows that coherent sets form weakly coupled clusters and can be regarded as a natural generalization of clusters in undirected graphs. $\triangle$

Zoom In Zoom Out Reset image size

This notion of coherence can be easily extended to time-evolving networks as shown in [23]. The structure of the clusters may then also change in time.

Transfer operators on graphs can be either directly expressed in terms of graph-related matrices (such as the transition probability matrix) or estimated from random walk data. We will start with a formal definition of transfer operators. In what follows, $\mathbb{U} = \{f \colon \mathcal{V} \to \mathbb{R} \}$ denotes the set of all real-valued functions defined on the set of vertices $\mathcal{V}$ of a graph $\mathscr{G}$.

Definition 2.5 (Perron–Frobenius & Koopman operators). Let S be the transition probability matrix associated with the graph $\mathscr{G}$.

• (i)  
  For a function $\rho \in \mathbb{U}$, the Perron–Frobenius operator $\mathcal{P} \colon \mathbb{U} \to \mathbb{U}$ is defined by

  $\begin{align*} \mathcal{P} \rho\left({\mathscr{v}_i}\right) = \sum_{j = 1}^n s_{ji} \hspace{0.1em} \rho\left({\mathscr{v}_j}\right). \end{align*}$
• (ii)  
  Similarly, for a function $f \in \mathbb{U}$, the Koopman operator $\mathcal{K} \colon \mathbb{U} \to \mathbb{U}$ is given by

  $\begin{align*} \mathcal{K} f\left({\mathscr{v}_i}\right) = \sum_{j = 1}^n s_{ij} \hspace{0.1em} f\left({\mathscr{v}_j}\right). \end{align*}$

The Koopman operator is the adjoint of the Perron–Frobenius operator with respect to the standard inner product.

Definition 2.6 (probability density). A probability density is a function $\mu \in \mathbb{U}$ that satisfies $\mu({\mathscr{v}_i}) \unicode{x2A7E} 0$ and

$\begin{align*} \sum_{i = 1}^n \mu\left({\mathscr{v}_i}\right) = 1. \end{align*}$

This allows us to define the Perron–Frobenius operator with respect to different probability densities. We define the µ-weighted inner product by

$\begin{align*} \left\langle \,f,\, g \right\rangle_{\mu} = \sum_{i = 1}^n \mu\left({\mathscr{v}_i}\right) \hspace{0.1em} f\left({\mathscr{v}_i}\right) \hspace{0.1em} g\left({\mathscr{v}_i}\right). \end{align*}$

Definition 2.7 (reversibility). The system is called reversible if there exists a probability density π that satisfies the detailed balance condition $\pi({\mathscr{v}_i}) \hspace{0.1em} s_{ij} = \pi({\mathscr{v}_j}) \hspace{0.1em} s_{ji}$ for all ${\mathscr{v}_i}, {\mathscr{v}_j} \in \mathscr{V}$.

Random walks on undirected graphs, for instance, are reversible [34, 35]. We will use this property to analyze the relationships between conventional spectral clustering and our approach.

Definition 2.8 (reweighted Perron–Frobenius operator). Let µ be the initial density and $\nu = \mathcal{P} \mu$ the resulting image density, i.e.

$\begin{align*} \nu\left({\mathscr{v}_i}\right) = \sum_{j = 1}^n s_{ji} \hspace{0.1em} \mu\left({\mathscr{v}_j}\right). \end{align*}$

Furthermore, assume that µ and ν are strictly positive. Given a function $u \in \mathbb{U}$, we define the reweighted Perron–Frobenius operator $\mathcal{T} \colon \mathbb{U} \to \mathbb{U}$ by

$\begin{align*} \mathcal{T} u\left({\mathscr{v}_i}\right) = \frac{1}{\nu\left({\mathscr{v}_i}\right)} \sum_{j = 1}^n s_{ji} \hspace{0.1em} \mu\left({\mathscr{v}_j}\right) \hspace{0.1em} u\left({\mathscr{v}_j}\right). \end{align*}$

The reweighted Perron–Frobenius operator describes the evolution of functions with respect to the initial and final densities µ and ν. The assumption that $\mu({\mathscr{v}_i}) \gt 0$ implies that the probability that a random walker starts in ${\mathscr{v}_i}$ is non-zero. Correspondingly, $\nu({\mathscr{v}_i}) \gt 0$ means that ${\mathscr{v}_i}$ is reachable from at least one vertex, which could be ${\mathscr{v}_i}$ itself. Adding self-loops thus regularizes the problem.

Definition 2.9 (forward–backward & backward–forward operators). Let $f, u \in \mathbb{U}$ as before.

• (i)  
  We define the forward–backward operator $\mathcal{F} \colon \mathbb{U} \to \mathbb{U}$ by

  $\begin{align*} \mathcal{F} u\left({\mathscr{v}_i}\right) = \mathcal{K} \mathcal{T} u\left({\mathscr{v}_i}\right) = \sum_{j = 1}^n s_{ij} \frac{1}{\nu\left({\mathscr{v}_j}\right)} \sum_{k = 1}^n s_{kj} \hspace{0.1em} \mu\left({\mathscr{v}_{k}}\right) \hspace{0.1em} u\left({\mathscr{v}_{k}}\right). \end{align*}$
• (ii)  
  Analogously, the backward–forward operator $\mathcal{B} \colon \mathbb{U} \to \mathbb{U}$ is defined by

  $\begin{align*} \mathcal{B} f\left({\mathscr{v}_i}\right) = \mathcal{T} \mathcal{K} f\left({\mathscr{v}_i}\right) = \frac{1}{\nu\left({\mathscr{v}_i}\right)} \sum_{j = 1}^n s_{ji} \hspace{0.1em} \mu\left({\mathscr{v}_j}\right) \sum_{k = 1}^n s_{jk} \hspace{0.1em} f\left({\mathscr{v}_{k}}\right). \end{align*}$

In order to detect coherent sets, we want to find eigenfunctions $\varphi_{\ell}$ of the forward–backward operator $\mathcal{F}$ corresponding to eigenvalues $\lambda_{\ell} \approx 1$ such that $\mathcal{F} \varphi_{\ell} = \lambda_{\ell} \hspace{0.1em} \varphi_{\ell} \approx \varphi_{\ell}$. These functions are almost invariant under the forward–backward dynamics. Defining $\psi_{\ell} : = \frac{1}{\sqrt{\lambda_{\ell}}} \mathcal{T} \varphi_{\ell}$, this results in

$\begin{align*} \mathcal{K} \psi_{\ell} = \sqrt{\lambda_{\ell}} \hspace{0.1em} \varphi_{\ell}, \quad \mathcal{T} \varphi_{\ell} = \sqrt{\lambda_{\ell}} \hspace{0.1em} \psi_{\ell}, \quad \mathcal{F} \varphi_{\ell} = \lambda_{\ell} \hspace{0.1em} \varphi_{\ell}, \quad \mathcal{B} \psi_{\ell} = \lambda_{\ell} \hspace{0.1em} \psi_{\ell}. \end{align*}$

That is, the functions $\varphi_{\ell}$ and $\psi_{\ell}$ come in pairs.

In addition to these different types of transfer operators, we define (uncentered) covariance and cross-covariance operators.

Definition 2.10 (covariance operators). Let S be the transition probability matrix and µ and ν the densities defined above. Given functions $f, g \in \mathbb{U}$, we call $\mathcal{C}_{xx}, \mathcal{C}_{yy} \colon \mathbb{U} \to \mathbb{U}$, with

$\begin{align*} \mathcal{C}_{xx} \hspace{0.1em} f\left({\mathscr{v}_i}\right) = \mu\left({\mathscr{v}_i}\right) \hspace{0.1em} f\left({\mathscr{v}_i}\right) \quad \textrm{and} \quad \mathcal{C}_{yy} \hspace{0.1em} g\left({\mathscr{v}_i}\right) = \nu\left({\mathscr{v}_i}\right) \hspace{0.1em} g\left({\mathscr{v}_i}\right), \end{align*}$

covariance operators and $\mathcal{C}_{xy} \colon \mathbb{U} \to \mathbb{U}$, with

$\begin{align*} \mathcal{C}_{xy} \hspace{0.1em} g\left({\mathscr{v}_i}\right) = \sum_{j = 1}^n \mu\left({\mathscr{v}_i}\right) \hspace{0.1em} s_{ij} \hspace{0.1em} g\left({\mathscr{v}_j}\right), \end{align*}$

cross-covariance operator.

Let $X \sim \mu$ and $Y \sim \nu$ be random variables representing the initial position of the random walker and the position after one step, respectively. It follows that

$\begin{align*} \left\langle \,f,\, \mathcal{C}_{xx} \hspace{0.1em} f\, \right\rangle & = \left\langle \,f,\, f\, \right\rangle_{\mu} = \sum_{i = 1}^n \mu\left({\mathscr{v}_i}\right) \hspace{0.1em} f\left({\mathscr{v}_i}\right)^2 = \mathbb{E}_X\left[f\left(X\right)^2\right], \\ \left\langle g,\, \mathcal{C}_{yy} \hspace{0.1em} g \right\rangle & = \left\langle g,\, g \right\rangle_{\nu} = \sum_{i = 1}^n \nu\left({\mathscr{v}_i}\right) \hspace{0.1em} g\left({\mathscr{v}_i}\right)^2 = \mathbb{E}_Y\left[g\left(Y\right)^2\right], \end{align*}$

and, since the joint probability distribution is given by $\mathbb{P}[X = {\mathscr{v}_i}, Y = {\mathscr{v}_j}] = \mu({\mathscr{v}_i}) \hspace{0.1em} s_{ij}$,

$\begin{align*} \left\langle \,f,\, \mathcal{C}_{xy} \hspace{0.1em} g \right\rangle = \left\langle \,f,\, \mathcal{K} g \right\rangle_{\mu} = \sum_{i = 1}^n \sum_{j = 1}^n \mu\left({\mathscr{v}_i}\right) \hspace{0.1em} s_{ij} \hspace{0.1em} f\left({\mathscr{v}_i}\right) \hspace{0.1em} g\left({\mathscr{v}_j}\right) = \mathbb{E}_{XY}\left[f\left(X\right) \hspace{0.1em} g\left(Y\right)\right]. \end{align*}$

In order to simplify the notation, given any function $g \in \mathbb{U}$, let $\boldsymbol{g} \in \mathbb{R}^n$ be defined by

$\begin{align*} \boldsymbol{g} = \left[\, g\left({\mathscr{v}_1}\right), \dots, g\left({\mathscr{v}_n}\right)\, \right]^\top, \end{align*}$

i.e. $\boldsymbol{\rho}, \boldsymbol{f}, \boldsymbol{u}, \boldsymbol{\mu}, \boldsymbol{\nu} \in \mathbb{R}^n$ are the vector representations of the functions $\rho, f, u, \mu, \nu \in \mathbb{U}$, respectively. Additionally, we define diagonal matrices

$\begin{align*} D_{\mu} = \operatorname{diag}\left(\boldsymbol{\mu}\right) \quad \textrm{and} \quad D_{\nu} = \operatorname{diag}\left(\boldsymbol{\nu}\right) \end{align*}$

so that we can write

$\begin{align*} \mathcal{P} \boldsymbol{\rho} &: = P \hspace{0.1em} \boldsymbol{\rho} = S^\top \boldsymbol{\rho}, \\ \mathcal{K} \boldsymbol{f} &: = K \hspace{0.1em} \boldsymbol{f} = S \boldsymbol{f}, \\ \mathcal{T} \hspace{0.1em} \boldsymbol{u} &: = T \hspace{0.1em} \boldsymbol{u} = D_{\nu}^{-1} \hspace{0.1em} S^\top D_{\mu} \hspace{0.1em} \boldsymbol{u}, \\ \mathcal{F} \hspace{0.1em} \boldsymbol{u} &: = F \hspace{0.1em} \boldsymbol{u} = S \hspace{0.1em} D_{\nu}^{-1} \hspace{0.1em} S^\top D_{\mu} \hspace{0.1em} \boldsymbol{u}, \\ \mathcal{B} \hspace{0.1em} \boldsymbol{f} &: = B \hspace{0.1em} \boldsymbol{f} = D_{\nu}^{-1} \hspace{0.1em} S^\top D_{\mu} \hspace{0.1em} S \hspace{0.1em} \boldsymbol{f}, \end{align*}$

where the operators (with a slight abuse of notation) are applied component-wise. That is, $P, K, T, F, B$ are the matrix representations of the corresponding operators $\mathcal{P}, \mathcal{K}, \mathcal{T}, \mathcal{F}, \mathcal{B}$, respectively. The matrices D_µ and D_ν are invertible since we assumed µ and ν to be strictly positive. For the covariance and cross-covariance operators $\mathcal{C}_{xx}$, $\mathcal{C}_{yy}$, and $\mathcal{C}_{xy}$, we obtain the matrix representations $C_{xx} = D_{\mu}$, $C_{yy} = D_{\nu}$, and $C_{xy} = D_{\mu} \hspace{0.1em} S$ and can thus write

$\begin{align*} \left\langle \,f,\, \mathcal{C}_{xx} \hspace{0.1em} f\, \right\rangle = \mathbf{f}^\top D_{\mu} \hspace{0.1em} \mathbf{f}, \quad \left\langle g,\, \mathcal{C}_{yy} \hspace{0.1em} g \right\rangle = \mathbf{g}^\top D_{\nu} \hspace{0.1em} \mathbf{g}, \quad \textrm{and} \quad \left\langle \,f,\, \mathcal{C}_{xy} \hspace{0.1em} g \right\rangle = \mathbf{f}^\top D_{\mu} S \hspace{0.1em} \mathbf{g}. \end{align*}$

In what follows, we will use functions and their vector representations as well as operators and their matrix representations interchangeably.

Remark 2.11. Let $\unicode{x1D7D9} \in \mathbb{R}^n$ be the vector of all ones. In [23], we assumed that the random walkers are initially uniformly distributed and defined $\widetilde{F} = S \hspace{0.1em} \widetilde{D}_{\nu}^{-1} S^\top$, with $\widetilde{D}_{\nu} = \operatorname{diag}(S^\top \unicode{x1D7D9})$, whereas setting $\boldsymbol{\mu} = \frac{1}{n} \unicode{x1D7D9}$ results in $D_{\nu} = \frac{1}{n}\operatorname{diag}(S^\top \unicode{x1D7D9})$. However, for this special case $F = S \hspace{0.1em} D_{\nu}^{-1} S^\top D_{\mu} = S \hspace{0.1em} \widetilde{D}_{\nu}^{-1} S^\top = \widetilde{F}$ and the definitions are equivalent.

The properties of transfer operators are analyzed in detail in appendix. We in particular show that the eigenvalues of the operators $\mathcal{F}$ and $\mathcal{B}$ are contained in the interval $[0, 1]$.

How are the graph structures associated with the matrix representations F and B of the forward–backward and backward–forward operators $\mathcal{F}$ and $\mathcal{B}$ related to the original directed graph $\mathscr{G}$?

Lemma 2.12. It holds that the entry f_ij of the matrix F is nonzero if and only if there exists an index k such that $({\mathscr{v}_i}, {\mathscr{v}_{k}}) \in \mathscr{E}$ and $({\mathscr{v}_j}, {\mathscr{v}_{k}}) \in \mathscr{E}$. That is, a random walker can go forward from ${\mathscr{v}_i}$ to ${\mathscr{v}_{k}}$ and then backward to ${\mathscr{v}_j}$. Similarly, the entry b_ij of the matrix B is nonzero if and only if there exists an index k such that $({\mathscr{v}_{k}}, {\mathscr{v}_i}) \in \mathscr{E}$ and $({\mathscr{v}_{k}}, {\mathscr{v}_j}) \in \mathscr{E}$.

Proof. Since $s_{ij} \neq 0 ~\Longleftrightarrow~ ({\mathscr{v}_i}, {\mathscr{v}_j}) \in \mathscr{E}$, we have

$\begin{align*} f_{ij} = \sum_{k = 1}^n \frac{s_{ik} \hspace{0.1em} s_{jk} \hspace{0.1em} \mu\left({\mathscr{v}_j}\right)}{\nu\left({\mathscr{v}_{k}}\right)} \neq 0 ~\Longleftrightarrow~ \exists k: \left({\mathscr{v}_i}, {\mathscr{v}_{k}}\right) \in \mathscr{E} \land \left({\mathscr{v}_j}, {\mathscr{v}_{k}}\right) \in \mathscr{E} \end{align*}$

and

$\begin{align*} b_{ij} = \sum_{k = 1}^n \frac{s_{ki} \hspace{0.1em} s_{kj} \hspace{0.1em} \mu\left({\mathscr{v}_{k}}\right)}{\nu\left({\mathscr{v}_i}\right)} \neq 0 ~\Longleftrightarrow~ \exists k: \left({\mathscr{v}_{k}}, {\mathscr{v}_i}\right) \in \mathscr{E} \land \left({\mathscr{v}_{k}}, {\mathscr{v}_j}\right) \in \mathscr{E}. \end{align*}$

 □

The difference between the two matrix representations of the operators is illustrated in figure 2. Note that the probability of each forward–backward random walk via ${\mathscr{v}_{k}}$ is divided by $\nu({\mathscr{v}_{k}})$, which represents the sum of the probabilities to go from any vertex to ${\mathscr{v}_{k}}$. This is intuitively clear since the more edges enter ${\mathscr{v}_{k}}$, the more options the random walker has to go backward. On the other hand, the more forward–backward random walks from ${\mathscr{v}_i}$ to ${\mathscr{v}_j}$ exist, the larger the weight f_ij .

Remark 2.13. This is strongly related to the bibliometric graph clustering approach described in [22]. Here, the bibliographic coupling matrix $A A^\top$ and the co-citation matrix $A^\top A$ represent the common out-links and in-links between a pair of nodes, respectively. Using the sum of these two matrices gives a symmetrization of the adjacency matrix that accounts for both of these types of common links. As in the forward–backward approach, the weights can be divided by the degrees of the nodes in order to give an effective measure of similarity for nodes that accounts for the number of possible routes or links into and out of the node. Alternatively, in [31] a complex-valued Hermitian matrix $H = i(A-A^\top)$ is constructed so that $H^2 = A A^\top + A^\top A - A^2 - (A^\top)^2$. This can then also be viewed as a symmetrization scheme that utilizes the bibliographic coupling matrix $A A^\top$ and the co-citation matrix $A^\top A$.

Example 2.14. Let us consider the graph depicted in figure 3(a), which we already analyzed in example 2.4, and compute the associated matrices F and B. The resulting graph structures are shown in figures 3(b) and (c). The graph corresponding to the matrix F, for instance, now contains an edge connecting ${\mathscr{v}_4}$ and ${\mathscr{v}_{8}}$ since a random walker could move forward from ${\mathscr{v}_4}$ to ${\mathscr{v}_{5}}$ and backward to ${\mathscr{v}_{8}}$ or the other way around. The probability of this forward–backward random walk, however, is low due to the small edge weight of $({\mathscr{v}_4}, {\mathscr{v}_{5}})$. Note that without the self-loops all vertices except ${\mathscr{v}_4}$, ${\mathscr{v}_{8}}$, and ${\mathscr{v}_{12}}$ would become disconnected. Every forward–backward random walk starting in one of these vertices would then end up in the initial position (since these vertices have only one outgoing edge) and thus form clusters of size one. This illustrates that adding self-loops can be regarded as a form of regularization. $\triangle$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Spectral clustering methods for undirected graphs are based on the eigenvectors corresponding to the smallest eigenvalues of an associated graph Laplacian. The random-walk normalized Laplacian is defined by

$\begin{align*} L_{\mathcal{K}} = I - D_{\mathscr{o}}^{-1} A = I - K, \end{align*}$

where K is the matrix representation of the Koopman operator. Equivalently, we can also compute the eigenvectors corresponding to the largest eigenvalues of K. This shows that conventional spectral clustering computes the dominant eigenfunctions of the Koopman operator, which contain information about the metastable sets of the system. The eigenvalues of the Koopman operator associated with directed graphs, however, are in general complex-valued and standard spectral clustering techniques typically fail, see [23] for a detailed analysis. For directed (and also time-evolving) graphs we thus propose to compute the eigenvectors associated with the largest eigenvalues of the forward–backward or backward–forward operators in order to detect coherent sets. The eigenvalues of these operators are real-valued even if the graph is directed, see appendix.

Algorithm 2.15 (transfer operator-based spectral clustering algorithm). 

• 1.  
  Compute the k largest eigenvalues $\lambda_{\ell}$ and associated eigenvectors $\boldsymbol{\varphi}_{\ell}$ of F.
• 2.  
  Define $\boldsymbol{\Phi} = [\boldsymbol{\varphi}_1, \dots, \boldsymbol{\varphi}_k] \in \mathbb{R}^{n \times k}$ and let r_i denote the ith row of Φ.
• 3.  
  Cluster the points $\{\hspace{0.1em} r_i \hspace{0.1em} \}_{i = 1}^n$ using, e.g. k-means.

The number of eigenvalues close to 1 indicates the number of coherent sets. That is, k should be chosen in such a way that there is a spectral gap between λ_k and $\lambda_{k+1}$. In order to illustrate the spectral clustering algorithm, we first consider a simple graph that has a clearly defined cluster structure.

Example 2.16. We sample a benchmark problem $\mathscr{G} \in \mathbf{G}(4, 100, E)$ from the DSBM, with

$\begin{align*} E = \begin{bmatrix} q & p & q & q \\ q & q & q & p \\ q & q & p & q \\ p & q & q & q \end{bmatrix}, \end{align*}$

where p = 0.8 and q = 0.1. The resulting adjacency matrix is shown in figure 4(a). Note that there are two different types of clusters, namely dense blocks on the diagonal, representing groups of vertices that are tightly coupled to other vertices contained in the same cluster, and dense off-diagonal blocks, which have the property that the contained vertices are all tightly coupled (unidirectionally) to vertices in another cluster. For both types of clusters, however, the probability that a random walker going forward and then backward ends up in the same cluster is large compared to the escape probability. Applying algorithm 2.15 yields the results shown in figures 4(c) and (d). The dominant eigenvalues are $\lambda_1 = 1$, $\lambda_2 = 0.72$, $\lambda_3 = 0.70$, and $\lambda_4 = 0.69$, followed by a spectral gap. The eigenfunctions of the forward–backward operator can also be used for visualizing directed graphs ³ as shown in figure 4(b). By defining the position of the vertex ${\mathscr{v}_i}$ to be $(\varphi_2({\mathscr{v}_i}), \varphi_3({\mathscr{v}_i}))$, we obtain a two-dimensional embedding of the graph. The edges are colored according to the source nodes. The blue vertices form a cluster in the conventional sense since they are strongly connected but only weakly coupled to the other clusters. The remaining sets of vertices have the property that many directed edges are pointing to one of the other clusters but not vice versa. These clusters correspond to non-sparse off-diagonal blocks. $\triangle$

Zoom In Zoom Out Reset image size

A natural question now is whether conventional spectral clustering methods for undirected graphs can be regarded as a special case of the algorithm derived above.

Lemma 2.17. Assume that the system is reversible with respect to π. Then, setting $\mu = \pi$, it holds that $\mathcal{T} = \mathcal{K}$ and $\mathcal{F} = \mathcal{B} = \mathcal{K}^2$.

Proof. First, we rewrite the detailed balance condition as $D_{\pi} \hspace{0.1em} S = S^\top D_{\pi}$. Since $\mu = \pi$ implies $\nu = \pi$, we obtain $T = D_{\pi}^{-1} \hspace{0.1em} S^\top D_{\pi} = S = K$. The second part follows immediately from the definitions of the operators $\mathcal{F}$ and $\mathcal{B}$. □

If π satisfies the detailed balance condition, then it is also automatically an invariant density. This can be seen as follows:

$\begin{align*} \mathcal{P} \boldsymbol{\pi} = S^\top \hspace{0.1em} \boldsymbol{\pi} = S^\top D_{\pi} \hspace{0.1em} \unicode{x1D7D9} = D_{\pi} \hspace{0.1em} S \hspace{0.1em} \unicode{x1D7D9} = D_{\pi} \hspace{0.1em} \unicode{x1D7D9} = \boldsymbol{\pi}. \end{align*}$

Lemma 2.18. Let $\mathscr{G}$ now be an undirected graph with unique invariant density π. Assume that K has only distinct positive eigenvalues. Then conventional spectral clustering and algorithm 2.15 with $\mu = \pi$ compute the same clusters.

Proof. Let $\kappa_{\ell}$ be the eigenvalues of K and $\lambda_{\ell}$ the eigenvalues of F. Using lemma 2.17, we know that $\lambda_{\ell} = \kappa_{\ell}^2$ and the eigenvectors are identical. Since all eigenvalues of K are by assumption positive and distinct, the ordering of the eigenvalues and eigenvectors remains unchanged. □

Negative eigenvalues of K might affect the ordering of the eigenvectors and thus the spectral clustering results. However, positive eigenvalues can be enforced by turning the transition probability matrix into a lazy Markov chain [37].

Given an operator $\mathcal{L} \colon \mathbb{U} \to \mathbb{U}$ with matrix representation L, our goal now is to compute its Galerkin projection $\mathcal{L}_r$ onto the r-dimensional subspace $\mathbb{U}_r \subset \mathbb{U}$ with basis $\{\phi_i \}_{i = 1}^r$. Here, r could potentially be much smaller than n. The matrix representation $L_r \in \mathbb{R}^{r \times r}$ of $\mathcal{L}_r$ is given by $L_r = \big(G^{(0)}\big)^{-1} G^{(1)}$, with entries

$\begin{align*} g_{ij}^{\left(0\right)} = \left\langle \phi_i,\, \phi_j \right\rangle_{\mu} \quad \textrm{and} \quad g_{ij}^{\left(1\right)} = \left\langle \phi_i,\, \mathcal{L} \phi_j \right\rangle_{\mu}. \end{align*}$

By defining the vector-valued function $\phi({\mathscr{v}_i}) = [\phi_1({\mathscr{v}_i}), \dots, \phi_r({\mathscr{v}_i})]^\top \in \mathbb{R}^r$ and

$\begin{align*} \Phi_{\mathscr{V}} = \begin{bmatrix} \phi\left({\mathscr{v}_1}\right), \dots, \phi\left({\mathscr{v}_n}\right) \end{bmatrix} \in \mathbb{R}^{r \times n}, \end{align*}$

the matrices $G^{(0)}$ and $G^{(1)}$ can be compactly written as

$\begin{align*} G^{\left(0\right)} & = \sum_{i = 1}^n \mu\left({\mathscr{v}_i}\right) \hspace{0.1em} \phi\left({\mathscr{v}_i}\right) \hspace{0.1em} \phi\left({\mathscr{v}_i}\right)^\top = \Phi_{\mathscr{V}} D_{\mu} \Phi_{\mathscr{V}}^\top \end{align*}$

and

$\begin{align*} G^{\left(1\right)} & = \sum_{i = 1}^n \sum_{j = 1}^n \mu\left({\mathscr{v}_i}\right) \hspace{0.1em} l_{ij} \hspace{0.1em} \phi\left({\mathscr{v}_i}\right) \hspace{0.1em} \phi\left({\mathscr{v}_j}\right)^\top = \Phi_{\mathscr{V}} D_{\mu} \hspace{0.1em} L \hspace{0.1em} \Phi_{\mathscr{V}}^\top. \end{align*}$

Any function $f \in \mathbb{U}_r$ is given by

$\begin{align*} f\left({\mathscr{v}_i}\right) = \sum_{j = 1}^r c_j \hspace{0.1em} \phi_j\left({\mathscr{v}_i}\right) = c^\top \phi\left({\mathscr{v}_i}\right), \end{align*}$

where $c = [c_1, \dots, c_r]^\top \in \mathbb{R}^r$, such that

$\begin{align*} \mathcal{L}_r f\left({\mathscr{v}_i}\right) = \left(L_r \hspace{0.1em} c\right)^\top \phi\left({\mathscr{v}_i}\right). \end{align*}$

Suppose now $\xi_{\ell}$ is an eigenvector of L_r with associated eigenvalue $\lambda_{\ell}$, i.e. $L_r \hspace{0.1em} \xi_{\ell} = \lambda_{\ell} \hspace{0.1em} \xi_{\ell}$, then $\varphi_{\ell}({\mathscr{v}_i}) = \xi_{\ell}^\top \phi({\mathscr{v}_i})$ is an eigenfunction of $\mathcal{L}_r$ since

$\begin{align*} \mathcal{L}_r \varphi_{\ell}\left({\mathscr{v}_i}\right) = \left(L_r \hspace{0.1em} \xi_{\ell}\right)^\top \phi\left({\mathscr{v}_i}\right) = \lambda_{\ell} \hspace{0.1em} \xi_{\ell}^\top \phi\left({\mathscr{v}_i}\right) = \lambda_{\ell} \hspace{0.1em} \varphi_{\ell}\left({\mathscr{v}_i}\right). \end{align*}$

Alternatively, assuming the matrix $G^{(0)}$ is invertible, we can also solve the generalized eigenvalue problem $G^{(1)} \hspace{0.1em} \xi_{\ell} = \lambda_{\ell} \hspace{0.1em} G^{(0)} \hspace{0.1em} \xi_{\ell}$. This shows that eigenfunctions of $\mathcal{L}$ can be approximated by eigenfunctions of $\mathcal{L}_r$.

Example 3.1. Choosing r = n and indicator functions for each vertex, i.e.

$\begin{align*} \phi_j\left({\mathscr{v}_i}\right) = \unicode{x1D7D9}_{{\mathscr{v}_j}}\left({\mathscr{v}_i}\right) = \begin{cases} 1, & i = j, \\ 0, & \textrm{otherwise}, \end{cases} \end{align*}$

we obtain $\Phi_{\mathscr{V}} = I$ and hence $G^{(0)} = D_{\mu}$ and $G^{(1)} = D_{\mu} L$ so that $L_r = L$. That is, in this case we obtain the full matrix representation of the operator $\mathcal{L}$. $\triangle$

Theoretically, we could choose any set of basis functions. Our goal, however, is to define the basis functions in such a way that the reduced operator $\mathcal{L}_r$ has essentially the same dominant eigenvalues as the full operator $\mathcal{L}$ and thus retains the metastability (and cluster structure) of the system.

Remark 3.2. The optimal basis is given by the eigenvectors $\boldsymbol{\varphi}_{\ell}$ corresponding to the largest eigenvalues $\lambda_{\ell}$ of L since choosing

$\begin{align*} \Phi_{\mathscr{V}} = \begin{bmatrix} \boldsymbol{\varphi}_1^\top \; \\ \vdots \\ \boldsymbol{\varphi}_r^\top \end{bmatrix} \end{align*}$

results in $G^{(0)} = \Phi_{\mathscr{V}} \hspace{0.1em} D_{\mu} \hspace{0.1em} \Phi_{\mathscr{V}}^\top$ and $G^{(1)} = \Phi_{\mathscr{V}} \hspace{0.1em} D_{\mu} \hspace{0.1em} L \hspace{0.1em} \Phi_{\mathscr{V}}^\top = \Phi_{\mathscr{V}} \hspace{0.1em} D_{\mu} \hspace{0.1em} \Phi_{\mathscr{V}}^\top \operatorname{diag}(\lambda_1, \dots, \lambda_r)$. That is, $L_r = \operatorname{diag}(\lambda_1, \dots, \lambda_r)$ and the reduced operator has exactly the same dominant eigenvalues as the full operator. However, if we already knew the eigenvectors of L, then we could immediately compute the clusters by applying, e.g. k-means. Also, the number of metastable sets is in general not known in advance. The question then is whether it is possible to derive near-optimal basis functions without having to compute the eigenvectors in the first place. One possibility would be to construct basis functions based on random walks started in different parts of the graph since the random walkers would with a high probability first explore nodes within the clusters before moving to other clusters. The generation of suitable basis functions, however, is beyond the scope of this paper.

The operators introduced above can also be estimated from random walk data. Given m random walkers $x^{(i)}$ sampled from the distribution µ, each random walker takes a single step forward according to the probabilities given by the transition probability matrix S. Let the final positions of the random walkers be denoted by $y^{(i)}$. Alternatively, instead of considering many random walkers, the data can be extracted from one long random walk $x^{(1)}, \dots, x^{(m+1)}$ by defining $y^{(i)} = x^{(i+1)}$, for $i = 1, \dots, m$. The density µ is then determined by the random walk. If the graph admits a unique invariant density π, then µ will converge to π for $m \to \infty$.

We now define data matrices $\Phi_x, \Phi_y \in \mathbb{R}^{r \times m}$ by

$\begin{align*} \Phi_x = \begin{bmatrix} \phi\left(x^{\left(1\right)}\right) & \phi\left(x^{\left(2\right)}\right) & \dots & \phi\left(x^{\left(m\right)}\right) \end{bmatrix} \quad \textrm{and} \quad \Phi_y = \begin{bmatrix} \phi\left(y^{\left(1\right)}\right) & \phi\left(y^{\left(2\right)}\right) & \dots & \phi\left(y^{\left(m\right)}\right) \end{bmatrix}. \end{align*}$

Further, let $\widehat{G}_{xx}, \widehat{G}_{yy}, \widehat{G}_{xy} \in \mathbb{R}^{r \times r}$ be defined by

$\begin{align*} \widehat{G}_{xx} &: = \frac{1}{m} \Phi_x \hspace{0.1em} \Phi_x^\top = \frac{1}{m} \sum_{i = 1}^m \phi\left(x^{\left(i\right)}\right) \hspace{0.1em} \phi\left(x^{\left(i\right)}\right)^\top \underset{\scriptscriptstyle m \rightarrow \infty}{\longrightarrow} \sum_{k = 1}^n \mu\left({\mathscr{v}_{k}}\right) \hspace{0.1em} \phi\left({\mathscr{v}_{k}}\right) \hspace{0.1em} \phi\left({\mathscr{v}_{k}}\right)^\top = : G_{xx}, \\ \widehat{G}_{yy} &: = \frac{1}{m} \Phi_y \hspace{0.1em} \Phi_y^\top = \frac{1}{m} \sum_{i = 1}^m \phi\left(y^{\left(i\right)}\right) \hspace{0.1em} \phi\left(y^{\left(i\right)}\right)^\top \underset{\scriptscriptstyle m \rightarrow \infty}{\longrightarrow} \sum_{k = 1}^n \nu\left({\mathscr{v}_{k}}\right) \hspace{0.1em} \phi\left({\mathscr{v}_{k}}\right) \hspace{0.1em} \phi\left({\mathscr{v}_{k}}\right)^\top = : G_{yy}, \\ \widehat{G}_{xy} &: = \frac{1}{m} \Phi_x \hspace{0.1em} \Phi_y^\top = \frac{1}{m} \sum_{i = 1}^m \phi\left(x^{\left(i\right)}\right) \hspace{0.1em} \phi\left(y^{\left(i\right)}\right)^\top \underset{\scriptscriptstyle m \rightarrow \infty}{\longrightarrow} \sum_{k = 1}^n \sum_{l = 1}^n \mu\left({\mathscr{v}_{k}}\right) s_{kl} \hspace{0.1em} \phi\left({\mathscr{v}_{k}}\right) \hspace{0.1em} \phi\left({\mathscr{v}_{l}}\right)^\top = : G_{xy}. \end{align*}$

Thus, the entries of the matrices $G_{xx} = \Phi_{\mathscr{V}} \hspace{0.1em} D_{\mu} \hspace{0.1em} \Phi_{\mathscr{V}}^\top$, $G_{yy} = \Phi_{\mathscr{V}} \hspace{0.1em} D_{\nu} \hspace{0.1em} \Phi_{\mathscr{V}}^\top$, and $G_{xy} = \Phi_{\mathscr{V}} \hspace{0.1em} D_{\mu} \hspace{0.1em} S \hspace{0.1em} \Phi_{\mathscr{V}}^\top$ are given by

$\begin{align*} \left[g_{xx}\right]_{ij} = \left\langle \phi_i,\, \phi_j \right\rangle_{\mu}, \quad \left[g_{yy}\right]_{ij} = \left\langle \phi_i,\, \phi_j \right\rangle_{\nu}, \quad \left[g_{xy}\right]_{ij} = \left\langle \phi_i,\, \mathcal{K} \phi_j \right\rangle_{\mu}. \end{align*}$

Remark 3.3. Note that the matrices G_xx and G_xy are not the Gram matrices used in kernel-based methods. Here, G stands for Galerkin. For functions $f({\mathscr{v}_i}) = c_x^\top \phi({\mathscr{v}_i})$ and $g({\mathscr{v}_i}) = c_y^\top \phi({\mathscr{v}_i})$, it holds that $\left\langle \,f,\, \mathcal{C}_{xx} f\, \right\rangle = c_x^\top \hspace{0.1em} G_{xx} \hspace{0.1em} c_x$ and $\left\langle \,f,\, \mathcal{C}_{xy} g \right\rangle = c_x^\top \hspace{0.1em} G_{xy} \hspace{0.1em} c_y$. If we choose the basis functions defined in example 3.1, then $G_{xx} = C_{xx} = D_{\mu}$ and $G_{xy} = C_{xy} = D_{\mu} \hspace{0.1em} S$.

Lemma 3.4. It follows that:

• (i)  
  $G_{xx}^{-1} \hspace{0.1em} G_{xy}$ is a Galerkin approximation of the operator $\mathcal{K}$ w.r.t. $\left\langle \cdot,\, \cdot \right\rangle_{\mu}$,
• (ii)  
  $G_{yy}^{-1} \hspace{0.1em} G_{yx}$ is a Galerkin approximation of the operator $\mathcal{T}$ w.r.t. $\left\langle \cdot,\, \cdot \right\rangle_{\nu}$.

Proof. For the first part, compare G_xx and G_xy with the matrices $G^{(0)}$ and $G^{(1)}$ introduced above. Then, using proposition A.1, we have

$\begin{align*} \left[g_{xy}\right]_{ij} = \left\langle \mathcal{T} \phi_i,\, \phi_j \right\rangle_{\nu} ~~ \Longrightarrow ~~ \left[g_{yx}\right]_{ij} = \left\langle \phi_i,\, \mathcal{T} \phi_j \right\rangle_{\nu}, \end{align*}$

which proves the second part. □

If µ is the uniform distribution, then $\left\langle \phi_i,\, \mathcal{K} \phi_j \right\rangle = \left\langle \mathcal{P} \phi_i,\, \phi_j \right\rangle$, and $G_{xx}^{-1} \hspace{0.1em} G_{yx}$ is an approximation of the operator $\mathcal{P}$ with respect to the standard inner product $\left\langle \cdot,\, \cdot \right\rangle$. As the forward–backward and backward–forward operators are simply compositions of these operators (see also proposition A.3), all the transfer operators introduced above can be estimated from data. This, however, implies that we do not directly compute Galerkin approximations of $\mathcal{F}$ and $\mathcal{B}$ but represent them as compositions of Galerkin projections, which introduces additional errors.

Example 3.5. For the graph introduced in example 2.16, we define sets $\mathbb{A}_1 = \{{\mathscr{v}_1}, \dots, {\mathscr{v}_{100}} \}$, $\mathbb{A}_2 = \{{\mathscr{v}_{101}}, \dots, {\mathscr{v}_{200}}\}$, $\mathbb{A}_3 = \{{\mathscr{v}_{201}}, \dots, {\mathscr{v}_{300}} \}$, and $\mathbb{A}_4 = \{{\mathscr{v}_{301}}, \dots, {\mathscr{v}_{400}} \}$ and basis functions $\phi_j({\mathscr{v}_i}) = \unicode{x1D7D9}_{\mathbb{A}_j}({\mathscr{v}_i})$, with $j = 1, \dots, 4$, where $\unicode{x1D7D9}_{\mathbb{B}}$ denotes the indicator function for the set $\mathbb{B}$. The representation of the forward–backward operator $\mathcal{F}$ projected onto the space spanned by the four basis functions is then approximately a diagonal matrix and can be regarded as a coarse-grained counterpart of the full operator, where the off-diagonal entries correspond to transition probabilities between the clusters. The resulting eigenvalues and eigenvectors are good approximations of the dominant eigenvalues and eigenvectors of the full operator as shown in figure 5. $\triangle$

Zoom In Zoom Out Reset image size

The example illustrates that it is possible to approximate the dominant eigenvalues and eigenvectors by computing eigendecompositions of potentially much smaller matrices, provided that a suitable set of basis functions is chosen and that the graph exhibits metastability.

Remark 3.6. Assuming that only random-walk data is available, but not the underlying graph itself, the data-driven approach can thus also be used for inferring the network structure and the transition probabilities. This will be further investigated in future work.

We will now show that it is possible to rewrite the clustering approach for directed graphs as a constrained optimization problem. This gives rise to a different interpretation of clusters. Canonical correlation analysis (CCA) [38–40] was originally developed to maximize the correlation between random variables—potentially in infinite-dimensional feature spaces if the kernel-based formulation is used. It has been shown in [9] that CCA, when applied to time-series data, can be used to detect coherent sets. In order to illustrate the relationships with the forward–backward operator $\mathcal{F}$ (or its matrix representation F), we will briefly outline the derivation of CCA, which can be written as

$\begin{align*} \max_{f, g \in \mathbb{U}} \left\langle \,f,\, \mathcal{C}_{xy} \hspace{0.1em} g \right\rangle \quad \textrm{s.t.} \quad \left\langle \,f,\, \mathcal{C}_{xx} \hspace{0.1em} f\, \right\rangle = \left\langle g,\, \mathcal{C}_{yy} \hspace{0.1em} g \right\rangle = 1. \end{align*}$

Restricted to the subspace $\mathbb{U}_r$, defining $f({\mathscr{v}_i}) = c_x^\top \phi({\mathscr{v}_i})$ and $g({\mathscr{v}_i}) = c_y^\top \phi({\mathscr{v}_i})$, this yields

$\begin{align} \max_{c_x, c_y \in \mathbb{R}^r} c_x^\top G_{xy} \hspace{0.1em} c_y \quad \textrm{s.t.} \quad c_x^\top G_{xx} \hspace{0.1em} c_x = c_y^\top G_{yy} \hspace{0.1em} c_y = 1. \end{align} \tag{ 1 }$

The Lagrangian function corresponding to this optimization problem is

$\begin{align*} \mathfrak{L}\left(c_x, c_y, \kappa_x, \kappa_y\right) = c_x^\top G_{xy} \hspace{0.1em} c_y - \frac{\kappa_x}{2} \left(c_x^\top G_{xx} \hspace{0.1em} c_x - 1\right) - \frac{\kappa_y}{2} \left(c_y^\top G_{yy} \hspace{0.1em} c_y - 1\right), \end{align*}$

where κ_x and κ_y are Lagrange multipliers. The derivatives with respect to c_x and c_y are

$\begin{align*} \nabla_{\!c_x} \mathfrak{L}\left(c_x, c_y, \kappa_x, \kappa_y\right) = G_{xy} \hspace{0.1em} c_y - \kappa_x \hspace{0.1em} G_{xx} \hspace{0.1em} c_x \stackrel{!}{=} 0, \\ \nabla_{\!c_y} \mathfrak{L}\left(c_x, c_y, \kappa_x, \kappa_y\right) = G_{yx} \hspace{0.1em} c_x - \kappa_y \hspace{0.1em} G_{yy} \hspace{0.1em} c_y \stackrel{!}{=} 0. \end{align*}$

Multiplying the first equation by $c_x^\top$ and the second by $c_y^\top$, it follows that $\kappa_x = \kappa_y = : \kappa_{\ell}$. We thus obtain the (generalized) eigenvalue problem

$\begin{align*} \begin{bmatrix} 0 & G_{xy} \\ G_{yx} & 0 \end{bmatrix} \begin{bmatrix} c_x \\ c_y \end{bmatrix} = \kappa_{\ell} \begin{bmatrix} G_{xx} & 0 \\ 0 & G_{yy} \end{bmatrix} \begin{bmatrix} c_x \\ c_y \end{bmatrix} \quad \textrm{or} \quad \begin{bmatrix} 0 & G_{xx}^{-1} G_{xy} \\ G_{yy}^{-1} G_{yx} & 0 \end{bmatrix} \begin{bmatrix} c_x \\ c_y \end{bmatrix} = \kappa_{\ell} \begin{bmatrix} c_x \\ c_y \end{bmatrix}. \end{align*}$

If we now solve one of the equations for c_x or c_y and insert the resulting expression into the other, we have

$\begin{align*} G_{xx}^{-1} G_{xy} \hspace{0.1em} G_{yy}^{-1} G_{yx} \hspace{0.1em} c_x = \kappa_{\ell}^2 \hspace{0.1em} c_x \quad \textrm{and} \quad G_{yy}^{-1} G_{yx} \hspace{0.1em} G_{xx}^{-1} G_{xy} \hspace{0.1em} c_y = \kappa_{\ell}^2 \hspace{0.1em} c_y. \end{align*}$

Since $G_{xx}^{-1} G_{xy}$ is a Galerkin approximation of $\mathcal{K}$ and $G_{yy}^{-1} G_{yx}$ a Galerkin approximation of $\mathcal{T}$, the first expression can be used to compute eigenfunctions of $\mathcal{F}$ and the second to compute eigenfunctions of $\mathcal{B}$, where $\lambda_{\ell} = \kappa_{\ell}^2$ is the corresponding eigenvalue. This shows that the functions maximizing the correlation are the eigenfunctions of the forward–backward and backward–forward operators. We can again use the dominant eigenfunctions to compute coherent sets.

Example 3.7. We now generate graphs $\mathscr{G} \in \mathbf{G}\big(2, 100, \big[\begin{smallmatrix} p & q \\ q & p \end{smallmatrix} \big]\big)$ with varying p and q using the DSBM. There are two dense diagonal blocks if p is large and q is small and two dense off-diagonal blocks if p is small and q large. In order to analyze the quality of the clustering for different values of p and q, we plot the second-largest eigenvalue κ₂ (i.e. the correlation between ϕ₂ and ψ₂) and the adjusted Rand index [41]. The results are shown in figure 6. The clustering works for both scenarios, but, as expected, fails when the values of p and q are too similar since there are in that case no clearly defined clusters. $\triangle$

Zoom In Zoom Out Reset image size

This demonstrates that the spectral clustering approach for directed graphs works for a wide range of parameter settings. Alternatively, we can rewrite (1) as

$\begin{align} \max_{\substack{\left\lVert \widetilde{c}_x \right\rVert = 1 \\ \left\lVert \widetilde{c}_y \right\rVert = 1}} \widetilde{c}_x^\top G_{xx}^{{-1/2}} G_{xy} \hspace{0.1em} G_{yy}^{{-1/2}} \hspace{0.1em} \widetilde{c}_y, \end{align} \tag{ 2 }$

where $\widetilde{c}_x = G_{xx}^{{-1/2}} \hspace{0.1em} c_x$ and $\widetilde{c}_y = G_{yy}^{{-1/2}} \hspace{0.1em} c_y$. The maximum of the cost function (2) is the largest singular value σ₁ of the matrix $G_{xx}^{{-1/2}} G_{xy} \hspace{0.1em} G_{yy}^{{-1/2}}$, attained at $\widetilde{c}_x = u_1$ and $\widetilde{c}_y = v_1$, where u₁ and v₁ are the corresponding left and right singular vectors [42]. Thus, setting $c_x = G_{xx}^{{-1/2}} u_1$ and $c_y = G_{yy}^{{-1/2}} v_1$ solves the CCA problem. This can be extended to the subsequent singular values and vectors.

In order to write spectral clustering algorithms for undirected graphs as a solution of a constrained optimization problem, set $\mu = \pi$. We then solve

$\begin{align} \max_{c \in \mathbb{R}^r} c^\top G_{xy} \hspace{0.1em} c \quad \textrm{s.t.} \quad c^\top G_{xx} \hspace{0.1em} c = 1, \end{align} \tag{ 3 }$

resulting in the Lagrangian function

$\begin{align*} \mathfrak{L}\left(c, \kappa\right) = c^\top G_{xy} \hspace{0.1em} c - \kappa \hspace{0.1em} \left(c^\top G_{xx} \hspace{0.1em} c - 1\right) \end{align*}$

and hence

$\begin{align*} \nabla_{\!c} \hspace{0.1em} \mathfrak{L}\left(c, \kappa\right) = 2 \hspace{0.1em} G_{xy} \hspace{0.1em} c - 2 \hspace{0.1em} \kappa \hspace{0.1em} G_{xx} \hspace{0.1em} c \stackrel{!}{=} 0, \end{align*}$

where we used that G_xy is symmetric for $\mu = \pi$ since due to the detailed balance condition $C_{xy} = D_{\pi} \hspace{0.1em} S = S^\top D_{\pi} = C_{xy}^\top$ and hence also $G_{xy} = G_{xy}^\top$. Thus, we obtain the eigenvalue problem $G_{xx}^{-1} G_{xy} \hspace{0.1em} c = \kappa \hspace{0.1em} c$. This is a Galerkin approximation of the Koopman operator and the cost function is maximized by the eigenfunction associated with the largest eigenvalue. Subdominant eigenfunctions can again be used for detecting metastable sets.

Remark 3.8. If the matrix G_xy is not symmetric, we can still solve the optimization problem, but would then obtain the eigenvalue problem $\frac{1}{2} G_{xx}^{-1} (G_{xy} + G_{yx}) \hspace{0.1em} c = \kappa \hspace{0.1em} c$, which can be regarded as a simple form of symmetrization.

Defining $\widetilde{c} = G_{xx}^{{-1/2}} \hspace{0.1em} c$, we can also write (3) as

$\begin{align} \max_{\left\lVert \widetilde{c} \right\rVert = 1} \widetilde{c}^\top G_{xx}^{{-1/2}} G_{xy} \hspace{0.1em} G_{xx}^{{-1/2}} \hspace{0.1em} \widetilde{c}, \end{align} \tag{ 4 }$

which is maximized by the dominant eigenvector of the matrix $G_{xx}^{{-1/2}} G_{xy} \hspace{0.1em} G_{xx}^{{-1/2}}$. We then again obtain the eigenvalue problem $G_{xx}^{-1} G_{xy} \hspace{0.1em} c = \kappa \hspace{0.1em} c$.

We have shown in section 2 that our transfer operator-based spectral clustering algorithm—obtained by replacing the Koopman operator by the forward–backward operator—can be interpreted as a straightforward generalization of conventional spectral clustering algorithms for undirected graphs. This is also corroborated by the optimization problem formulation. Additionally, comparing the constrained optimization problems (1) and (3) for directed and undirected graphs, we see that the former is clearly more flexible and powerful since we are allowed to choose two functions, whereas the latter is restricted to one, which makes it impossible to detect dense off-diagonal blocks.

We generate 300 weighted graphs with heterogeneous node-degree distributions and community sizes as described in [32] using their publicly available C++code (Package 4). All the graphs comprise 500 vertices and have between four and six clusters. We consider three different test cases: (1) only dense blocks on the diagonal, (2) only off-diagonal dense blocks, and (3) a mix of both ⁴ . A typical adjacency matrix (case 3) is shown in figure 7.

Zoom In Zoom Out Reset image size

We apply the degree-discounted bibliometric symmetrization [22], the Hermitian clustering method proposed in [31] (both described in remark 2.13), and algorithm 2.15 to the three different test cases. We run each algorithm ten times (using the same k-means implementation and settings) for each graph and compute the average adjusted Rand index and the average number of misclassified vertices. The results, listed in table 1, show that the Hermitian clustering approach works well for graphs with dense off-diagonal blocks, but fails to identify 'conventional' clusters, i.e. groups of nodes that are strongly coupled by directed edges, since the method was designed to detect only imbalances in the orientation of the edges. Instead of the random-walk normalized matrix $A_{\textrm{rw}} = D_{\mathscr{o}}^{-1} A$ as proposed in [31], we use the normalized adjacency matrix $A_{\textrm{nn}} = D_{\mathscr{o}}^{{-1/2}} A D_{\mathscr{i}}^{{-1/2}}$, which seems to improve the results. Although the reweighting parameters for the degree-discounted bibliometric symmetrization were determined empirically, the method works for all test cases. The transfer operator-based clustering also works well for all considered graph types and leads to marginally better results. The advantage of the dynamical systems approach is that it offers a rigorous interpretation of clusters in terms of coherent sets (or CCA) and that it can be easily extended to dynamic graphs. For more complicated benchmark problems, it might also be beneficial to tune the probability density µ or to use a combination of the operators $\mathcal{F}$ and $\mathcal{B}$.

As a last example, we consider a 32-bit adder. The data set can be found on the Matrix Market website. Since the matrix contains negative values, we shift the entries so that all weights are positive. The circuit has a fairly simple structure and the spectral clustering algorithm detects 32 blocks of nearly the same size. Reordering the adjacency matrix based on the cluster numbers, we obtain a block matrix with very few nonzero entries in the off-diagonal blocks as shown in figure 8. The bandwidth of the matrix could be minimized by applying generalizations of the Cuthill–McKee algorithm to the block structure [43].

Zoom In Zoom Out Reset image size